Try using purely vertical spaces for column configurations

examples/Manifest.toml

@@ -353,9 +353,11 @@ weakdeps = ["CUDA", "MPI"]
 
 [[deps.ClimaCore]]
 deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"]
-git-tree-sha1 = "7209c2ed595f535db446709b9d19e22f8b000736"
+git-tree-sha1 = "137644370352036c05dd51d507df7921d680d34a"
+repo-rev = "tr/support-atmos-single-column"
+repo-url = "https://github.com/CliMA/ClimaCore.jl"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
-version = "0.14.22"
+version = "0.14.23"
 weakdeps = ["CUDA", "Krylov"]
 
     [deps.ClimaCore.extensions]
@@ -394,7 +396,9 @@ version = "0.7.6"
 
 [[deps.ClimaDiagnostics]]
 deps = ["Accessors", "ClimaComms", "ClimaCore", "Dates", "NCDatasets", "OrderedCollections", "SciMLBase"]
-git-tree-sha1 = "e9ac94af815dcae2a2ab24e54b53e76cca6258b7"
+git-tree-sha1 = "80456efa7642466e6bd967c6fd03f303df63d8b1"
+repo-rev = "tr/point-space-support"
+repo-url = "https://github.com/CliMA/ClimaDiagnostics.jl"
 uuid = "1ecacbb8-0713-4841-9a07-eb5aa8a2d53f"
 version = "0.2.11"
 

post_processing/ci_plots.jl

@@ -403,7 +403,7 @@ make_plots_generic(
     simulation_path,
     vars,
     time = LAST_SNAP,
-    x = 0.0, # Our columns are still 3D objects...
+    x = 0.0,
     y = 0.0,
     more_kwargs = YLINEARSCALE,
 )
@@ -419,7 +419,7 @@ make_plots_generic(
     simulation_path,
     vars,
     time = LAST_SNAP,
-    x = 0.0, # Our columns are still 3D objects...
+    x = 0.0,
     y = 0.0,
     more_kwargs = YLINEARSCALE,
 )
@@ -481,13 +481,16 @@ function make_plots(::ColumnPlots, output_paths::Vector{<:AbstractString})
     simdirs = SimDir.(output_paths)
     short_names = ["ta", "wa"]
     vars = map_comparison(get, simdirs, short_names)
-
+    slice_indices = ()
+    # columns were previously 3d, so they may need to be sliced for reproducibilty testing
+    vars = map(vars) do var
+        length(var.dims) > 2 ? slice(var, x = 0.0, y = 0.0) : var
+    end
     make_plots_generic(
         output_paths,
-        vars,
+        vars;
         time = LAST_SNAP,
-        x = 0.0, # Our columns are still 3D objects...
-        y = 0.0,
+        slice_indices...,
         MAX_NUM_COLS = length(simdirs),
         more_kwargs = YLINEARSCALE,
     )
