Merge #2091

2091: Use weak divergence in scalar advection r=simonbyrne a=simonbyrne



Co-authored-by: Simon Byrne <simonbyrne@gmail.com>

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) )
 ```

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

regression_tests/ref_counter.jl

@@ -1 +1 @@
-125
+126

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,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