Merge pull request #649 from CliMA/aj/new_epsilon

Change to sqrt(ftmin)

Author: trontrytel

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.28.1"
+version = "0.29.0"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

src/CloudDiagnostics.jl

@@ -12,6 +12,7 @@ import ..Parameters as CMP
 import ..Microphysics1M as CM1
 import ..Microphysics2M as CM2
 import ..DistributionTools as DT
+import ..Common as CO
 
 """
     radar_reflectivity_1M(precip, q, ρ)
@@ -119,7 +120,7 @@ function effective_radius_2M(
     M2_c = notvalid(Bc) ? FT(0) : DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_lcl, n_mass) / C^(n_mass)
     M2_r = notvalid(Br) ? FT(0) : DT.generalized_gamma_Mⁿ(νr, μr, Br, N_rai, n_mass) / C^(n_mass)
 
-    return M2_c + M2_r <= eps(FT) ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r)
+    return M2_c + M2_r <= CO.ϵ_numerics(FT) ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r)
 end
 
 """
@@ -147,7 +148,7 @@ function effective_radius_Liu_Hallet_97(
 
     k = FT(0.8)
     r_vol =
-        ((N_lcl + N_rai) < eps(FT)) ? FT(0) :
+        ((N_lcl + N_rai) < CO.ϵ_numerics(FT)) ? FT(0) :
         (
             (FT(3) * (q_lcl + q_rai) * ρ_air) /
             (FT(4) * π * ρw * (N_lcl + N_rai))

src/Common.jl

@@ -21,6 +21,12 @@ export Chen2022_monodisperse_pdf
 export Chen2022_exponential_pdf
 export ventilation_factor
 
+"""
+    Smallest number that is different than zero for the purpose of microphysics
+    computations.
+"""
+ϵ_numerics(FT) = sqrt(floatmin(FT))
+
 """
     G_func_liquid(air_props, tps, T)
 
@@ -89,9 +95,9 @@ function logistic_function(x::FT, x_0::FT, k::FT) where {FT}
     @assert x_0 >= 0
     x = max(0, x)
 
-    if abs(x) < eps(FT)
+    if abs(x) < ϵ_numerics(FT)
         return FT(0)
-    elseif abs(x_0) < eps(FT)
+    elseif abs(x_0) < ϵ_numerics(FT)
         return FT(1)
     end
 
@@ -114,9 +120,9 @@ function logistic_function_integral(x::FT, x_0::FT, k::FT) where {FT}
     @assert x_0 >= 0
     x = max(0, x)
 
-    if abs(x) < eps(FT)
+    if abs(x) < ϵ_numerics(FT)
         return FT(0)
-    elseif abs(x_0) < eps(FT)
+    elseif abs(x_0) < ϵ_numerics(FT)
         return x
     end
 

src/Microphysics1M.jl

@@ -94,7 +94,7 @@ function lambda_inverse(
     (; r0, m0, me, Δm, χm) = mass
 
     λ_inv = FT(0)
-    if q > FT(0) && ρ > FT(0)
+    if q > CO.ϵ_numerics(FT) && ρ > CO.ϵ_numerics(FT)
         λ_inv = exp(1 / (me + Δm + 1) *
             log(
             ρ * q * exp((me + Δm) * log(r0)) / χm / m0 / n0 / SF.gamma(me + Δm + FT(1)),
@@ -162,7 +162,7 @@ function terminal_velocity(
     ρ::FT,
     q::FT,
 ) where {FT}
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         # terminal_velocity(size)
         (; χv, ve, Δv) = vel
         v0 = get_v0(vel, ρ)
@@ -184,7 +184,7 @@ function terminal_velocity(
     q::FT,
 ) where {FT}
     fall_w = FT(0)
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         # coefficients from Table B1 from Chen et. al. 2022
         aiu, bi, ciu = CO.Chen2022_vel_coeffs(vel, ρₐ)
         # size distribution parameter
@@ -208,7 +208,7 @@ function terminal_velocity(
     # We assume the B4 table coeffs for snow and B2 table coeffs for cloud ice.
     # Instead we should do partial integrals
     # from D=125um to D=625um using B2 and D=625um to inf using B4.
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         # coefficients from Table B4 from Chen et. al. 2022
         aiu, bi, ciu = CO.Chen2022_vel_coeffs(vel, ρₐ, ρᵢ)
         # size distribution parameter
@@ -235,7 +235,7 @@ function terminal_velocity(
 ) where {FT}
     fall_w = FT(0)
     # see comments above about B2 vs B4 coefficients
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         # coefficients from Table B4 from Chen et. al. 2022
         aiu, bi, ciu = CO.Chen2022_vel_coeffs(vel, ρₐ, ρᵢ)
         # size distribution parameter
@@ -325,7 +325,7 @@ function conv_q_icl_to_q_sno(
     acnv_rate = FT(0)
     S = TDI.supersaturation_over_ice(tps, q_tot, q_lcl + q_rai, q_icl + q_sno, ρ, T)
 
-    if (q_icl > FT(0) && S > FT(0))
+    if (q_icl > CO.ϵ_numerics(FT) && S > FT(0))
         (; me, Δm) = mass
         G = CO.G_func_ice(aps, tps, T)
         n0 = get_n0(pdf)
@@ -364,7 +364,7 @@ function accretion(
 ) where {FT}
 
     accr_rate = FT(0)
-    if (q_clo > FT(0) && q_pre > FT(0))
+    if (q_clo > CO.ϵ_numerics(FT) && q_pre > CO.ϵ_numerics(FT))
 
         n0::FT = get_n0(precip.pdf, q_pre, ρ)
         v0::FT = get_v0(vel, ρ)
@@ -407,7 +407,7 @@ function accretion_rain_sink(
     ρ::FT,
 ) where {FT}
     accr_rate = FT(0)
-    if (q_icl > FT(0) && q_rai > FT(0))
+    if (q_icl > CO.ϵ_numerics(FT) && q_rai > CO.ϵ_numerics(FT))
 
         n0_ice = get_n0(ice.pdf)
         λ_ice_inv = lambda_inverse(ice.pdf, ice.mass, q_icl, ρ)
@@ -459,7 +459,7 @@ function accretion_snow_rain(
 ) where {FT}
 
     accr_rate = FT(0)
-    if (q_i > FT(0) && q_j > FT(0))
+    if (q_i > CO.ϵ_numerics(FT) && q_j > CO.ϵ_numerics(FT))
 
         n0_i = get_n0(type_i.pdf, q_i, ρ)
         n0_j = get_n0(type_j.pdf, q_j, ρ)
@@ -524,7 +524,7 @@ function evaporation_sublimation(
     evap_subl_rate = FT(0)
     S = TDI.supersaturation_over_liquid(tps, q_tot, q_lcl + q_rai, q_icl + q_sno, ρ, T)
 
-    if (q_rai > FT(0) && S < FT(0))
+    if (q_rai > CO.ϵ_numerics(FT) && S < FT(0))
 
         (; ν_air, D_vapor) = aps
         G = CO.G_func_liquid(aps, tps, T)
@@ -565,7 +565,7 @@ function evaporation_sublimation(
     T::FT,
 ) where {FT}
     evap_subl_rate = FT(0)
-    if q_sno > FT(0)
+    if q_sno > CO.ϵ_numerics(FT)
         (; ν_air, D_vapor) = aps
 
         S = TDI.supersaturation_over_ice(tps, q_tot, q_lcl + q_rai, q_icl + q_sno, ρ, T)
@@ -618,7 +618,7 @@ function snow_melt(
 ) where {FT}
     snow_melt_rate = FT(0)
 
-    if (q_sno > FT(0) && T > T_freeze)
+    if (q_sno > CO.ϵ_numerics(FT) && T > T_freeze)
         (; ν_air, D_vapor, K_therm) = aps
 
         L = TDI.Lf(tps, T)

src/Microphysics2M.jl

@@ -365,7 +365,7 @@ function autoconversion(
     q_lcl, q_rai, ρ, N_lcl,
 ) where {FT}
 
-    if q_lcl < eps(FT) || N_lcl < eps(FT)
+    if q_lcl < CO.ϵ_numerics(FT) || N_lcl < CO.ϵ_numerics(FT)
         return LclRaiRates{FT}()
     end
 
@@ -411,7 +411,7 @@ Compute accretion rate
 """
 function accretion((; accr)::CMP.SB2006{FT}, q_lcl, q_rai, ρ, N_lcl) where {FT}
 
-    if q_lcl < eps(FT) || q_rai < eps(FT) || N_lcl < eps(FT)
+    if q_lcl < CO.ϵ_numerics(FT) || q_rai < CO.ϵ_numerics(FT) || N_lcl < CO.ϵ_numerics(FT)
         return LclRaiRates{FT}()
     end
 
@@ -456,7 +456,7 @@ function cloud_liquid_self_collection(
     q_lcl, ρ, dN_lcl_dt_au,
 ) where {FT}
 
-    if q_lcl < eps(FT)
+    if q_lcl < CO.ϵ_numerics(FT)
         return FT(0)
     end
     (; kcc, ρ0) = acnv
@@ -518,7 +518,7 @@ function rain_self_collection(
     q_rai, ρ, N_rai,
 ) where {FT}
 
-    if q_rai < eps(FT) || N_rai < eps(FT)
+    if q_rai < CO.ϵ_numerics(FT) || N_rai < CO.ϵ_numerics(FT)
         return FT(0)
     end
 
@@ -554,7 +554,7 @@ function rain_breakup(
     q_rai, ρ, N_rai, dN_rai_dt_sc,
 ) where {FT}
 
-    if q_rai < eps(FT) || N_rai < eps(FT)
+    if q_rai < CO.ϵ_numerics(FT) || N_rai < CO.ϵ_numerics(FT)
         return FT(0)
     end
     (; Deq, Dr_th, kbr, κbr) = brek
@@ -624,7 +624,7 @@ function cloud_terminal_velocity(
     q_liq, ρₐ, N_liq,
 ) where {FT}
 
-    if N_liq < eps(FT) || q_liq < eps(FT)
+    if N_liq < CO.ϵ_numerics(FT) || q_liq < CO.ϵ_numerics(FT)
         return (FT(0), FT(0))
     end
 
@@ -633,10 +633,10 @@ function cloud_terminal_velocity(
 
     terminal_velocity_prefactor = FT(1 / 18) * (6 / ρw / π)^(2 // 3) * (ρw / ρₐ - 1) * grav / ν_air
     vt0 =
-        N_liq < eps(FT) ? FT(0) :
+        N_liq < CO.ϵ_numerics(FT) ? FT(0) :
         terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(2 / 3)) / N_liq
     vt1 =
-        q_liq < eps(FT) ? FT(0) :
+        q_liq < CO.ϵ_numerics(FT) ? FT(0) :
         terminal_velocity_prefactor * DT.generalized_gamma_Mⁿ(νc, μc, Bc, N_liq, FT(5 / 3)) / ρₐ / q_liq
 
     return (vt0, vt1)
@@ -673,10 +673,10 @@ function rain_terminal_velocity(
         _sb_rain_terminal_velocity_helper(pdf_r, 1 / Dr_mean, aR, bR, cR)
 
     vt0 =
-        N_rai < eps(FT) ? FT(0) :
+        N_rai < CO.ϵ_numerics(FT) ? FT(0) :
         max(FT(0), sqrt(ρ0 / ρ) * (aR * _pa0 - bR * _pb0 / (1 + cR * Dr_mean)))
     vt1 =
-        q_rai < eps(FT) ? FT(0) :
+        q_rai < CO.ϵ_numerics(FT) ? FT(0) :
         max(FT(0), sqrt(ρ0 / ρ) * (aR * _pa1 - bR * _pb1 / (1 + cR * Dr_mean)^4))
     return (vt0, vt1)
 end
@@ -757,7 +757,7 @@ function rain_evaporation(
     evap_rate_1 = FT(0)
     S = TDI.supersaturation_over_liquid(tps, q_tot, q_lcl + q_rai, q_icl + q_sno, ρ, T)
 
-    if (N_rai > eps(FT) && S < FT(0))
+    if (N_rai > CO.ϵ_numerics(FT) && S < FT(0))
 
         (; ν_air, D_vapor) = aps
         (; av, bv, α, β, ρ0) = evap
@@ -917,7 +917,7 @@ function conv_q_lcl_to_q_rai(
     N_d,
     smooth_transition = false,
 ) where {FT}
-    if q_lcl <= eps(FT)
+    if q_lcl <= CO.ϵ_numerics(FT)
         return FT(0)
     else
         # Mean volume radius in microns (assuming spherical cloud droplets)

src/MicrophysicsNonEq.jl

@@ -138,7 +138,7 @@ function terminal_velocity(
     q::FT,
 ) where {FT}
     fall_w = FT(0)
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         # TODO: Coefficients from Table B1 from Chen et. al. 2022 are only valid
         # for D > 100mm. We should look for a different parameterization
         # that is more suited for cloud droplets. For now I'm just multiplying
@@ -162,7 +162,7 @@ function terminal_velocity(
     q::FT,
 ) where {FT}
     fall_w = FT(0)
-    if q > FT(0)
+    if q > CO.ϵ_numerics(FT)
         v_term = CO.particle_terminal_velocity(vel, ρₐ, ρᵢ)
         # See the comment for liquid droplets above
         N = FT(500 * 1e6)

src/P3_terminal_velocity.jl

@@ -3,7 +3,7 @@
     ice_particle_terminal_velocity(velocity_params, ρₐ, state::P3State; [use_aspect_ratio])
     ice_particle_terminal_velocity(velocity_params, ρₐ, params::CMP.ParametersP3, F_rim, ρ_rim; [use_aspect_ratio])
 
-Returns the terminal velocity of a single ice particle as a function of its size 
+Returns the terminal velocity of a single ice particle as a function of its size
     (maximum dimension, `D`) using the Chen 2022 parametrization.
 
 # Arguments
@@ -57,6 +57,7 @@ function ice_terminal_velocity_number_weighted(
     use_aspect_ratio = true, p = 1e-6, ∫kwargs...,
 )
     (; N_ice, L_ice) = state
+    # TODO - do we want to swicth to ϵ_numerics(FT)
     if N_ice < eps(one(N_ice)) || L_ice < eps(one(L_ice))
         return zero(N_ice)
     end
@@ -101,6 +102,7 @@ function ice_terminal_velocity_mass_weighted(
     use_aspect_ratio = true, p = 1e-6, ∫kwargs...,
 )
     (; N_ice, L_ice) = state
+    # TODO - do we want to swicth to ϵ_numerics(FT)
     if N_ice < eps(one(N_ice)) || L_ice < eps(one(L_ice))
         return zero(L_ice)
     end

test/cloud_diagnostics.jl

@@ -46,27 +46,40 @@ function test_cloud_diagnostics(FT)
         #setup
         ρₐ = FT(1)
 
-        q_lcl = [FT(2.128e-4), FT(2.128e-20), FT(1.6e-12), FT(0), FT(1.037e-25)]
-        N_lcl = [FT(15053529), FT(3), FT(5512), FT(0), FT(5.225e-12)]
-        q_rai = [FT(1.573e-4), FT(1.573e-4), FT(1.9e-15), FT(0), FT(2.448e-27)]
-        N_rai = [FT(510859), FT(510859), FT(0), FT(0), FT(5.136e-18)]
+        q_lcl = [FT(2.128e-4), FT(2.128e-20), FT(1.6e-12), FT(0)]
+        N_lcl = [FT(15053529), FT(3), FT(5512), FT(0)]
+        q_rai = [FT(1.573e-4), FT(1.573e-4), FT(1.9e-15), FT(0)]
+        N_rai = [FT(510859), FT(510859), FT(0), FT(0)]
 
         # reference values
-        rr = [FT(-12.559725319858543), FT(-12.579899), FT(-150), FT(-150), FT(-150)]
-        reff = [FT(2.319383e-5), FT(6.91594e-5), FT(0), FT(0), FT(0)]
+        rr = [FT(-12.559725319858543), FT(-12.579899), FT(-150), FT(-150)]
+        reff = [FT(2.319383e-5), FT(6.91594e-5), FT(0), FT(0)]
 
         for (qₗ, Nₗ, qᵣ, Nᵣ, rₑ, Z) in zip(q_lcl, N_lcl, q_rai, N_rai, reff, rr)
             for SB in [SB2006, SB2006_no_limiters]
-
                 #action
                 Z_val = CMD.radar_reflectivity_2M(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
                 rₑ_val = CMD.effective_radius_2M(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
-
                 #test
                 TT.@test rₑ_val ≈ rₑ atol = FT(1e-6)
                 TT.@test Z_val ≈ Z atol = FT(1e-4)
             end
         end
+
+        # Additional test for small numbers
+        qₗ = FT(1.037e-25)
+        Nₗ = FT(5.225e-12)
+        qᵣ = FT(2.448e-27)
+        Nᵣ = FT(5.136e-18)
+        Z = FT(-150)
+        for SB in [SB2006, SB2006_no_limiters]
+            rₑ = (FT == Float32 || SB == SB2006) ? FT(0) : FT(8e-5)
+            Z_val = CMD.radar_reflectivity_2M(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
+            rₑ_val = CMD.effective_radius_2M(SB, qₗ, qᵣ, Nₗ, Nᵣ, ρₐ)
+            #test
+            TT.@test rₑ_val ≈ rₑ atol = FT(1e-6)
+            TT.@test Z_val ≈ Z atol = FT(1e-4)
+        end
     end
 
     TT.@testset "Effective radius - '1/3' power law from Liu and Hallett (1997)" begin

test/microphysics1M_tests.jl

@@ -82,6 +82,9 @@ function test_microphysics1M(FT)
         TT.@test v_bigger > vt_sno
     end
 
+    TT.@testset "1M_microphysics - 1M snow terminal velocity Nans" begin
+        TT.@test !isnan(CM1.terminal_velocity(snow, blk1mvel.snow, FT(0.2439843), FT(3.0f-45)))
+    end
 
     TT.@testset "RainAutoconversion" begin