Add limiter option for NonEq in parcel

Author: oalcabes

parcel/Example_NonEq.jl

@@ -11,6 +11,7 @@ include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 
 # choosing a value for the ice relaxation timescale
 τ = FT(10)
+limiter = true
 @info("using τ value ", τ)
 override_file = Dict(
     "condensation_evaporation_timescale" =>
@@ -72,6 +73,7 @@ ax6 = MK.Axis(fig[3, 2], xlabel = "Time [s]", ylabel = "q_ice [g/kg]")
 params = parcel_params{FT}(
     condensation_growth = condensation_growth,
     deposition_growth = deposition_growth,
+    limiter = limiter,
     const_dt = const_dt,
     w = w,
     liquid = liquid,

parcel/ParcelModel.jl

@@ -21,6 +21,7 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     deposition_growth = "None"
     liq_size_distribution = "Monodisperse"
     ice_size_distribution = "Monodisperse"
+    limiter = false
     aerosol = Empty{FT}()
     aero_σ_g = FT(0)
     wps = CMP.WaterProperties(FT)
@@ -293,7 +294,7 @@ function run_parcel(IC, t_0, t_end, pp)
     elseif pp.condensation_growth == "NonEq_Condensation_simple"
         ce_params = NonEqCondParams_simple{FT}(pp.tps, pp.liquid)
     elseif pp.condensation_growth == "NonEq_Condensation"
-        ce_params = NonEqCondParams{FT}(pp.tps, pp.liquid, pp.const_dt)
+        ce_params = NonEqCondParams{FT}(pp.tps, pp.liquid, pp.limiter, pp.const_dt)
     else
         throw("Unrecognized condensation growth mode")
     end
@@ -306,7 +307,7 @@ function run_parcel(IC, t_0, t_end, pp)
     elseif pp.deposition_growth == "NonEq_Deposition_simple"
         ds_params = NonEqDepParams_simple{FT}(pp.tps, pp.ice)
     elseif pp.deposition_growth == "NonEq_Deposition"
-        ds_params = NonEqDepParams{FT}(pp.tps, pp.ice, pp.const_dt)
+        ds_params = NonEqDepParams{FT}(pp.tps, pp.ice, pp.limiter, pp.const_dt)
     else
         throw("Unrecognized deposition growth mode")
     end

parcel/ParcelParameters.jl

@@ -90,6 +90,7 @@ end
 struct NonEqCondParams{FT} <: CMP.ParametersType{FT}
     tps::TDP.ThermodynamicsParameters{FT}
     liquid::CMP.CloudLiquid{FT}
+    limiter::Bool
     dt::FT
 end
 
@@ -107,5 +108,6 @@ end
 struct NonEqDepParams{FT} <: CMP.ParametersType{FT}
     tps::TDP.ThermodynamicsParameters{FT}
     ice::CMP.CloudIce{FT}
+    limiter::Bool
     dt::FT
 end

parcel/ParcelTendencies.jl

@@ -12,6 +12,11 @@ function limit(q, dt, n::Int)
     return q / dt / n
 end
 
+# the triangle inequality limiter used in ClimaAtmos
+function triangle_inequality_limiter(force, limit)
+    return force + limit - sqrt(force^2 + limit^2)
+end
+
 function aerosol_activation(::Empty, state)
     FT = eltype(state)
     return FT(0)
@@ -249,23 +254,28 @@ end
 function condensation(params::NonEqCondParams, PSD, state, ρ_air)
     FT = eltype(state)
     (; T, qₗ, qᵥ, qᵢ) = state
-    (; tps, liquid, dt) = params
 
-    if qᵥ + qₗ > FT(0)
+    (; tps, liquid, limiter, dt) = params
 
-        qₜ = qᵥ + qₗ + qᵢ
+    qₜ = qᵥ + qₗ + qᵢ
 
-        cond_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
+    cond_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
 
-        # Using same limiter as ClimaAtmos for now
+    # Using same limiter as ClimaAtmos
+    if limiter
         # Not sure why, but without intermediate storing of the tendencies for the
         # if/else branch this code segfaults on julia v1.11 (works fine on v1.10)
+
         cond_limit = min(cond_rate, limit(qᵥ, dt, 1))
         evap_limit = min(abs(cond_rate), limit(qₗ, dt, 1))
-        ret = ifelse(cond_rate > FT(0), cond_limit, -1 * evap_limit)
-        return ret
+
+        return ifelse(
+            cond_rate > FT(0),
+            triangle_inequality_limiter(cond_rate, cond_limit),
+            -triangle_inequality_limiter(abs(cond_rate), evap_limit),
+        )
     else
-        return FT(0)
+        return cond_rate
     end
 end
 
@@ -307,21 +317,24 @@ function deposition(params::NonEqDepParams, PSD, state, ρ_air)
     FT = eltype(state)
     (; T, qₗ, qᵥ, qᵢ) = state
 
-    (; tps, ice, dt) = params
+    (; tps, ice, limiter, dt) = params
 
-    if qᵥ + qᵢ > FT(0)
-        qₜ = qᵥ + qₗ + qᵢ
+    qₜ = qᵥ + qₗ + qᵢ
 
-        dep_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
+    dep_rate = MNE.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, qₜ, qₗ, qᵢ, FT(0), FT(0), ρ_air, T)
 
-        # Using same limiter as ClimaAtmos for now
+    # using same limiter as ClimaAtmos for now
+    if limiter
         # Not sure why, but without intermediate storing of the tendencies for the
         # if/else branch this code segfaults on julia v1.11 (works fine on v1.10)
         dep_limit = min(dep_rate, limit(qᵥ, dt, 1))
         sub_limit = min(abs(dep_rate), limit(qᵢ, dt, 1))
-        ret = ifelse(dep_rate > FT(0), dep_limit, -1 * sub_limit)
-        return ret
+        return ifelse(
+            dep_rate > FT(0),
+            triangle_inequality_limiter(dep_rate, dep_limit),
+            -triangle_inequality_limiter(abs(dep_rate), sub_limit),
+        )
     else
-        return FT(0)
+        return dep_rate
     end
 end