Add test for conservation of kinetic energy with topography

.buildkite/pipeline.yml

@@ -181,6 +181,12 @@ steps:
           --config_file $CONFIG_PATH/plane_density_current_test.yml
         artifact_paths: "plane_density_current_test/output_active/*"
 
+      - label: ":computer: Kinetic energy conservation check for Schar mountain experiment (stretched grid)"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $CONFIG_PATH/plane_schar_mountain_kinetic_energy_check_stretched.yml
+        artifact_paths: "plane_schar_mountain_kinetic_energy_check_stretched/output_active/*"
+
 
   - group: "Sphere Examples (Dycore)"
     steps:

config/default_configs/default_config.yml

@@ -226,6 +226,9 @@ check_conservation:
 check_precipitation:
   help: "Sanity checks for 1-moment precipitation [`false` (default), `true`]"
   value: false
+check_kinetic_energy:
+  help: "Whether to disable the effects of pressure, gravity, and the Coriolis force on velocity, and also check that kinetic energy is conserved at the end of the simulation"
+  value: false
 ls_adv:
   help: "Large-scale advection [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `ARM_SGP`, `GATE_III`]"
   value: ~

config/model_configs/plane_schar_mountain_kinetic_energy_check_stretched.yml

@@ -0,0 +1,21 @@
+check_kinetic_energy: true
+dt_save_state_to_disk: "3600secs"
+initial_condition: "ModifiedScharProfile"
+x_max: 120000.0
+z_elem: 45
+dt: "0.5secs"
+t_end: "120secs"
+dz_top: 1000.0
+x_elem: 100
+dz_bottom: 200.0
+use_reference_state: false
+ode_algo: "SSP33ShuOsher"
+config: "plane"
+hyperdiff: "false"
+z_max: 25000.0
+topography: "Schar"
+job_id: "plane_schar_mountain_kinetic_energy_check_stretched"
+toml: [toml/plane_schar_mountain_test_stretched.toml]
+diagnostics:
+  - short_name: wa
+    period: 4hours

examples/hybrid/driver.jl

@@ -162,6 +162,16 @@ if config.parsed_args["check_conservation"]
     end
 end
 
+# Kinetic energy conservation check
+if config.parsed_args["check_kinetic_energy"]
+    Y_0 = sol.u[1]
+    Y_1 = sol.u[end]
+    ρK_0 = sum(Y_0.c.ρ .* CA.compute_kinetic!(similar(Y_0.c.ρ), Y_0))
+    ρK_1 = sum(Y_1.c.ρ .* CA.compute_kinetic!(similar(Y_1.c.ρ), Y_1))
+    @info "    Relative kinetic energy change: $((ρK_1 - ρK_0) / ρK_0)"
+    @test abs(ρK_1 - ρK_0) <= 2 * eps(ρK_0)
+end
+
 # Precipitation characteristic checks
 if config.parsed_args["check_precipitation"]
     # run some simple tests based on the output

src/cache/cache.jl

@@ -159,6 +159,11 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
         ᶠf¹² = nothing
     end
 
+    if atmos.check_kinetic_energy
+        @. ᶜf³ *= 0
+        isnothing(ᶠf¹²) || @. ᶠf¹² *= 0
+    end
+
     quadrature_style =
         Spaces.quadrature_style(Spaces.horizontal_space(axes(Y.c)))
     do_dss = quadrature_style isa Quadratures.GLL

src/initial_conditions/initial_conditions.jl

