From c3768522ef4d36c3b38e63b738de12d3f3a32514 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 | 10 ++++++++ src/cache/cache.jl | 5 ++++ src/initial_conditions/initial_conditions.jl | 10 +++++--- src/prognostic_equations/advection.jl | 24 +++++++++++++------ .../implicit/implicit_tendency.jl | 13 ++++++---- src/solver/type_getters.jl | 11 ++++++++- src/solver/types.jl | 1 + 10 files changed, 88 insertions(+), 16 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 d376ffe44a..80b5c6e175 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 d54016af88..c3b525a206 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 0000000000..5a61ed8eaf --- /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: "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 diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 688748d7d1..8c85e1ab6c 100644 --- a/examples/hybrid/driver.jl +++ b/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 diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 5e6d1da627..53b005c746 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/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index b3f6367af9..ca6e357b49 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/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 diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index ffa9d84ec1..9186dc9499 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 aeed0f2cfd..9942c67711 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 a8e8c19417..88011c4e6b 100644 --- a/src/solver/type_getters.jl +++ b/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"])", diff --git a/src/solver/types.jl b/src/solver/types.jl index 563ee7c696..5424c3a99b 100644 --- a/src/solver/types.jl +++ b/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)