From bd8a1b832bff8330878d20298aed8cf418d0f90e Mon Sep 17 00:00:00 2001 From: Gabriele Bozzola Date: Wed, 2 Oct 2024 08:30:42 -0700 Subject: [PATCH] Add information about deep grid in restart files This commit adds some information that was lost in saving to a checkpoint: the `deep` option. It also adds tests for types too, making the tests more stringent --- NEWS.md | 11 ++ src/InputOutput/readers.jl | 4 +- src/InputOutput/writers.jl | 5 + .../hybrid3dcubedsphere_topography.jl | 111 ++++++++++-------- 4 files changed, 79 insertions(+), 52 deletions(-) diff --git a/NEWS.md b/NEWS.md index a6f3d88395..49972a122e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -18,6 +18,17 @@ Due to a bug, `==` was not recursively checking `FieldVector`s with different types, which resulted in false positives. This is now fixed and `FieldVector`s with different types are always considered different. +### ![][badge-🐛bugfix] Fix restarting simulations from `Space`s with `deep = true` + +Prior to this change, the `ClimaCore.InputOutput` module did not save whether a +`Space` was constructed with `deep = true`. This meant that restarting a +simulation from a HDF5 file led to inconsistent and incorrect spaces and +`Field`s. This affected only extruded 3D spectral spaces. + +We now expect `Space`s read from a file to be bitwise identical to the original +one. + + PR [#2021](https://github.com/CliMA/ClimaCore.jl/pull/2021). v0.14.16 diff --git a/src/InputOutput/readers.jl b/src/InputOutput/readers.jl index 6affc65d43..2fa02f09e1 100644 --- a/src/InputOutput/readers.jl +++ b/src/InputOutput/readers.jl @@ -351,6 +351,7 @@ function read_grid_new(reader, name) vertical_grid = read_grid(reader, attrs(group)["vertical_grid"]) horizontal_grid = read_grid(reader, attrs(group)["horizontal_grid"]) hypsography_type = get(attrs(group), "hypsography_type", "Flat") + deep = get(attrs(group), "deep", false) if hypsography_type == "Flat" hypsography = Grids.Flat() elseif hypsography_type == "LinearAdaption" @@ -370,7 +371,8 @@ function read_grid_new(reader, name) return Grids.ExtrudedFiniteDifferenceGrid( horizontal_grid, vertical_grid, - hypsography, + hypsography; + deep, ) elseif type == "LevelGrid" full_grid = read_grid(reader, attrs(group)["full_grid"]) diff --git a/src/InputOutput/writers.jl b/src/InputOutput/writers.jl index 493c82f340..68576886d5 100644 --- a/src/InputOutput/writers.jl +++ b/src/InputOutput/writers.jl @@ -396,6 +396,11 @@ function write_new!( write!(writer, hypsography.surface, "_z_surface/$name"), ) end + write_attribute( + group, + "deep", + grid.global_geometry isa Geometry.DeepSphericalGlobalGeometry, + ) return name end diff --git a/test/InputOutput/hybrid3dcubedsphere_topography.jl b/test/InputOutput/hybrid3dcubedsphere_topography.jl index fece7d85d5..3cc8301a54 100644 --- a/test/InputOutput/hybrid3dcubedsphere_topography.jl +++ b/test/InputOutput/hybrid3dcubedsphere_topography.jl @@ -10,7 +10,8 @@ using ClimaCore: Fields, DataLayouts, Hypsography, - InputOutput + InputOutput, + Grids using ClimaComms const comms_ctx = ClimaComms.context(ClimaComms.CPUSingleThreaded()) @@ -20,64 +21,72 @@ if ClimaComms.iamroot(comms_ctx) @info "Comms context" comms_ctx nprocs filename end -@testset "HDF5 restart test for 3d hybrid cubed sphere" begin - FT = Float32 - R = FT(6.371229e6) +@testset "HDF5 restart test for 3d hybrid deep cubed sphere with topography and deep" begin + for deep in (false, true) + FT = Float32 + R = FT(6.371229e6) - npoly = 4 - z_max = FT(30e3) - z_elem = 10 - h_elem = 4 - device = ClimaComms.device(comms_ctx) - # horizontal space - domain = Domains.SphereDomain(R) - horizontal_mesh = Meshes.EquiangularCubedSphere(domain, h_elem) - topology = Topologies.Topology2D( - comms_ctx, - horizontal_mesh, - Topologies.spacefillingcurve(horizontal_mesh), - ) - quad = Quadratures.GLL{npoly + 1}() - h_space = Spaces.SpectralElementSpace2D(topology, quad) - # vertical space - z_domain = Domains.IntervalDomain( - Geometry.ZPoint(zero(z_max)), - Geometry.ZPoint(z_max); - boundary_names = (:bottom, :top), - ) + npoly = 4 + z_max = FT(30e3) + z_elem = 10 + h_elem = 4 + device = ClimaComms.device(comms_ctx) + # horizontal space + domain = Domains.SphereDomain(R) + horizontal_mesh = Meshes.EquiangularCubedSphere(domain, h_elem) + topology = Topologies.Topology2D( + comms_ctx, + horizontal_mesh, + Topologies.spacefillingcurve(horizontal_mesh), + ) + quad = Quadratures.GLL{npoly + 1}() + h_space = Spaces.SpectralElementSpace2D(topology, quad) + # vertical space + z_domain = Domains.IntervalDomain( + Geometry.ZPoint(zero(z_max)), + Geometry.ZPoint(z_max); + boundary_names = (:bottom, :top), + ) - z_surface = - Geometry.ZPoint.( - z_max / 8 .* ( - cosd.(Fields.coordinate_field(h_space).lat) .+ - cosd.(Fields.coordinate_field(h_space).long) .+ 1 + z_surface = + Geometry.ZPoint.( + z_max / 8 .* ( + cosd.(Fields.coordinate_field(h_space).lat) .+ + cosd.(Fields.coordinate_field(h_space).long) .+ 1 + ) ) + + z_mesh = Meshes.IntervalMesh(z_domain, nelems = z_elem) + z_topology = Topologies.IntervalTopology(device, z_mesh) + z_space = Spaces.CenterFiniteDifferenceSpace(z_topology) + # Extruded 3D space + h_grid = Spaces.grid(h_space) + z_grid = Spaces.grid(z_space) + hypsography = Hypsography.LinearAdaption(z_surface) + grid = Grids.ExtrudedFiniteDifferenceGrid( + h_grid, + z_grid, + hypsography; + deep, ) - z_mesh = Meshes.IntervalMesh(z_domain, nelems = z_elem) - z_topology = Topologies.IntervalTopology(device, z_mesh) - z_space = Spaces.CenterFiniteDifferenceSpace(z_topology) - # Extruded 3D space - center_space = Spaces.ExtrudedFiniteDifferenceSpace( - h_space, - z_space, - Hypsography.LinearAdaption(z_surface), - ) - face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(center_space) + center_space = Spaces.CenterExtrudedFiniteDifferenceSpace(grid) + face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(grid) - ᶜlocal_geometry = Fields.local_geometry_field(center_space) - ᶠlocal_geometry = Fields.local_geometry_field(face_space) + ᶜlocal_geometry = Fields.local_geometry_field(center_space) + ᶠlocal_geometry = Fields.local_geometry_field(face_space) - Y = Fields.FieldVector(c = ᶜlocal_geometry, f = ᶠlocal_geometry) + Y = Fields.FieldVector(c = ᶜlocal_geometry, f = ᶠlocal_geometry) - # write field vector to hdf5 file - writer = InputOutput.HDF5Writer(filename, comms_ctx) - InputOutput.write!(writer, Y, "Y") - close(writer) + # write field vector to hdf5 file + writer = InputOutput.HDF5Writer(filename, comms_ctx) + InputOutput.write!(writer, Y, "Y") + close(writer) - reader = InputOutput.HDF5Reader(filename, comms_ctx) - restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file - close(reader) - @test restart_Y == Y # test if restart is exact + reader = InputOutput.HDF5Reader(filename, comms_ctx) + restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file + close(reader) + @test restart_Y == Y # test if restart is exact + end end