From 139879f866255dc8a482efe92294a7aa86a2c436 Mon Sep 17 00:00:00 2001 From: Dennis Yatunin Date: Thu, 25 Apr 2024 17:06:51 -0700 Subject: [PATCH] Add test for conservation of kinetic energy with topography --- .buildkite/pipeline.yml | 6 +++++ config/default_configs/default_config.yml | 3 +++ ...ountain_kinetic_energy_check_stretched.yml | 21 ++++++++++++++++ examples/hybrid/driver.jl | 12 ++++++++++ src/cache/cache.jl | 5 ++++ src/prognostic_equations/advection.jl | 24 +++++++++++++------ .../implicit/implicit_tendency.jl | 13 ++++++---- src/solver/type_getters.jl | 6 +++++ src/solver/types.jl | 1 + 9 files changed, 79 insertions(+), 12 deletions(-) create mode 100644 config/model_configs/plane_schar_mountain_kinetic_energy_check_stretched.yml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 3294f35610b..72d3d33c0b8 100644 --- a/.buildkite/pipeline.yml +++ b/.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: diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 341ae0739c3..5dc80d691f7 100644 --- a/config/default_configs/default_config.yml +++ b/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: ~ diff --git a/config/model_configs/plane_schar_mountain_kinetic_energy_check_stretched.yml b/config/model_configs/plane_schar_mountain_kinetic_energy_check_stretched.yml new file mode 100644 index 00000000000..4b64972eea3 --- /dev/null +++ b/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: "ScharProfile" +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 diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index a895d11eb28..5193a79b810 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -163,6 +163,18 @@ if config.parsed_args["check_conservation"] end end +# Kinetic energy conservation check +if config.parsed_args["check_kinetic_energy"] + ᶜρ_init = sol.u[1].c.ρ + ᶜρ_final = sol.u[end].c.ρ + + ᶜK_init = CA.compute_kinetic!(similar(ᶜρ_init), sol.u[1]) + ᶜK_final = CA.compute_kinetic!(similar(ᶜρ_init), sol.u[end]) + + FT = eltype(ᶜρ_init) + @test sum(ᶜρ_init .* ᶜK_init) ≈ sum(ᶜρ_final .* ᶜK_final) rtol = 2 * eps(FT) +end + # Precipitation characteristic checks if config.parsed_args["check_precipitation"] # run some simple tests based on the output diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 5e6d1da6273..53b005c7469 100644 --- a/src/cache/cache.jl +++ b/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 diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index ffa9d84ec1d..9186dc94993 100644 --- a/src/prognostic_equations/advection.jl +++ b/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 @@ -230,11 +233,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] += diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index aeed0f2cfdd..9942c677117 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -160,6 +160,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 @@ -213,11 +214,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 diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 5689fe3ba49..5c413c8d770 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -96,9 +96,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 diff --git a/src/solver/types.jl b/src/solver/types.jl index a1757835464..fa86ebe7eac 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -393,6 +393,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)