Increment ref_counter. Add triangle inequality limiter based on Horn 2012 and apply it in microphysics

Author: trontrytel

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-240
+241
 
 # **README**
 #
@@ -20,6 +20,9 @@
 
 
 #=
+241
+- Add triangle inequality limiter from Horn (2012) and apply it in microphysics
+
 240
 - Fix a bug in viscous sponge
 

src/parameterized_tendencies/microphysics/microphysics_wrappers.jl

@@ -29,7 +29,7 @@ function qₚ(q_rain_snow)
     return max(FT(0), q_rain_snow)
 end
 
-# Helper function to limit the tendency
+# Helper function to compute the limit of the tendency
 function limit(q, dt, n::Int)
     return q / float(dt) / n
 end
@@ -85,11 +85,10 @@ function cloud_sources(
 
     S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cm_params, thp, q, ρ, Tₐ(thp, ts))
 
-    # keeping the same limiter for now
     return ifelse(
         S > FT(0),
-        min(S, limit(qᵥ(thp, ts), dt, 2)),
-        -min(abs(S), limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 2)),
+        triangle_inequality_limiter(S, limit(qᵥ(thp, ts), dt, 2)),
+        -triangle_inequality_limiter(abs(S), limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 2)),
     )
 end
 function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, qₛ, dt) where {FT}
@@ -99,11 +98,10 @@ function cloud_sources(cm_params::CMP.CloudIce{FT}, thp, ts, qₛ, dt) where {FT
 
     S = CMNe.conv_q_vap_to_q_liq_ice_MM2015(cm_params, thp, q, ρ, Tₐ(thp, ts))
 
-    # keeping the same limiter for now
     return ifelse(
         S > FT(0),
-        min(S, limit(qᵥ(thp, ts), dt, 2)),
-        -min(abs(S), limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 2)),
+        triangle_inequality_limiter(S, limit(qᵥ(thp, ts), dt, 2)),
+        -triangle_inequality_limiter(abs(S), limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 2)),
     )
 end
 