@@ -521,10 +524,12 @@ function make_plots(
     simdir = simdirs[1]
 
     short_names = ["hus", "clw", "cli", "husra", "hussn", "ta"]
-    vars = [
-        slice(get(simdir; short_name), x = 0.0, y = 0.0) for
-        short_name in short_names
-    ]
+    is_purely_vertical = true
+    vars = [get(simdir; short_name) for short_name in short_names]
+    # make plotting backwards compatible with 3d columns
+    vars = map(vars) do var
+        length(var.dims) > 2 ? slice(var, x = 0.0, y = 0.0) : var
+    end
 
     # We first prepare the axes with all the nice labels with ClimaAnalysis, then we use
     # CairoMakie to add the additional lines.
@@ -565,11 +570,15 @@ function make_plots(
 
     # surface_precipitation
     surface_precip = read_var(simdir.variable_paths["pr"]["inst"]["10s"])
-    viz.line_plot1D!(
-        fig,
-        slice(surface_precip, x = 0.0, y = 0.0);
-        p_loc = [3, 1:3],
-    )
+    if length(surface_precip.dims) > 2
+        viz.line_plot1D!(
+            fig,
+            slice(surface_precip, x = 0.0, y = 0.0);
+            p_loc = [3, 1:3],
+        )
+    else
+        viz.line_plot1D!(fig, surface_precip; p_loc = [3, 1:3])
+    end
 
     file_path = joinpath(output_paths[1], "summary.pdf")
     CairoMakie.save(file_path, fig)
@@ -1325,13 +1334,17 @@ 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 [
-                slice(
-                    get(simdir; short_name, reduction, period),
-                    x = 0.0,
-                    y = 0.0,
-                ) for short_name in name_tuple
-            ]
+            grouped_vars = []
+            for short_name in name_tuple
+                var = get(simdir; short_name, reduction, period)
+                # If the var has less than two dimensions, it is a column (no horiontal space)
+                if length(var.dims) > 2
+                    push!(grouped_vars, slice(var, x = 0.0, y = 0.0))
+                else
+                    push!(grouped_vars, var)
+                end
+            end
+            return grouped_vars
         end
 
     var_groups_z = [

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -328,7 +328,9 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_do_integral!(
     @. ᶠΦ = CAP.grav(params) * ᶠz
     ᶜ∇Φ³ = p.scratch.ᶜtemp_CT3
     @. ᶜ∇Φ³ = CT3(ᶜgradᵥ(ᶠΦ))
-    @. ᶜ∇Φ³ += CT3(gradₕ(ᶜΦ))
+    if !iscolumn(axes(Y.c))
+        @. ᶜ∇Φ³ += CT3(gradₕ(ᶜΦ))
+    end
     ᶜ∇Φ₃ = p.scratch.ᶜtemp_C3
     @. ᶜ∇Φ₃ = ᶜgradᵥ(ᶠΦ)
 

src/callbacks/get_callbacks.jl

@@ -175,6 +175,23 @@ function get_diagnostics(
             )...,
             diagnostics...,
         ]
+        if !iscolumn(axes(Y.c))
+            # The topography is only defined when there is a horizontal grid.
+            # 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
+            diagnostics = [
+                CAD.ScheduledDiagnostic(;
+                    variable = CAD.get_diagnostic_variable("orog"),
+                    output_schedule_func = (integrator) -> false,
+                    compute_schedule_func = (integrator) -> false,
+                    output_writer = netcdf_writer,
+                    output_short_name = "orog_inst",
+                ),
+                diagnostics...,
+            ]
+
+        end
     end
     diagnostics = collect(diagnostics)
 

src/diagnostics/default_diagnostics.jl

@@ -167,16 +167,6 @@ function core_default_diagnostics(output_writer, duration, start_date)
     end
 
     return [
-        # 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),

src/parameterized_tendencies/les_sgs_models/smagorinsky_lilly.jl

@@ -52,11 +52,14 @@ function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
     # Gradients
     ## cell centers
     ∇ᶜu_uvw = @. ᶜtemp_UVWxUVW = Geometry.project(axis_uvw, ᶜgradᵥ(ᶠu_uvw))  # vertical component
-    @. ∇ᶜu_uvw += Geometry.project(axis_uvw, gradₕ(ᶜu_uvw))  # horizontal component
+    if !iscolumn(axes(Y.c))
+        @. ∇ᶜu_uvw += Geometry.project(axis_uvw, gradₕ(ᶜu_uvw))  # horizontal component
+    end
     ## cell faces
     ∇ᶠu_uvw = @. ᶠtemp_UVWxUVW = Geometry.project(axis_uvw, ᶠgradᵥ_uvw(ᶜu_uvw))  # vertical component
-    @. ∇ᶠu_uvw += Geometry.project(axis_uvw, gradₕ(ᶠu_uvw))  # horizontal component
-
+    if !iscolumn(axes(Y.c))
+        @. ∇ᶠu_uvw += Geometry.project(axis_uvw, gradₕ(ᶠu_uvw))  # horizontal component
+    end
     # Strain rate tensor
     ᶜS = @. ᶜtemp_strain = (∇ᶜu_uvw + adjoint(∇ᶜu_uvw)) / 2
     ᶠS = @. ᶠtemp_strain = (∇ᶠu_uvw + adjoint(∇ᶠu_uvw)) / 2
@@ -72,10 +75,14 @@ function set_smagorinsky_lilly_precomputed_quantities!(Y, p)
     ᶜfb = @. ᶜtemp_scalar = ifelse(ᶜRi ≤ 0, 1, max(0, 1 - ᶜRi / Pr_t)^(1 / 4))
 
     # filter scale
-    h_space = Spaces.horizontal_space(axes(Y.c))
-    Δ_xy = Spaces.node_horizontal_length_scale(h_space)^2 # Δ_x * Δ_y
-    ᶜΔ_z = Fields.Δz_field(Y.c)
-    ᶜΔ = @. ᶜtemp_scalar = ∛(Δ_xy * ᶜΔ_z) * ᶜfb
+    if !iscolumn(axes(Y.c))
+        h_space = Spaces.horizontal_space(axes(Y.c))
+        Δ_xy = Spaces.node_horizontal_length_scale(h_space)^2 # Δ_x * Δ_y
+        ᶜΔ_z = Fields.Δz_field(Y.c)
+        ᶜΔ = @. ᶜtemp_scalar = ∛(Δ_xy * ᶜΔ_z) * ᶜfb
+    else
+        ᶜΔ = @. ᶜtemp_scalar = FT(0)
+    end
 
     # Smagorinsky-Lilly eddy viscosity
     ᶜνₜ = @. ᶜtemp_scalar = c_smag^2 * ᶜΔ^2 * ᶜS_norm

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -278,7 +278,9 @@ function compute_precipitation_heating!(
 
     # compute full temperature gradient
     @. ᶜ∇T = CT123(ᶜgradᵥ(ᶠinterp(Tₐ(thp, ᶜts))))
-    @. ᶜ∇T += CT123(gradₕ(Tₐ(thp, ᶜts)))
+    if !iscolumn(axes(ᶜSeₜᵖ))
+        @. ᶜ∇T += CT123(gradₕ(Tₐ(thp, ᶜts)))
+    end
     # dot product with effective velocity of precipitation
     # (times q and specific heat)
     @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwᵣ)))) * cᵥₗ(thp) * ᶜqᵣ

src/parameterized_tendencies/radiation/radiation.jl

@@ -259,7 +259,8 @@ function radiation_model_cache(
             rrtmgp_params,
             data_loader,
             context;
-            ncol = length(Spaces.all_nodes(axes(Spaces.level(Y.c, 1)))),
+            ncol = iscolumn(axes(Y.c)) ? 1 :
+                   length(Spaces.all_nodes(axes(Spaces.level(Y.c, 1)))),
             domain_nlay = Spaces.nlevels(axes(Y.c)),
             radiation_mode,
             interpolation,

src/prognostic_equations/advection.jl

@@ -91,17 +91,20 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
 
     if point_type <: Geometry.Abstract3DPoint
         @. ᶜω³ = curlₕ(Y.c.uₕ)
-    elseif point_type <: Geometry.Abstract2DPoint
+    elseif point_type <: Geometry.Abstract2DPoint ||
+           point_type <: Geometry.Abstract1DPoint
         @. ᶜω³ = zero(ᶜω³)
     end
 
     @. ᶠω¹² = ᶠcurlᵥ(Y.c.uₕ)
     for j in 1:n
         @. ᶠω¹²ʲs.:($$j) = ᶠω¹²
     end
-    @. ᶠω¹² += CT12(curlₕ(Y.f.u₃))
-    for j in 1:n
-        @. ᶠω¹²ʲs.:($$j) += CT12(curlₕ(Y.f.sgsʲs.:($$j).u₃))
+    if !iscolumn(axes(Y.c))
+        @. ᶠω¹² += CT12(curlₕ(Y.f.u₃))
+        for j in 1:n
+            @. ᶠω¹²ʲs.:($$j) += CT12(curlₕ(Y.f.sgsʲs.:($$j).u₃))
+        end
     end
     # Without the CT12(), the right-hand side would be a CT1 or CT2 in 2D space.
 

src/prognostic_equations/forcing/external_forcing.jl

@@ -103,9 +103,7 @@ function external_forcing_cache(Y, external_forcing::GCMForcing, params)
             parent(cc_field) .= ds.group[cfsite_number]["coszen"][1]
         end
 
-        zc_forcing = gcm_height(ds.group[cfsite_number])
-        Fields.bycolumn(axes(Y.c)) do colidx
-
+        function set_column_data(colidx)
             zc_gcm = Fields.coordinate_field(Y.c).z[colidx]
 
             setvar!(ᶜdTdt_hadv, "tntha", colidx, zc_gcm, zc_forcing)
@@ -139,6 +137,15 @@ function external_forcing_cache(Y, external_forcing::GCMForcing, params)
             @. ᶜinv_τ_wind[colidx] = 1 / (6 * 3600)
             @. ᶜinv_τ_scalar[colidx] = compute_gcm_driven_scalar_inv_τ(zc_gcm)
         end
+
+        zc_forcing = gcm_height(ds.group[cfsite_number])
+        if iscolumn(axes(Y.c))
+            set_column_data(:)
+        else
+            Fields.bycolumn(axes(Y.c)) do colidx
+                set_column_data(colidx)
+            end
+        end
     end
 
     return (;

src/prognostic_equations/hyperdiffusion.jl

@@ -13,8 +13,6 @@ hyperdiffusion_cache(Y, atmos) =
 hyperdiffusion_cache(Y, hyperdiff::Nothing, _) = (;)
 
 function hyperdiffusion_cache(Y, hyperdiff::ClimaHyperdiffusion, turbconv_model)
-    quadrature_style =
-        Spaces.quadrature_style(Spaces.horizontal_space(axes(Y.c)))
     FT = eltype(Y)
     n = n_mass_flux_subdomains(turbconv_model)
 

src/prognostic_equations/remaining_tendency.jl

@@ -13,17 +13,22 @@ end
 NVTX.@annotate function remaining_tendency!(Yₜ, Yₜ_lim, Y, p, t)
     Yₜ_lim .= zero(eltype(Yₜ_lim))
     Yₜ .= zero(eltype(Yₜ))
-    horizontal_tracer_advection_tendency!(Yₜ_lim, Y, p, t)
-    fill_with_nans!(p)
-    horizontal_advection_tendency!(Yₜ, Y, p, t)
-    hyperdiffusion_tendency!(Yₜ, Yₜ_lim, Y, p, t)
+    # only calulate horizontal tendencies if there is a horizontal space
+    if !iscolumn(axes(Y.c))
+        horizontal_tracer_advection_tendency!(Yₜ_lim, Y, p, t)
+        fill_with_nans!(p)
+        horizontal_advection_tendency!(Yₜ, Y, p, t)
+        hyperdiffusion_tendency!(Yₜ, Yₜ_lim, Y, p, t)
+    end
     explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     additional_tendency!(Yₜ, Y, p, t)
     return Yₜ
 end
 
 NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
-    viscous_sponge_tendency!(Yₜ, Y, p, t, p.atmos.viscous_sponge)
+    if !iscolumn(axes(Y.c))
+        viscous_sponge_tendency!(Yₜ, Y, p, t, p.atmos.viscous_sponge)
+    end
 
     # Vertical tendencies
     rayleigh_sponge_tendency!(Yₜ, Y, p, t, p.atmos.rayleigh_sponge)
@@ -87,7 +92,9 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     pressure_work_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 
     sl = p.atmos.smagorinsky_lilly
-    horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl)
+    if !iscolumn(axes(Y.c))
+        horizontal_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl)
+    end
     vertical_smagorinsky_lilly_tendency!(Yₜ, Y, p, t, sl)
 
     # NOTE: This will zero out all momentum tendencies in the edmfx advection test

src/solver/model_getters.jl

@@ -43,6 +43,10 @@ end
 function get_hyperdiffusion_model(parsed_args, ::Type{FT}) where {FT}
     hyperdiff_name = parsed_args["hyperdiff"]
     if hyperdiff_name in ("ClimaHyperdiffusion", "true", true)
+        if parsed_args["config"] == "column"
+            @warn "Hyperdiffusion should not be used with column configuration\
+            The `hyperdiffusion_tendency!` function will not be called."
+        end
         ν₄_vorticity_coeff =
             FT(parsed_args["vorticity_hyperdiffusion_coefficient"])
         ν₄_scalar_coeff = FT(parsed_args["scalar_hyperdiffusion_coefficient"])
@@ -60,6 +64,10 @@ function get_hyperdiffusion_model(parsed_args, ::Type{FT}) where {FT}
         @info "Using CAM_SE hyperdiffusion. vorticity_hyperdiffusion_coefficient, \
                scalar_hyperdiffusion_coefficient and divergence_damping_factor in the config \
                will be ignored."
+        if parsed_args["config"] == "column"
+            @warn "Hyperdiffusion should not be used with column configuration. \
+            The `hyperdiffusion_tendency!` function will not be called."
+        end
         ν₄_vorticity_coeff = FT(0.150 * 1.238)
         ν₄_scalar_coeff = FT(0.751 * 1.238)
         divergence_damping_factor = FT(5)