address most comments

Author: nefrathenrici

.buildkite/PrecompileCI/src/PrecompileCI.jl

@@ -25,21 +25,26 @@ import ClimaParams
 
         sphere_grid = CA.SphereGrid(
             FT;
-            context, radius, h_elem, nh_poly, z_elem, z_max, z_stretch,
-            dz_bottom, bubble,
-            topography,
+            context,
+            radius, h_elem, nh_poly,
+            z_elem, z_max, z_stretch, dz_bottom,
+            bubble, topography,
         )
         box_grid = CA.BoxGrid(
-            FT; context, x_elem, x_max, y_elem, y_max, z_elem, z_max, nh_poly,
-            z_stretch, dz_bottom, bubble,
-            periodic_x = true, periodic_y = true, topography,
+            FT;
+            context,
+            x_elem, x_max, y_elem, y_max, nh_poly, periodic_x = true, periodic_y = true,
+            z_elem, z_max, z_stretch, dz_bottom,
+            bubble, topography,
         )
         plane_grid = CA.PlaneGrid(
-            FT; context, x_elem, x_max, z_elem, z_max, nh_poly, z_stretch,
-            dz_bottom, bubble, periodic_x = true,
-            topography,
+            FT;
+            context,
+            x_elem, x_max, nh_poly, periodic_x = true,
+            z_elem, z_max, z_stretch, dz_bottom,
+            bubble, topography,
         )
