Add default_num_points

`default_num_points` is a new function that automatically computes a
tuple with the number of points along each direction of a `space`.
`default_num_points` is now the default value for the `NetCDF` writer,
so that users will automatically get a reasonably sampled output.

Author: Sbozzolo

NEWS.md

@@ -14,6 +14,14 @@ the interval `[0.0, 10.0]`. If one knows that the data represents a time
 average, then the time of `10.0` is the time average over the interval
 `[0.0, 10.0]`.
 
+### NetCDF writer now defaults to a reasonable number of points
+
+`ClimaDiagnostics.Writers.NetCDF` now has a new default argument that depends on
+the input Space. With this new default, obtained by calling the
+`ClimaDiagnostics.Writers.default_num_points(space)` function, the output
+diagnostics will be sampled with approximately the same resolution as the given
+`space`.
+
 ### Support for `lazy`
 
 Starting version `0.2.13`, `ClimaDiagnostics` supports diagnostic variables

docs/src/api.md

@@ -50,5 +50,6 @@ ClimaDiagnostics.Writers.NetCDFWriter
 ClimaDiagnostics.Writers.HDF5Writer
 ClimaDiagnostics.Writers.interpolate_field!
 ClimaDiagnostics.Writers.write_field!
+ClimaDiagnostics.Writers.default_num_points
 Base.close
 ```

src/Writers.jl

@@ -15,6 +15,8 @@ Writers can implement:
 """
 module Writers
 
+import ClimaCore
+
 import ..AbstractWriter, ..ScheduledDiagnostic
 import ..ScheduledDiagnostics: output_short_name, output_long_name
 

src/netcdf_writer.jl

@@ -121,7 +121,7 @@ Keyword arguments
 function NetCDFWriter(
     space::Spaces.AbstractSpace,
     output_dir;
-    num_points = (180, 90, 50),
+    num_points = default_num_points(space),
     compression_level = 9,
     sync_schedule = ClimaComms.device(space) isa ClimaComms.CUDADevice ?
                     EveryStepSchedule() : nothing,
@@ -149,7 +149,7 @@ function NetCDFWriter(
                 Spaces.horizontal_space(space) isa
                 Spaces.SpectralElementSpace1D ? 1 : 2
 
-            num_vpts = Meshes.nelements(Grids.vertical_topology(space).mesh)
+            num_vpts = Spaces.nlevels(space)
 
             # For any configuration, it is reasonable to assume that the last
             # value of `num_pts` is the number of vertical points
@@ -214,16 +214,17 @@ end
 function NetCDFWriter(
     space::Spaces.Spaces.FiniteDifferenceSpace,
     output_dir;
-    num_points = (180, 90, 50),
+    num_points = default_num_points(space),
     compression_level = 9,
     sync_schedule = ClimaComms.device(space) isa ClimaComms.CUDADevice ?
                     EveryStepSchedule() : nothing,
     z_sampling_method = LevelsMethod(),
     start_date = nothing,
 )
     if z_sampling_method isa LevelsMethod
-        num_vpts = Meshes.nelements(Grids.vertical_topology(space).mesh)
-        @warn "Disabling vertical interpolation, the provided number of points is ignored (using $num_vpts)"
+        num_vpts = Spaces.nlevels(ClimaCore.Spaces.center_space(space))
+        num_vpts == last(num_points) ||
+            @warn "Disabling vertical interpolation, the provided number of points is ignored (using $num_vpts)"
         num_points = (num_vpts,)
     end
     vpts = target_coordinates(space, num_points, z_sampling_method)

src/netcdf_writer_coordinates.jl

@@ -661,3 +661,90 @@ Prepare the matrix of horizontal coordinates with the correct type according to
 and `domain` (e.g., `ClimaCore.Geometry.LatLongPoint`s).
 """
 function hcoords_from_horizontal_space(space, domain, hpts) end
+
+
+"""
+    default_num_points(space)
+
+Return a tuple with number of points that are optimally suited to interpolate the given
+`space`.
+
+"Optimally suited" here means approximately the same as the number of points as the given
+`space`.
+"""
+function default_num_points(
+    space::ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace,
+)
+    horizontal_space = ClimaCore.Spaces.horizontal_space(space)
+    num_horz = default_num_points(horizontal_space)
+
+    vertical_space = ClimaCore.Spaces.FiniteDifferenceSpace(
+        Spaces.vertical_topology(space),
+        Spaces.staggering(space),
+    )
+    num_vert = default_num_points(vertical_space)
+    return (num_horz..., num_vert...)
+end
+
+# 2D sphere
+function default_num_points(
+    space::ClimaCore.Spaces.CubedSphereSpectralElementSpace2D,
+)
+    # A cubed sphere has 4 panels to cover the range of longitudes, each panel has
+    # `num_elements_per_panel` elements, each with `unique_degrees_of_freedom` points. Same
+    # for latitudes, except that we need 2 panels to cover from 0 to 180.
+
+    unique_degrees_of_freedom = ClimaCore.Quadratures.unique_degrees_of_freedom(
+        ClimaCore.Grids.quadrature_style(space),
+    )
+    num_elements_per_panel = ClimaCore.Meshes.n_elements_per_panel_direction(
+        ClimaCore.Spaces.topology(space).mesh,
+    )
+    num_lat = 2 * num_elements_per_panel * unique_degrees_of_freedom
+    num_lon = 2num_lat
+    return (num_lon, num_lat)
+end
+
+# TODO: Maybe move to ClimaCore?
+const RectilinearSpectralElementSpace1D =
+    ClimaCore.Spaces.SpectralElementSpace1D{
+        <:ClimaCore.Grids.SpectralElementGrid1D{
+            <:ClimaCore.Topologies.IntervalTopology,
+        },
+    }
+
+# 1D box
+function default_num_points(space::RectilinearSpectralElementSpace1D)
+    unique_degrees_of_freedom = ClimaCore.Quadratures.unique_degrees_of_freedom(
+        ClimaCore.Grids.quadrature_style(space),
+    )
+    return (
+        unique_degrees_of_freedom *
+        ClimaCore.Meshes.nelements(ClimaCore.Spaces.topology(space).mesh),
+    )
+end
+
+# 2D box
+function default_num_points(
+    space::ClimaCore.Spaces.RectilinearSpectralElementSpace2D,
+)
+    unique_degrees_of_freedom = ClimaCore.Quadratures.unique_degrees_of_freedom(
+        ClimaCore.Grids.quadrature_style(space),
+    )
+
+    return (
+        unique_degrees_of_freedom * ClimaCore.Meshes.nelements(
+            ClimaCore.Spaces.topology(space).mesh.intervalmesh1,
+        ),
+        unique_degrees_of_freedom * ClimaCore.Meshes.nelements(
+            ClimaCore.Spaces.topology(space).mesh.intervalmesh2,
+        ),
+    )
+end
+
+# Column
+function default_num_points(space::Spaces.FiniteDifferenceSpace)
+    # We always want the center space for interpolation
+    cspace = Spaces.center_space(space)
+    return (ClimaCore.Spaces.nlevels(cspace),)
+end

test/writers.jl

@@ -10,6 +10,7 @@ import ClimaCore
 import ClimaCore.Fields
 import ClimaCore.Spaces
 import ClimaCore.Geometry
+import ClimaCore.CommonSpaces
 import ClimaComms
 
 import ClimaDiagnostics
@@ -38,6 +39,51 @@ output_dir = mktempdir(pwd())
 end
 
 @testset "NetCDFWriter" begin
+    @testset "default_num_points" begin
+        @test Writers.default_num_points(
+            CommonSpaces.ExtrudedCubedSphereSpace(;
+                z_elem = 10,
+                z_min = 0,
+                z_max = 1,
+                radius = 10,
+                h_elem = 10,
+                n_quad_points = 4,
+                staggering = CommonSpaces.CellCenter(),
+            ),
+        ) == (120, 60, 10)
+
+        @test Writers.default_num_points(
+            CommonSpaces.SliceXZSpace(;
+                z_elem = 10,
+                x_min = 0,
+                x_max = 1,
+                z_min = 0,
+                z_max = 1,
+                periodic_x = false,
+                n_quad_points = 4,
+                x_elem = 4,
+                staggering = CommonSpaces.CellCenter(),
+            ),
+        ) == (12, 10)
+        @test Writers.default_num_points(
+            CommonSpaces.Box3DSpace(;
+                z_elem = 10,
+                x_min = 0,
+                x_max = 1,
+                y_min = 0,
+                y_max = 1,
+                z_min = 0,
+                z_max = 10,
+                periodic_x = false,
+                periodic_y = false,
+                n_quad_points = 4,
+                x_elem = 3,
+                y_elem = 4,
+                staggering = CommonSpaces.CellCenter(),
+            ),
+        ) == (9, 12, 10)
+    end
+
     space = SphericalShellSpace()
     field = Fields.coordinate_field(space).z