Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use weak divergence in scalar advection #2091

Merged
merged 3 commits into from
Sep 19, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 3 additions & 4 deletions docs/src/equations.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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) ).
```
Expand Down Expand Up @@ -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) )
```
Expand Down
2 changes: 1 addition & 1 deletion examples/hybrid/driver.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion regression_tests/ref_counter.jl
Original file line number Diff line number Diff line change
@@ -1 +1 @@
125
126
18 changes: 9 additions & 9 deletions src/prognostic_equations/advection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 + ᶜΦ))
Expand All @@ -59,15 +59,15 @@ 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
for j in 1:n
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
Expand Down
Loading