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

@@ -13,6 +13,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
 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

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"])

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
 

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