Merge pull request #3975 from CliMA/as/update-pressure-ref

Update PGF + modify hyperdiffusion reference state.

Author: ray-chew

NEWS.md

@@ -4,6 +4,10 @@ ClimaAtmos.jl Release Notes
 main
 -------
 
+v0.31.5
+-------
+PR [#3975](https://github.com/CliMA/ClimaAtmos.jl/pull/3975) updates the pressure gradient formulation to subtract a reference state and use the Exner pressure.
+
 v0.31.4
 -------
 

config/model_configs/box_hydrostatic_balance.yml

@@ -1,6 +1,6 @@
 t_end: "1days"
 config: "box"
-hyperdiff: false
+hyperdiff: true
 dt: "20secs"
 perturb_initstate: false
 disable_surface_flux_tendency: true

docs/src/equations.md

@@ -61,36 +61,32 @@ We make use of the following operators
      ```math
      \boldsymbol{\Omega} = \Omega \sin(\phi) \boldsymbol{e}^v
      ```
-     where ``\phi`` is latitude, and ``\Omega`` is the planetary rotation rate in rads/sec (for Earth, ``7.29212 \times 10^{-5} s^{-1}``) and ``\boldsymbol{e}^v`` is the unit radial basis vector. This implies that the horizontal contravariant component ``\boldsymbol{\Omega}^h`` is zero.
+      where ``\phi`` is latitude, and ``\Omega`` is the planetary rotation rate in rads/sec (for Earth, ``7.29212 \times 10^{-5} s^{-1}``) and ``\boldsymbol{e}^v`` is the unit radial basis vector. This implies that the horizontal contravariant component ``\boldsymbol{\Omega}^h`` is zero.
+
   *  a _deep atmosphere_, with
      ```math
      \boldsymbol{\Omega} = (0, 0, \Omega)
      ```
      i.e. aligned with Earth's rotational axis.
 * ``\tilde{\boldsymbol{u}}`` is the mass-weighted reconstruction of velocity at the interfaces:
   by interpolation of contravariant components
-    ```math
-  \tilde{\boldsymbol{u}} = WI^f(\rho J, \boldsymbol{u}_h) + \boldsymbol{u}_v
+  ```math
+
   ```
 * ``\bar{\boldsymbol{u}}`` is the reconstruction of velocity at cell-centers, carried out by linear interpolation of the covariant vertical component:
   ```math
   \bar{\boldsymbol{u}} = \boldsymbol{u}_h + I_{c}(\boldsymbol{u}_v)
   ```
 
 * ``\Phi = g z`` is the geopotential, where ``g`` is the gravitational acceleration rate and ``z`` is altitude above the mean sea level.
-* ``\boldsymbol{b}`` is the reduced gravitational acceleration
-  ```math
-  \boldsymbol{b} = - \frac{\rho - \rho_{\text{ref}}}{\rho} \nabla \Phi
-  ```
-* ``\rho_{\text{ref}}`` is the reference state density
 * ``K = \tfrac{1}{2} \|\boldsymbol{u}\|^2 `` is the specific kinetic energy (J/kg), reconstructed at cell centers by
   ```math
   K = \tfrac{1}{2} (\boldsymbol{u}_{h} \cdot \boldsymbol{u}_{h} + 2 \boldsymbol{u}_{h} \cdot I_{c} (\boldsymbol{u}_{v}) + I_{c}(\boldsymbol{u}_{v} \cdot \boldsymbol{u}_{v})),
   ```
   where ``\boldsymbol{u}_{h}`` is defined on cell-centers, ``\boldsymbol{u}_{v}`` is defined on cell-faces, and ``I_{c} (\boldsymbol{u}_{v})`` is interpolated using covariant components.
 
 * ``p`` is air pressure, derived from the thermodynamic state, reconstructed at cell centers.
-* ``p_{\text{ref}}`` is the reference state pressure. It is related to the reference state density by analytical hydrostatic balance: ``\nabla p_{\text{ref}} = - \rho_{\text{ref}} \nabla \Phi``.
+* ``\Pi = (\frac{p}{p_0})^{\frac{R_d}{c_{pd}}}`` is the Exner function evaluated with dry-air constants.
 * ``\boldsymbol{F}_R`` are the radiative fluxes: these are assumed to align vertically (i.e. the horizontal contravariant components are zero), and are constructed at cell faces from [RRTMGP.jl](https://github.com/CliMA/RRTMGP.jl).
 
 * ``\nu_u``, ``\nu_h``, and ``\nu_\chi`` are hyperdiffusion coefficients, and ``c`` is the divergence damping factor.
@@ -126,20 +122,31 @@ term treated implicitly (check this)
 
 Uses the advective form equation
 ```math
-\frac{\partial}{\partial t} \boldsymbol{u}  = - (2 \boldsymbol{\Omega} + \nabla \times \boldsymbol{u}) \times \boldsymbol{u} - \frac{1}{\rho} \nabla (p - p_{\text{ref}})  + \boldsymbol{b} - \nabla K
+\frac{\partial}{\partial t} \boldsymbol{u}  = - (2 \boldsymbol{\Omega} + \nabla \times \boldsymbol{u}) \times \boldsymbol{u} - c_{pd} (\theta_v - \theta_{v, r}) \nabla_h \Pi  - \nabla_h [(\Phi - \Phi_r) + K]. 
+```
+Here, we use the Exner function to compute pressure gradients and are subtracting a hydrostatic reference state
+```math
+- \frac{1}{\rho} \nabla p = - c_{pd} \theta_v \Pi
+```
+where ``\theta_v`` is the virtual potential temperature. ``\theta_{v,r} = T_r / \Pi`` is a reference virtual potential temperature (with reference temperature ``T_r``), and 
+```math
+\Phi_r = -c_{pd} \left[ T_\text{min} \log(\Pi) + \frac{(T_\text{sfc} - T_\text{min})}{n_s} (\Pi^{n_s} - 1) \right],
 ```
+is a reference geopotential, which satisfies the hydrostatic balance equation $c_{pd} \theta_{v,r} \nabla \Pi + \nabla \Phi_r = 0$ for any $\Pi$.
+We use the reference temperature profile ``T_r = T_\text{min} + (T_\text{sfc} - T_\text{min}) \Pi^{n_s}``, with constants ``T_\text{min} = 215\,K``, ``T_\text{sfc}= 288\,K``, and ``n_s = 7``.
 
 #### Horizontal momentum
 
-By breaking the curl and cross product terms into horizontal and vertical contributions, and removing zero terms (e.g. ``\nabla_v  \times \boldsymbol{u}_v = 0``), we obtain
-
+By breaking the curl and cross product terms into horizontal and vertical contributions, and removing zero terms (e.g. ``\nabla_v  \times \boldsymbol{u}_v = 0``), we obtain 
 ```math
 \frac{\partial}{\partial t} \boldsymbol{u}_h  =
   - (2 \boldsymbol{\Omega}^h + \nabla_v \times \boldsymbol{u}_h +  \nabla_h \times \boldsymbol{u}_v) \times \boldsymbol{u}^v
   - (2 \boldsymbol{\Omega}^v + \nabla_h \times \boldsymbol{u}_h) \times \boldsymbol{u}^h
-  - \frac{1}{\rho} \nabla_h (p - p_{\text{ref}})  - \nabla_h (\Phi + K),
+  - c_{pd} (\theta_v - \theta_{v, r}) \nabla_h \Pi  - \nabla_h [(\Phi - \Phi_r) + K],
 ```
-where ``\boldsymbol{u}^h`` and ``\boldsymbol{u}^v`` are the horizontal and vertical _contravariant_ vectors. The effect of topography is accounted for through the computation of the contravariant velocity components (projections from the covariant velocity representation) prior to computing the cross-product contributions.
+where ``\boldsymbol{u}^h`` and ``\boldsymbol{u}^v`` are the horizontal and vertical _contravariant_ vectors. 
+
+The effect of topography is accounted for through the computation of the contravariant velocity components (projections from the covariant velocity representation) prior to computing the cross-product contributions.
 
 This is stabilized with the addition of 4th-order vector hyperviscosity
 ```math
@@ -163,9 +170,9 @@ The ``(2 \boldsymbol{\Omega}^v + \nabla_h \times \boldsymbol{u}_h) \times \bolds
 ```math
 (2 \boldsymbol{\Omega}^v + \mathcal{C}_h[\boldsymbol{u}_h]) \times \boldsymbol{u}^h
 ```
-and the ``\frac{1}{\rho} \nabla_h (p - p_h)  + \nabla_h (\Phi + K)`` as
+and the ``c_{pd} (\theta_v - \theta_{v,r}) \nabla_h \Pi + \nabla_h (\Phi - \Phi_r + K)`` term is discretized as
 ```math
-\frac{1}{\rho} \mathcal{G}_h[p - p_{\text{ref}}] + \mathcal{G}_h[\Phi + K] ,
+c_{pd} (\theta_v - \theta_{v,r}) \mathcal{G}_h[\Pi] + \mathcal{G}_h[\Phi - \Phi_r + K] ,
 ```
 where all these terms are treated explicitly.
 
@@ -183,18 +190,22 @@ Similarly for vertical velocity
 ```math
 \frac{\partial}{\partial t} \boldsymbol{u}_v  =
   - (2 \boldsymbol{\Omega}^h + \nabla_v \times \boldsymbol{u}_h + \nabla_h \times \boldsymbol{u}_v) \times \boldsymbol{u}^h
-  - \frac{1}{\rho} \nabla_v (p - p_{\text{ref}}) - \frac{\rho - \rho_{\text{ref}}}{\rho} \nabla_v \Phi - \nabla_v K .
+  -c_{pd} (\theta_v - \theta_{v, r}) \nabla_v \Pi  - \nabla_v [(\Phi - \Phi_r)].
 ```
 
 The ``(2 \boldsymbol{\Omega}^h + \nabla_v \times \boldsymbol{u}_h + \nabla_h \times \boldsymbol{u}_v) \times \boldsymbol{u}^h`` term is discretized as
 ```math
 (2 \boldsymbol{\Omega}^h + \mathcal{C}^f_v[\boldsymbol{u}_h] + \mathcal{C}_h[\boldsymbol{u}_v]) \times I^f(\boldsymbol{u}^h) ,
 ```
-and the ``\frac{1}{\rho} \nabla_v (p - p_{\text{ref}}) - \frac{\rho - \rho_{\text{ref}}}{\rho} \nabla_v \Phi - \nabla_v K`` term as
+The ``\nabla_v K`` term is discretized as
+```math
+\mathcal{G}^f_v[K],
+```
+The ``c_{pd} (\theta_v - \theta_{v,r}) \nabla_v \Pi + \nabla_v (\Phi - \Phi_r)`` term is discretized as
 ```math
-\frac{1}{I^f(\rho)} \mathcal{G}^f_v[p - p_{\text{ref}}] - \frac{I^f(\rho - \rho_{\text{ref}})}{I^f(\rho)} \mathcal{G}^f_v[\Phi] - \mathcal{G}^f_v[K] ,
+I^f[c_{pd} (\theta_v - \theta_{v, r} ) ] \mathcal{G}^f_v[\Pi] - \mathcal{G}^f_v[\Phi - \Phi_r],
 ```
-with the latter treated implicitly.
+and is treated implicitly.
 
 This is stabilized with the addition of 4th-order vector hyperviscosity
 ```math

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-266
+267
 
 # **README**
 #
@@ -20,6 +20,10 @@
 
 
 #=
+267
+- Changed the pressure-gradient formulation to use the Exner function
+  and virtual potential temperature.
+
 266
 - No behavior change, but some new diagnostics are added, and the reference 
   simulations don't have those diagnostics for plotting

src/ClimaAtmos.jl

@@ -47,6 +47,7 @@ include(joinpath("cache", "surface_albedo.jl"))
 include(joinpath("initial_conditions", "InitialConditions.jl"))
 include(joinpath("surface_conditions", "SurfaceConditions.jl"))
 include(joinpath("utils", "discrete_hydrostatic_balance.jl"))
+include(joinpath("utils", "refstate_thermodynamics.jl"))
 
 include(joinpath("prognostic_equations", "pressure_work.jl"))
 include(joinpath("prognostic_equations", "zero_velocity.jl"))

src/cache/cache.jl

@@ -129,42 +129,6 @@ function build_cache(
     sfc_local_geometry =
         Fields.level(Fields.local_geometry_field(Y.f), Fields.half)
 
-    # Reference enthalpy for hyperdiffusion computation
-    ᶜz = Fields.coordinate_field(Y.c).z
-    T_ref = similar(ᶜz)
-    p_ref = similar(ᶜz)
-    q_tot_ref = similar(ᶜz)
-    cv_m = similar(ᶜz)
-    thermo_params = CAP.thermodynamics_params(params)
-    decay_scale_height = FT(8000)
-
-    # Reference State
-    temp_profile = TD.TemperatureProfiles.DecayingTemperatureProfile{FT}(
-        thermo_params,
-        FT(300),  # surface temperature
-        FT(220),  # top temperature
-        decay_scale_height,  # decay height scale
-    )
-
-    rel_hum_ref =
-        atmos.moisture_model isa DryModel ? FT(0) :
-        @. FT(0.5) * exp(- ᶜz / decay_scale_height)
-
-    T_0 = CAP.T_0(params)
-    grav = CAP.grav(params)
-    @. T_ref = first(temp_profile(thermo_params, ᶜz))
-    @. p_ref = last(temp_profile(thermo_params, ᶜz))
-    @. q_tot_ref =
-        TD.q_vap_from_RH_liquid(thermo_params, p_ref, T_ref, rel_hum_ref)
-    @. cv_m = TD.cv_m(thermo_params, TD.PhasePartition(q_tot_ref, FT(0), FT(0)))
-
-    ᶜh_ref = @. TD.total_specific_enthalpy(
-        thermo_params,
-        cv_m * (T_ref - T_0) + grav * ᶜz,
-        T_ref,
-        TD.PhasePartition(q_tot_ref, FT(0), FT(0)),
-    )
-
     core = (
         ᶜΦ,
         ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(ᶜΦ),
@@ -175,7 +139,6 @@ function build_cache(
         surface_ct3_unit = CT3.(
             unit_basis_vector_data.(CT3, sfc_local_geometry)
         ),
-        ᶜh_ref,
     )
     external_forcing = external_forcing_cache(Y, atmos, params, start_date)
     sfc_setup = surface_setup(params)

src/prognostic_equations/advection.jl

@@ -39,6 +39,7 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t)
     (; ᶜu, ᶜK, ᶜp, ᶜts) = p.precomputed
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)
+    cp_d = thermo_params.cp_d
 
     if p.atmos.turbconv_model isa PrognosticEDMFX
         (; ᶜuʲs) = p.precomputed
@@ -82,7 +83,12 @@ NVTX.@annotate function horizontal_dynamics_tendency!(Yₜ, Y, p, t)
 
     end
 
-    @. Yₜ.c.uₕ -= C12(gradₕ(ᶜp) / Y.c.ρ + gradₕ(ᶜK + ᶜΦ))
+    # This is equivalent to grad_h(Φ + K) + grad_h(p) / ρ
+    ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜts))
+    ᶜθ_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.
     return nothing
 end

src/prognostic_equations/hyperdiffusion.jl

@@ -106,12 +106,16 @@ end
 # dss_hyperdiffusion_tendency_pairs
 NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     (; hyperdiff, turbconv_model) = p.atmos
+    (; params) = p
+    (; ᶜΦ) = p.core
+    (; ᶜts) = p.precomputed
+    thermo_params = CAP.thermodynamics_params(params)
+
     isnothing(hyperdiff) && return nothing
 
     n = n_mass_flux_subdomains(turbconv_model)
     diffuse_tke = use_prognostic_tke(turbconv_model)
     (; ᶜp) = p.precomputed
-    (; ᶜh_ref) = p.core
     (; ᶜ∇²u, ᶜ∇²specific_energy) = p.hyperdiff
     if turbconv_model isa PrognosticEDMFX
         (; ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²mseʲs) = p.hyperdiff
@@ -122,6 +126,8 @@ NVTX.@annotate function prep_hyperdiffusion_tendency!(Yₜ, Y, p, t)
         C123(wgradₕ(divₕ(p.precomputed.ᶜu))) -
         C123(wcurlₕ(C123(curlₕ(p.precomputed.ᶜu))))
 
+    ᶜh_ref = @. lazy(h_dr(thermo_params, ᶜts, ᶜΦ))
+
     @. ᶜ∇²specific_energy =
         wdivₕ(gradₕ(specific(Y.c.ρe_tot, Y.c.ρ) + ᶜp / Y.c.ρ - ᶜh_ref))
 

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -151,13 +151,14 @@ vertical_advection(ᶠu³, ᶜχ, ::Val{:third_order}) =
 function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; moisture_model, turbconv_model, rayleigh_sponge, microphysics_model) =
         p.atmos
-    (; dt) = p
+    (; params, dt) = p
     n = n_mass_flux_subdomains(turbconv_model)
     ᶜJ = Fields.local_geometry_field(axes(Y.c)).J
     ᶠJ = Fields.local_geometry_field(axes(Y.f)).J
     (; ᶠgradᵥ_ᶜΦ) = p.core
     (; ᶠu³, ᶜp, ᶜts) = p.precomputed
-    thermo_params = CAP.thermodynamics_params(p.params)
+    thermo_params = CAP.thermodynamics_params(params)
+    cp_d = CAP.cp_d(params)
     ᶜh_tot = @. lazy(
         TD.total_specific_enthalpy(
             thermo_params,
@@ -245,7 +246,13 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     # TODO - decide if this needs to be explicit or implicit
     #vertical_advection_of_water_tendency!(Yₜ, Y, p, t)
 
-    @. Yₜ.f.u₃ -= ᶠgradᵥ(ᶜp) / ᶠinterp(Y.c.ρ) + ᶠgradᵥ_ᶜΦ
+    # This is equivalent to grad_v(Φ) + grad_v(p) / ρ
+    ᶜΦ_r = @. lazy(phi_r(thermo_params, ᶜts))
+    ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts))
+    ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts))
+    ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts))
+    @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) +
+                  cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ)
 
     if rayleigh_sponge isa RayleighSponge
         ᶠz = Fields.coordinate_field(Y.f).z

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -462,11 +462,20 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
 
     ∂ᶠu₃_err_∂ᶜρ = matrix[@name(f.u₃), @name(c.ρ)]
     ∂ᶠu₃_err_∂ᶜρe_tot = matrix[@name(f.u₃), @name(c.ρe_tot)]
+
+    ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts))
+    ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts))
+    # In implicit tendency, we use the new pressure-gradient formulation (PGF) and gravitational acceleration: 
+    #              grad(p) / ρ + grad(Φ)  =  cp_d * θ_v * grad(Π) + grad(Φ).
+    # Here below, we use the old formulation of (grad(Φ) + grad(p) / ρ).
+    # This is because the new formulation would require computing the derivative of θ_v.
+    # The only exception is:
+    # We are rewriting grad(p) / ρ from the expansion of ∂ᶠu₃_err_∂ᶜρ with the new PGF.
     @. ∂ᶠu₃_err_∂ᶜρ =
         dtγ * (
             ᶠp_grad_matrix ⋅
             DiagonalMatrixRow(ᶜkappa_m * (T_0 * cp_d - ᶜK - ᶜΦ)) +
-            DiagonalMatrixRow(ᶠgradᵥ(ᶜp) / abs2(ᶠinterp(ᶜρ))) ⋅
+            DiagonalMatrixRow(cp_d * ᶠinterp(ᶜθ_v) * ᶠgradᵥ(ᶜΠ) / ᶠinterp(ᶜρ)) ⋅
             ᶠinterp_matrix()
         )
     @. ∂ᶠu₃_err_∂ᶜρe_tot = dtγ * ᶠp_grad_matrix ⋅ DiagonalMatrixRow(ᶜkappa_m)