test split divergence

Author: nefrathenrici

.buildkite/pipeline.yml

@@ -28,6 +28,7 @@ steps:
       - julia --check-bounds=yes -e 'using TOML; @assert TOML.parsefile(".buildkite/Manifest-v1.11.toml")["deps"]["ClimaAtmos"][1]["path"] == ".."'
 
       - echo "--- Instantiate .buildkite"
+      - julia --project=.buildkite -e 'using Pkg; Pkg.add(PackageSpec(name="ClimaCore", rev="ne-ts/splitform"))'
       - "julia --project=.buildkite -e 'using Pkg; Pkg.instantiate(;verbose=true); Pkg.precompile(;strict=true); using CUDA; CUDA.precompile_runtime(); Pkg.status()'"
 
     agents:

src/prognostic_equations/advection.jl

@@ -38,14 +38,16 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t)
     (; ᶜΦ) = p.core
     (; ᶜu, ᶜK, ᶜp, ᶜts) = p.precomputed
     (; params) = p
+    FT = eltype(params)
     thermo_params = CAP.thermodynamics_params(params)
     cp_d = thermo_params.cp_d
 
     if p.atmos.turbconv_model isa PrognosticEDMFX
         (; ᶜuʲs) = p.precomputed
     end
 
-    @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu)
+    @. Yₜ.c.ρ -= split_divₕ(Y.c.ρ * ᶜu, FT(1))  # This one
+    
     if p.atmos.turbconv_model isa PrognosticEDMFX
         for j in 1:n
             @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j))
@@ -54,7 +56,7 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t)
 
     ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, Y.c.ρ))
     ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot))
-    @. Yₜ.c.ρe_tot -= wdivₕ(Y.c.ρ * ᶜh_tot * ᶜu)
+    @. Yₜ.c.ρe_tot -= split_divₕ(Y.c.ρ * ᶜu, ᶜh_tot) # this one
 
     if p.atmos.turbconv_model isa PrognosticEDMFX
         for j in 1:n
@@ -73,8 +75,16 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t)
     ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts))
     ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts))
     ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts))
-    @. Yₜ.c.uₕ -= C12(gradₕ(ᶜK + ᶜΦ - ᶜΦ_r) + cp_d * (ᶜθ_v - ᶜθ_vr) * gradₕ(ᶜΠ))
-    # Without the C12(), the right-hand side would be a C1 or C2 in 2D space.
+    ᶜθ_v_diff = @. lazy(ᶜθ_v - ᶜθ_vr)
+    # PG = 0.5 * cp_d * [θv ∇Π + ∇(θv Π) - Π∇θv]
+    @. Yₜ.c.uₕ -= C12(
+        gradₕ(ᶜK + ᶜΦ - ᶜΦ_r) +
+        FT(0.5) * cp_d * (
+            ᶜθ_v_diff * gradₕ(ᶜΠ) +  # θv ∇Π
+            gradₕ(ᶜθ_v_diff * ᶜΠ) -  # ∇(θv Π)
+            ᶜΠ * gradₕ(ᶜθ_v_diff)    # Π∇θv
+        ),
+    )    # Without the C12(), the right-hand side would be a C1 or C2 in 2D space.
     return nothing
 end
 
@@ -109,8 +119,9 @@ NVTX.@annotate function horizontal_tracer_advection_tendency!(Yₜ, Y, p, t)
         (; ᶜuʲs) = p.precomputed
     end
 
-    for ρχ_name in filter(is_tracer_var, propertynames(Y.c))
-        @. Yₜ.c.:($$ρχ_name) -= wdivₕ(Y.c.:($$ρχ_name) * ᶜu)
+    for ρχ_name in filter(is_tracer_var, propertynames(Y.c)) # this one
+        ᶜχ = @. lazy(specific(Y.c.:($$ρχ_name), Y.c.ρ))
+        @. Yₜ.c.:($$ρχ_name) -= split_divₕ(Y.c.ρ * ᶜu, ᶜχ)
     end
 
     if p.atmos.turbconv_model isa PrognosticEDMFX

src/utils/abbreviations.jl

@@ -16,6 +16,7 @@ const UVW = Geometry.UVWVector
 
 const divₕ = Operators.Divergence()
 const wdivₕ = Operators.WeakDivergence()
+const split_divₕ = Operators.SplitDivergence()
 const gradₕ = Operators.Gradient()
 const wgradₕ = Operators.WeakGradient()
 const curlₕ = Operators.Curl()