-        column_grid = CA.ColGrid(
+        column_grid = CA.ColumnGrid(
             FT; context, z_elem, z_max, z_stretch, dz_bottom,
         )
         all_grids = (sphere_grid, box_grid, plane_grid, column_grid)

docs/src/api.md

@@ -52,7 +52,7 @@ ClimaAtmos.InitialConditions.ColumnInterpolatableField
 ## Grids
 
 ```@docs
-ClimaAtmos.ColGrid
+ClimaAtmos.ColumnGrid
 ClimaAtmos.SphereGrid
 ClimaAtmos.PlaneGrid
 ClimaAtmos.BoxGrid

examples/topography_spectra.jl

@@ -14,11 +14,57 @@ import ClimaCoreSpectra: power_spectrum_2d
 
 const AA = AtmosArtifacts
 
+using ClimaCore:
+    Geometry, Domains, Meshes, Topologies, Spaces, Grids, Hypsography, Fields
+using ClimaComms
+using ClimaUtilities: SpaceVaryingInputs.SpaceVaryingInput
 
-function mask(x::FT) where {FT}
-    return x * FT(x > 0)
+# h_elem is the number of elements per side of every panel (6 panels in total)
+function cubed_sphere_mesh(; radius, h_elem)
+    domain = Domains.SphereDomain(radius)
+    return Meshes.EquiangularCubedSphere(domain, h_elem)
 end
 
+function make_horizontal_space(
+    mesh,
+    quad,
+    comms_ctx::ClimaComms.SingletonCommsContext,
+    bubble,
+)
+
+    space = if mesh isa Meshes.AbstractMesh1D
+        topology = Topologies.IntervalTopology(comms_ctx, mesh)
+        Spaces.SpectralElementSpace1D(topology, quad)
+    elseif mesh isa Meshes.AbstractMesh2D
+        topology = Topologies.Topology2D(
+            comms_ctx,
+            mesh,
+            Topologies.spacefillingcurve(mesh),
+        )
+        Spaces.SpectralElementSpace2D(topology, quad; enable_bubble = bubble)
+    end
+    return space
+end
+
+function make_horizontal_space(mesh, quad, comms_ctx, bubble)
+    if mesh isa Meshes.AbstractMesh1D
+        error("Distributed mode does not work with 1D horizontal spaces.")
+    elseif mesh isa Meshes.AbstractMesh2D
+        topology = Topologies.DistributedTopology2D(
+            comms_ctx,
+            mesh,
+            Topologies.spacefillingcurve(mesh),
+        )
+        space = Spaces.SpectralElementSpace2D(
+            topology,
+            quad;
+            enable_bubble = bubble,
+        )
+    end
+    return space
+end
+
+
 """
     generate_spaces(; h_elem=16)
 For a given number of elements `h_elem`, this function generates a 
@@ -37,10 +83,10 @@ function generate_spaces(;
 )
     FT = Float32
     cubed_sphere_mesh =
-        CA.cubed_sphere_mesh(; radius = FT(planet_radius), h_elem)
+        cubed_sphere_mesh(; radius = FT(planet_radius), h_elem)
     quad = Quadratures.GLL{4}()
     comms_ctx = ClimaComms.context()
-    h_space = CA.make_horizontal_space(cubed_sphere_mesh, quad, comms_ctx, true)
+    h_space = make_horizontal_space(cubed_sphere_mesh, quad, comms_ctx, true)
     Δh_scale = Spaces.node_horizontal_length_scale(h_space)
     @assert h_space isa CC.Spaces.SpectralElementSpace2D
     coords = CC.Fields.coordinate_field(h_space)
@@ -49,7 +95,7 @@ function generate_spaces(;
         "z",
         h_space,
     )
-    elev_from_file = @. mask(elev_from_file)
+    elev_from_file = @. max(elev_from_file, 0)
     # Fractional damping of smallest resolved scale
     # Approximated as k₀ ≈ 1/Δx, with n_attenuation 
     # the factor by which we wish to damp wavenumber
@@ -58,7 +104,7 @@ function generate_spaces(;
     κ = diff_courant * Δh_scale^2
     maxiter = Int(round(log(n_attenuation) / diff_courant))
     diffuse_surface_elevation!(elev_from_file; κ, dt = FT(1), maxiter)
-    elev_from_file = @. mask(elev_from_file)
+    elev_from_file = @. max(elev_from_file, 0)
     return elev_from_file
 end
 

post_processing/ci_plots.jl

@@ -1494,112 +1494,24 @@ function make_plots(
     sim_type::Union{
         EDMFBoxPlots,
         DiagEDMFBoxPlotsWithPrecip,
-    },
-    output_paths::Vector{<:AbstractString},
-)
-    simdirs = SimDir.(output_paths)
-
-    if sim_type isa DiagEDMFBoxPlotsWithPrecip
-        precip_names = ("husra", "hussn", "husraup", "hussnup")
-    else
-        precip_names = ()
-    end
-
-    short_names = [
-        "wa",
-        "waup",
-        "ta",
-        "taup",
-        "hus",
-        "husup",
-        "arup",
-        "tke",
-        "ua",
-        "thetaa",
-        "thetaaup",
-        "ha",
-        "haup",
-        "hur",
-        "hurup",
-        "lmix",
-        "cl",
-        "clw",
-        "clwup",
-        "cli",
-        "cliup",
-        precip_names...,
-    ]
-    reduction = "inst"
-
-    available_periods = ClimaAnalysis.available_periods(
-        simdirs[1];
-        short_name = short_names[1],
-        reduction,
-    )
-    # choose the shortest available period
-    available_periods = collect(available_periods) # ensure vector for indexing
-    period = available_periods[argmin(CA.time_to_seconds.(available_periods))]
-
-    short_name_tuples = pair_edmf_names(short_names)
-    var_groups_zt =
-        map_comparison(simdirs, short_name_tuples) do simdir, name_tuple
-            return [
-                slice(
-                    get(simdir; short_name, reduction, period),
-                    x = 0.0,
-                    y = 0.0,
-                ) for short_name in name_tuple
-            ]
-        end
-
-    var_groups_z = [
-        ([slice(v, time = LAST_SNAP) for v in group]...,) for
-        group in var_groups_zt
-    ]
-
-    tmp_file = make_plots_generic(
-        output_paths,
-        output_name = "tmp",
-        var_groups_z;
-        plot_fn = plot_edmf_vert_profile!,
-        MAX_NUM_COLS = 2,
-        MAX_NUM_ROWS = 4,
-    )
-
-    make_plots_generic(
-        output_paths,
-        vcat(var_groups_zt...),
-        plot_fn = plot_parsed_attribute_title!,
-        summary_files = [tmp_file],
-        MAX_NUM_COLS = 2,
-        MAX_NUM_ROWS = 4,
-    )
-end
-
-function make_plots(
-    sim_type::Union{
         EDMFColumnPlots,
         EDMFColumnPlotsWithPrecip,
     },
     output_paths::Vector{<:AbstractString},
 )
     simdirs = SimDir.(output_paths)
 
-    if sim_type isa EDMFColumnPlotsWithPrecip
+    # Determine if this is a box or column type
+    is_box = sim_type isa Union{EDMFBoxPlots, DiagEDMFBoxPlotsWithPrecip}
+
+    # Determine precipitation names based on type
+    if sim_type isa DiagEDMFBoxPlotsWithPrecip
+        precip_names = ("husra", "hussn", "husraup", "hussnup")
+    elseif sim_type isa EDMFColumnPlotsWithPrecip
         if sim_type isa Val{:prognostic_edmfx_rico_column_2M}
             precip_names = (
-                "husra",
-                "hussn",
-                "husraup",
-                "hussnup",
-                "husraen",
-                "hussnen",
-                "cdnc",
-                "ncra",
-                "cdncup",
-                "ncraup",
-                "cdncen",
-                "ncraen",
+                "husra", "hussn", "husraup", "hussnup", "husraen", "hussnen",
+                "cdnc", "ncra", "cdncup", "ncraup", "cdncen", "ncraen",
             )
         else
             precip_names =
@@ -1610,27 +1522,9 @@ function make_plots(
     end
 
     short_names = [
-        "wa",
-        "waup",
-        "ta",
-        "taup",
-        "hus",
-        "husup",
-        "arup",
-        "tke",
-        "ua",
-        "thetaa",
-        "thetaaup",
-        "ha",
-        "haup",
-        "hur",
-        "hurup",
-        "lmix",
-        "cl",
-        "clw",
-        "clwup",
-        "cli",
-        "cliup",
+        "wa", "waup", "ta", "taup", "hus", "husup", "arup", "tke", "ua",
+        "thetaa", "thetaaup", "ha", "haup", "hur", "hurup", "lmix",
+        "cl", "clw", "clwup", "cli", "cliup",
         precip_names...,
     ]
     reduction = "inst"
@@ -1647,9 +1541,13 @@ function make_plots(
     short_name_tuples = pair_edmf_names(short_names)
     var_groups_zt =
         map_comparison(simdirs, short_name_tuples) do simdir, name_tuple
-            return [
-                get(simdir; short_name, reduction, period) for short_name in name_tuple
-            ]
+            vars = map(short_name -> get(simdir; short_name, reduction, period), name_tuple)
+            # For box types, slice to a point (x=0, y=0)
+            if is_box
+                return map(var -> slice(var, x = 0.0, y = 0.0), vars)
+            else
+                return vars
+            end
         end
 
     var_groups_z = [

reproducibility_tests/ref_counter.jl

@@ -4,7 +4,7 @@
 #
 # What is the ref_counter?
 #
-# The ref_counter is part of reproducability tests. The counter identifies a
+# The ref_counter is part of reproducibility tests. The counter identifies a
 # particular snapshot of our code, a "reference". Each PR is tested against this
 # reference to check if it reproduces the expected behavior. This allows us to
 # check that PRs that are expected to not modify the behavior do indeed preserve

src/ClimaAtmos.jl

@@ -17,7 +17,6 @@ import .Parameters as CAP
 
 include(joinpath("utils", "abbreviations.jl"))
 include(joinpath("utils", "gpu_compat.jl"))
-include(joinpath("utils", "common_spaces.jl"))
 include(joinpath("solver", "types.jl"))
 include(joinpath("solver", "cli_options.jl"))
 include(joinpath("utils", "utilities.jl"))

src/diagnostics/default_diagnostics.jl

@@ -169,31 +169,26 @@ function core_default_diagnostics(output_writer, duration, start_date; topograph
         min_func = (args...; kwargs...) -> hourly_min(FT, args...; kwargs...)
         max_func = (args...; kwargs...) -> hourly_max(FT, args...; kwargs...)
     end
-    core_diagnostics = if topography
-        [
-            # We need to compute the topography at the beginning of the
-            # simulation (and only at the beginning), so we set
-            # output/compute_schedule_func to false. It is still computed at the
-            # very beginning
-            ScheduledDiagnostic(;
-                variable = get_diagnostic_variable("orog"),
-                output_schedule_func = (integrator) -> false,
-                compute_schedule_func = (integrator) -> false,
-                output_writer,
-                output_short_name = "orog_inst",
-            ),
-            average_func(core_diagnostics...; output_writer, start_date)...,
-            min_func("ts"; output_writer, start_date),
-            max_func("ts"; output_writer, start_date)]
-    else
-        [
-            average_func(core_diagnostics...; output_writer, start_date)...,
-            min_func("ts"; output_writer, start_date),
-            max_func("ts"; output_writer, start_date),
-        ]
+    # Common diagnostics for all cases
+    common_diagnostics = [
+        average_func(core_diagnostics...; output_writer, start_date)...,
+        min_func("ts"; output_writer, start_date),
+        max_func("ts"; output_writer, start_date),
+    ]
+
+    # Prepend orography diagnostic if topography is enabled
+    if topography
+        orog_diagnostic = ScheduledDiagnostic(;
+            variable = get_diagnostic_variable("orog"),
+            output_schedule_func = (integrator) -> false,
+            compute_schedule_func = (integrator) -> false,
+            output_writer,
+            output_short_name = "orog_inst",
+        )
+        pushfirst!(common_diagnostics, orog_diagnostic)
     end
 
-    return core_diagnostics
+    return common_diagnostics
 end
 
 ##################

src/prognostic_equations/advection.jl

@@ -240,7 +240,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
 
     if point_type <: Geometry.Abstract3DPoint
         @. ᶜω³ = wcurlₕ(Y.c.uₕ)
-    elseif point_type <: Union{Geometry.Abstract2DPoint, Geometry.ZPoint}
+    else
         @. ᶜω³ = zero(ᶜω³)
     end