From 01f7455ffd78a6cc4b063ad95a7c9b9f10e452b3 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Fri, 15 Sep 2023 10:28:28 -0700 Subject: [PATCH 1/3] Use weak divergence in scalar advection --- docs/src/equations.md | 7 +++---- src/prognostic_equations/advection.jl | 18 +++++++++--------- 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/docs/src/equations.md b/docs/src/equations.md index 181eb2b851..62f073a34f 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -33,7 +33,6 @@ We make use of the following operators ### Differential operators -- ``\mathcal{D}_h`` is the [discrete horizontal spectral divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.Divergence). - ``\hat{\mathcal{D}}_h`` is the [discrete horizontal spectral weak divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.WeakDivergence). - ``\mathcal{D}^c_v`` is the [face-to-center vertical divergence](https://clima.github.io/ClimaCore.jl/stable/operators/#ClimaCore.Operators.DivergenceF2C). !!! todo @@ -107,7 +106,7 @@ Follows the continuity equation This is discretized using the following ```math \frac{\partial}{\partial t} \rho -= - \mathcal{D}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right] += - \hat{\mathcal{D}}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right] ``` with the @@ -210,7 +209,7 @@ which is stabilized with the addition of a 4th-order hyperdiffusion term on tota This is discretized using ```math \frac{\partial}{\partial t} \rho e \approx -- \mathcal{D}_h[ (\rho e + p) \bar{\boldsymbol{u}} ] +- \hat{\mathcal{D}}_h[ (\rho e + p) \bar{\boldsymbol{u}} ] - \mathcal{D}^c_v \left[ WI^f(J,\rho) \, \tilde{\boldsymbol{u}} \, I^f \left(\frac{\rho e + p}{\rho} \right) + \boldsymbol{F}_R \right] - \nu_h \hat{\mathcal{D}}_h( \rho \mathcal{G}_h(\psi) ). ``` @@ -247,7 +246,7 @@ This is stabilized with the addition of a 4th-order hyperdiffusion term This is discretized using ```math \frac{\partial}{\partial t} \rho \chi \approx -- \mathcal{D}_h[ \rho \chi \bar{\boldsymbol{u}}] +- \hat{\mathcal{D}}_h[ \rho \chi \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[ WI^f(J,\rho) \, U^f\left( \tilde{\boldsymbol{u}}, \frac{\rho \chi}{\rho} \right) \right] - \nu_\chi \hat{\mathcal{D}}_h ( \rho \, \mathcal{G}_h (\psi) ) ``` diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index aff688bd18..9e73101eaa 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -16,34 +16,34 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) (; ᶜuʲs) = p end - @. Yₜ.c.ρ -= divₕ(Y.c.ρ * ᶜu) + @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu) if p.atmos.turbconv_model isa EDMFX for j in 1:n - @. Yₜ.c.sgsʲs.:($$j).ρa -= divₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j)) + @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j)) end end if :ρθ in propertynames(Y.c) - @. Yₜ.c.ρθ -= divₕ(Y.c.ρθ * ᶜu) + @. Yₜ.c.ρθ -= wdivₕ(Y.c.ρθ * ᶜu) elseif :ρe_tot in propertynames(Y.c) (; ᶜh_tot) = p - @. Yₜ.c.ρe_tot -= divₕ(Y.c.ρ * ᶜh_tot * ᶜu) + @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu) end if p.atmos.turbconv_model isa EDMFX for j in 1:n if :ρθ in propertynames(Y.c) @. Yₜ.c.sgsʲs.:($$j).ρaθ -= - divₕ(Y.c.sgsʲs.:($$j).ρaθ * ᶜuʲs.:($$j)) + wdivₕ(Y.c.sgsʲs.:($$j).ρaθ * ᶜuʲs.:($$j)) elseif :ρe_tot in propertynames(Y.c) (; ᶜh_totʲs) = p @. Yₜ.c.sgsʲs.:($$j).ρae_tot -= - divₕ(Y.c.sgsʲs.:($$j).ρa * ᶜh_totʲs.:($$j) * ᶜuʲs.:($$j)) + wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜh_totʲs.:($$j) * ᶜuʲs.:($$j)) end end end if use_prognostic_tke(p.atmos.turbconv_model) - @. Yₜ.c.sgs⁰.ρatke -= divₕ(Y.c.sgs⁰.ρatke * ᶜu⁰) + @. Yₜ.c.sgs⁰.ρatke -= wdivₕ(Y.c.sgs⁰.ρatke * ᶜu⁰) end @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp - ᶜp_ref) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ)) @@ -59,7 +59,7 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t) end for ρχ_name in filter(is_tracer_var, propertynames(Y.c)) - @. Yₜ.c.:($$ρχ_name) -= divₕ(Y.c.:($$ρχ_name) * ᶜu) + @. Yₜ.c.:($$ρχ_name) -= wdivₕ(Y.c.:($$ρχ_name) * ᶜu) end if p.atmos.turbconv_model isa EDMFX @@ -67,7 +67,7 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t) for ρaχ_name in filter(is_tracer_var, propertynames(Y.c.sgsʲs.:($j))) @. Yₜ.c.sgsʲs.:($$j).:($$ρaχ_name) -= - divₕ(Y.c.sgsʲs.:($$j).:($$ρaχ_name) * ᶜuʲs.:($$j)) + wdivₕ(Y.c.sgsʲs.:($$j).:($$ρaχ_name) * ᶜuʲs.:($$j)) end end end From 3e7a848b0ce98da43d26a190c7571cbff37749e7 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Mon, 18 Sep 2023 13:01:28 -0700 Subject: [PATCH 2/3] bump ref --- regression_tests/ref_counter.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/regression_tests/ref_counter.jl b/regression_tests/ref_counter.jl index d136d6a714..0a3e7b0422 100644 --- a/regression_tests/ref_counter.jl +++ b/regression_tests/ref_counter.jl @@ -1 +1 @@ -125 +126 From ce41d13a23dd483f14b1279445f0366161b55689 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Mon, 18 Sep 2023 17:00:25 -0700 Subject: [PATCH 3/3] increase conservation threshold --- examples/hybrid/driver.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index a119a1f578..7e9db26786 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -173,5 +173,5 @@ if config.parsed_args["check_conservation"] @test sum(sol.u[1].c.ρ) ≈ sum(sol.u[end].c.ρ) rtol = 50 * eps(FT) @test sum(sol.u[1].c.ρe_tot) + (p.net_energy_flux_sfc[][] - p.net_energy_flux_toa[][]) ≈ - sum(sol.u[end].c.ρe_tot) rtol = 30 * eps(FT) + sum(sol.u[end].c.ρe_tot) rtol = 100 * eps(FT) end