Fix InputOutput when enable_bubble is true

It turns out that we were not saving and restoring a piece of
information that is required to fully reconstruct a `2DSpectralSpace`:
the bubble correction.

When enabled, the bubble correction recomputes the Jacobian to improve
their accuracy. `InputOutput` was not saving this piece of information
into the HDF5 file, so spaces read from file always assumed that it was
false. This led to discrepancies in the space, and these discrepancies
led to problems in the evolution.

Interestingly, for the low-resolution cases we consider in the Atmos CI,
`enable_bubble` has very small effects on the value of the Jacobian,
with a relative difference compared to the case without bubble smaller
than the floating point precision for Float32. In spite of this, it
affected the evolution in more noticeable way, in particular in vectors
and with hyperdiffusion.

Now, `bubble` is saved as attribute so that it can be read back.

This commit also bumps the version of ClimaCore to 0.14.6

NEWS.md

@@ -4,6 +4,22 @@ ClimaCore.jl Release Notes
 main
 -------
 
+v0.14.16
+-------
+
+### ![][badge-🐛bugfix] Fix restarting simulations from `Space`s with `enable_bubble = true`
+
+Prior to this change, the `ClimaCore.InputOutput` module did not save whether a
+`Space` was constructed with `enable_bubble = true`. This meant that restarting
+a simulation from a HDF5 file led to inconsistent and incorrect spaces and
+`Field`s. This affected only 2D spectral spaces (and extruded ones that have
+this type of horizontal space).
+
+We now expect `Space`s read from a file to be bitwise identical to the original
+one.
+
+PR [#1999](https://github.com/CliMA/ClimaCore.jl/pull/1999).
+
 v0.14.15
 -------
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ClimaCore"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
 authors = ["CliMA Contributors <[email protected]>"]
-version = "0.14.15"
+version = "0.14.16"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

src/Grids/spectralelement.jl

@@ -129,6 +129,7 @@ mutable struct SpectralElementGrid2D{
     local_dss_weights::D
     internal_surface_geometry::IS
     boundary_surface_geometries::BS
+    enable_bubble::Bool
 end
 
 local_geometry_type(
@@ -467,6 +468,7 @@ function _SpectralElementGrid2D(
         dss_local_weights,
         internal_surface_geometry,
         boundary_surface_geometries,
+        enable_bubble,
     )
 end
 

src/InputOutput/readers.jl

@@ -333,10 +333,15 @@ function read_grid_new(reader, name)
         quadrature_style =
             _scan_quadrature_style(attrs(group)["quadrature_type"], npts)
         topology = read_topology(reader, attrs(group)["topology"])
+        enable_bubble = get(attrs(group), "bubble", "false") == "true"
         if type == "SpectralElementGrid1D"
             return Grids.SpectralElementGrid1D(topology, quadrature_style)
         else
-            return Grids.SpectralElementGrid2D(topology, quadrature_style)
+            return Grids.SpectralElementGrid2D(
+                topology,
+                quadrature_style;
+                enable_bubble,
+            )
         end
     elseif type == "FiniteDifferenceGrid"
         topology = read_topology(reader, attrs(group)["topology"])

src/InputOutput/writers.jl

@@ -343,6 +343,7 @@ function write_new!(
         "quadrature_num_points",
         Quadratures.degrees_of_freedom(Spaces.quadrature_style(grid)),
     )
+    write_attribute(group, "bubble", grid.enable_bubble ? "true" : "false")
     write_attribute(group, "topology", write!(writer, Spaces.topology(grid)))
     return name
 end

test/InputOutput/spectralelement2d.jl

@@ -71,20 +71,27 @@ end
     @info "Using device" device
     context = ClimaComms.SingletonCommsContext(device)
     grid_topology = Topologies.Topology2D(context, mesh)
-    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
 
-    y0 = init_state.(Fields.local_geometry_field(space), Ref(parameters))
-    Y = Fields.FieldVector(y0 = y0)
+    for enable_bubble in (true, false)
+        space =
+            Spaces.SpectralElementSpace2D(grid_topology, quad; enable_bubble)
 
-    # write field vector to hdf5 file
-    filename = tempname(pwd())
-    writer = InputOutput.HDF5Writer(filename, context)
-    InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
-    close(writer)
-    reader = InputOutput.HDF5Reader(filename, context)
-    restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
-    close(reader)
-    ClimaComms.allowscalar(device) do
-        @test restart_Y == Y # test if restart is exact
+        y0 = init_state.(Fields.local_geometry_field(space), Ref(parameters))
+        Y = Fields.FieldVector(y0 = y0)
+
+        # write field vector to hdf5 file
+        filename = tempname(pwd())
+        writer = InputOutput.HDF5Writer(filename, context)
+        InputOutput.write!(writer, "Y" => Y) # write field vector from hdf5 file
+        close(writer)
+        reader = InputOutput.HDF5Reader(filename, context)
+        restart_Y = InputOutput.read_field(reader, "Y") # read fieldvector from hdf5 file
+        close(reader)
+        ClimaComms.allowscalar(device) do
+            @test restart_Y == Y # test if restart is exact
+        end
+        ClimaComms.allowscalar(device) do
+            @test axes(restart_Y) == axes(Y) # test if restart is exact for space
+        end
     end
 end