@@ -197,9 +197,12 @@ end
 
 An `InitialCondition` with a prescribed Brunt-Vaisala Frequency
 """
-Base.@kwdef struct ScharProfile <: InitialCondition end
+Base.@kwdef struct ScharProfile <: InitialCondition
+    constant_velocity::Bool = true
+end
 
 function (initial_condition::ScharProfile)(params)
+    (; constant_velocity) = initial_condition
     function local_state(local_geometry)
         FT = eltype(params)
 
@@ -219,13 +222,14 @@ function (initial_condition::ScharProfile)(params)
         T = π_exner * θ # temperature
         ρ = p₀ / (R_d * T) * (π_exner)^(cp_d / R_d)
         p = ρ * R_d * T
-        velocity = Geometry.UVVector(FT(10), FT(0))
+        u = constant_velocity ? FT(10) : FT(0.5) * sin(z)
+        v = constant_velocity ? FT(0) : FT(500) * sin(x / 1000)
 
         return LocalState(;
             params,
             geometry = local_geometry,
             thermo_state = TD.PhaseDry_pT(thermo_params, p, T),
-            velocity = velocity,
+            velocity = Geometry.UVVector(u, v),
         )
     end
     return local_state

src/prognostic_equations/advection.jl

@@ -39,8 +39,11 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
         @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰)
     end
 
-    @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp - ᶜp_ref) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ))
-    # Without the C12(), the right-hand side would be a C1 or C2 in 2D space.
+    if p.atmos.check_kinetic_energy
+        @. Yₜ.c.uₕ -= C12(gradₕ(ᶜK))
+    else
+        @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp - ᶜp_ref) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ))
+    end # Without the C12(), the right-hand side would be a C1 or C2 in 2D space.
     return nothing
 end
 
@@ -225,11 +228,18 @@ function edmfx_sgs_vertical_advection_tendency!(
             ᶜleft_bias(ᶠKᵥʲs.:($$j)[colidx]),
             ᶜright_bias(ᶠKᵥʲs.:($$j)[colidx]),
         )
-        # For the updraft u_3 equation, we assume the grid-mean to be hydrostatic
-        # and calcuate the buoyancy term relative to the grid-mean density.
-        @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
-            (ᶠinterp(ᶜρʲs.:($$j)[colidx] - Y.c.ρ[colidx]) * ᶠgradᵥ_ᶜΦ[colidx]) /
-            ᶠinterp(ᶜρʲs.:($$j)[colidx]) + ᶠgradᵥ(ᶜKᵥʲ[colidx])
+
+        if p.atmos.check_kinetic_energy
+            @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -= ᶠgradᵥ(ᶜKᵥʲ[colidx])
+        else
+            # Assume the grid-mean to be hydrostatic and calcuate the buoyancy
+            # term relative to the grid-mean density.
+            @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] -=
+                (
+                    ᶠinterp(ᶜρʲs.:($$j)[colidx] - Y.c.ρ[colidx]) *
+                    ᶠgradᵥ_ᶜΦ[colidx]
+                ) / ᶠinterp(ᶜρʲs.:($$j)[colidx]) + ᶠgradᵥ(ᶜKᵥʲ[colidx])
+        end
 
         # buoyancy term in mse equation
         @. Yₜ.c.sgsʲs.:($$j).mse[colidx] +=

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -135,6 +135,7 @@ vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:third_order}) =
 
 function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     (; moisture_model, turbconv_model, rayleigh_sponge, precip_model) = p.atmos
+    (; check_kinetic_energy) = p.atmos
     (; dt) = p
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜJ = Fields.local_geometry_field(Y.c).J
@@ -185,11 +186,13 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
         )
     end
 
-    @. Yₜ.f.u₃[colidx] +=
-        -(
-            ᶠgradᵥ(ᶜp[colidx] - ᶜp_ref[colidx]) +
-            ᶠinterp(Y.c.ρ[colidx] - ᶜρ_ref[colidx]) * ᶠgradᵥ_ᶜΦ[colidx]
-        ) / ᶠinterp(Y.c.ρ[colidx])
+    if !check_kinetic_energy
+        @. Yₜ.f.u₃[colidx] +=
+            -(
+                ᶠgradᵥ(ᶜp[colidx] - ᶜp_ref[colidx]) +
+                ᶠinterp(Y.c.ρ[colidx] - ᶜρ_ref[colidx]) * ᶠgradᵥ_ᶜΦ[colidx]
+            ) / ᶠinterp(Y.c.ρ[colidx])
+    end
 
     if rayleigh_sponge isa RayleighSponge
         (; ᶠβ_rayleigh_w) = p.rayleigh_sponge

src/solver/type_getters.jl

@@ -95,9 +95,15 @@ function get_atmos(config::AtmosConfig, params)
         surface_model = get_surface_model(parsed_args),
         surface_albedo = get_surface_albedo_model(parsed_args, params, FT),
         numerics = get_numerics(parsed_args),
+        check_kinetic_energy = parsed_args["check_kinetic_energy"],
     )
     @assert !@any_reltype(atmos, (UnionAll, DataType))
 
+    if atmos.check_kinetic_energy
+        @assert isnothing(atmos.hyperdiff)
+        @assert isnothing(atmos.rayleigh_sponge)
+    end
+
     @info "AtmosModel: \n$(summary(atmos))"
     return atmos
 end
@@ -141,7 +147,8 @@ function get_spaces(parsed_args, params, comms_ctx)
     bubble = parsed_args["bubble"]
     deep = parsed_args["deep_atmosphere"]
 
-    @assert topography in ("NoWarp", "DCMIP200", "Earth", "Agnesi", "Schar")
+    @assert topography in
+            ("NoWarp", "DCMIP200", "Earth", "Agnesi", "Tall Schar", "Schar")
     if topography == "DCMIP200"
         warp_function = topography_dcmip200
     elseif topography == "Agnesi"
@@ -354,6 +361,8 @@ function get_initial_condition(parsed_args)
         "PrecipitatingColumn",
     ]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
+    elseif parsed_args["initial_condition"] == "ModifiedScharProfile"
+        return ICs.ScharProfile(; constant_velocity = false)
     else
         error(
             "Unknown `initial_condition`: $(parsed_args["initial_condition"])",

src/solver/types.jl

@@ -399,6 +399,7 @@ Base.@kwdef struct AtmosModel{
     surface_model::SM = nothing
     surface_albedo::SA = nothing
     numerics::NUM = nothing
+    check_kinetic_energy::Bool
 end
 
 Base.broadcastable(x::AtmosModel) = tuple(x)