@@ -119,9 +117,9 @@ Returns the qₜ source term due to precipitation formation
 defined as Δm_tot / (m_dry + m_tot) for the 0-moment scheme
 """
 function q_tot_0M_precipitation_sources(thp, cmp::CMP.Parameters0M, dt, qₜ, ts)
-    return -min(
-        max(qₜ, 0) / float(dt),
+    return -triangle_inequality_limiter(
         -CM0.remove_precipitation(cmp, PP(thp, ts)),
+        max(qₜ, 0) / float(dt),
     )
 end
 
@@ -190,77 +188,77 @@ function compute_precipitation_sources!(
         CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts, qₚ(qᵣ)), true),
         CM2.conv_q_liq_to_q_rai(mp.var, qₗ(thp, ts, qₚ(qᵣ)), ρ, mp.Ndp),
     )
-    @. Sᵖ = min(limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5), Sᵖ)
+    @. Sᵖ = triangle_inequality_limiter(Sᵖ, limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5))
     @. Sqₗᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
 
     # snow autoconversion assuming no supersaturation: q_ice -> q_snow
-    @. Sᵖ = min(
-        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts, qₚ(qₛ)), true),
+        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
     )
     @. Sqᵢᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
 
     # accretion: q_liq + q_rain -> q_rain
-    @. Sᵖ = min(
-        limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts, qₚ(qᵣ)), qₚ(qᵣ), ρ),
+        limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
     )
     @. Sqₗᵖ -= Sᵖ
     @. Sqᵣᵖ += Sᵖ
 
     # accretion: q_ice + q_snow -> q_snow
-    @. Sᵖ = min(
-        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qₛ), ρ),
+        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
     )
     @. Sqᵢᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
 
     # accretion: q_liq + q_sno -> q_sno or q_rai
     # sink of cloud water via accretion cloud water + snow
-    @. Sᵖ = min(
-        limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts, qₚ(qᵣ)), qₚ(qₛ), ρ),
+        limit(qₗ(thp, ts, qₚ(qᵣ)), dt, 5),
     )
     # if T < T_freeze cloud droplets freeze to become snow
     # else the snow melts and both cloud water and snow become rain
     α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze)
     @. Sᵖ_snow = ifelse(
         Tₐ(thp, ts) < mp.ps.T_freeze,
         Sᵖ,
-        FT(-1) * min(Sᵖ * α(thp, ts), limit(qₚ(qₛ), dt, 5)),
+        FT(-1) * triangle_inequality_limiter(Sᵖ * α(thp, ts), limit(qₚ(qₛ), dt, 5)),
     )
     @. Sqₛᵖ += Sᵖ_snow
     @. Sqₗᵖ -= Sᵖ
     @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow)
 
     # accretion: q_ice + q_rai -> q_sno
-    @. Sᵖ = min(
-        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qᵣ), ρ),
+        limit(qᵢ(thp, ts, qₚ(qₛ)), dt, 5),
     )
     @. Sqᵢᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
     # sink of rain via accretion cloud ice - rain
-    @. Sᵖ = min(
-        limit(qₚ(qᵣ), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₚ(qₛ)), qₚ(qᵣ), ρ),
+        limit(qₚ(qᵣ), dt, 5),
     )
     @. Sqᵣᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
 
     # accretion: q_rai + q_sno -> q_rai or q_sno
     @. Sᵖ = ifelse(
         Tₐ(thp, ts) < mp.ps.T_freeze,
-        min(
-            limit(qₚ(qᵣ), dt, 5),
+        triangle_inequality_limiter(
             CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₚ(qₛ), qₚ(qᵣ), ρ),
+            limit(qₚ(qᵣ), dt, 5),
         ),
-        -min(
-            limit(qₚ(qₛ), dt, 5),
+        -triangle_inequality_limiter(
             CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qₚ(qᵣ), qₚ(qₛ), ρ),
+            limit(qₚ(qₛ), dt, 5),
         ),
     )
     @. Sqₛᵖ += Sᵖ
@@ -302,16 +300,16 @@ function compute_precipitation_sinks!(
 
     #! format: off
     # evaporation: q_rai -> q_vap
-    @. Sᵖ = -min(
-        limit(qₚ(qᵣ), dt, 5),
+    @. Sᵖ = -triangle_inequality_limiter(
         -CM1.evaporation_sublimation(rps..., PP(thp, ts), qₚ(qᵣ), ρ, Tₐ(thp, ts)),
+        limit(qₚ(qᵣ), dt, 5),
     )
     @. Sqᵣᵖ += Sᵖ
 
     # melting: q_sno -> q_rai
-    @. Sᵖ = min(
-        limit(qₚ(qₛ), dt, 5),
+    @. Sᵖ = triangle_inequality_limiter(
         CM1.snow_melt(sps..., qₚ(qₛ), ρ, Tₐ(thp, ts)),
+        limit(qₚ(qₛ), dt, 5),
     )
     @. Sqᵣᵖ += Sᵖ
     @. Sqₛᵖ -= Sᵖ
@@ -320,8 +318,8 @@ function compute_precipitation_sinks!(
     @. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₚ(qₛ), ρ, Tₐ(thp, ts))
     @. Sᵖ = ifelse(
         Sᵖ > FT(0),
-        min(limit(qᵥ(thp, ts), dt, 5), Sᵖ),
-        -min(limit(qₚ(qₛ), dt, 5), FT(-1) * Sᵖ),
+        triangle_inequality_limiter(Sᵖ, limit(qᵥ(thp, ts), dt, 5)),
+        -triangle_inequality_limiter(FT(-1) * Sᵖ, limit(qₚ(qₛ), dt, 5)),
     )
     @. Sqₛᵖ += Sᵖ
     #! format: on

src/prognostic_equations/limited_tendencies.jl

@@ -1,4 +1,5 @@
 
+# Interacting with ClimaCore quasimonotone limiter for advection
 NVTX.@annotate function limiters_func!(Y, p, t, ref_Y)
     (; limiter) = p.numerics
     if !isnothing(limiter)
@@ -9,3 +10,9 @@ NVTX.@annotate function limiters_func!(Y, p, t, ref_Y)
     end
     return nothing
 end
+
+# Helper function to limit source term tendencies the trianlge inequality-based
+# limiter from Horn (2012, doi:10.5194/gmd-5-345-2012)
+function triangle_inequality_limiter(force, limit)
+    return force + limit - sqrt(force^2 + limit^2)
+end