Add information about deep grid in restart files

NEWS.md

@@ -12,6 +12,25 @@ v0.14.17
  - Fixed some type instabilities PR [#2004](https://github.com/CliMA/ClimaCore.jl/pull/2004)
  - More fixes to higher resolution column cases for the GPU [#1854](https://github.com/CliMA/ClimaCore.jl/pull/1854)
 
+### ![][badge-🐛bugfix] Fix equality for `FieldVector`s with different type
+
+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
 -------
 

src/Fields/fieldvector.jl

@@ -487,10 +487,18 @@ end
 
 Recursively compare given fieldvectors via `==`.
 Returns `true` if `x == y` recursively.
+
+FieldVectors with different types are considered different.
 """
 rcompare(x::T, y::T) where {T <: Union{FieldVector, NamedTuple}} =
     _rcompare(true, x, y)
 
-# Define == to call rcompare for two fieldvectors of the same
-# exact type.
-Base.:(==)(x::T, y::T) where {T <: FieldVector} = rcompare(x, y)
+rcompare(x::T, y::T) where {T <: FieldVector} = _rcompare(true, x, y)
+
+rcompare(x::T, y::T) where {T <: NamedTuple} = _rcompare(true, x, y)
+
+# FieldVectors with different types are always different
+rcompare(x::FieldVector, y::FieldVector) = false
+
+# Define == to call rcompare for two fieldvectors
+Base.:(==)(x::FieldVector, y::FieldVector) = rcompare(x, y)

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