fix doc error + DistributionTools

Author: haakon-e

docs/src/plots/RainEvapoartionSB2006.jl

@@ -30,10 +30,10 @@ function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T)
         x_star = SB2006.pdf_r.xr_min
         G = CO.G_func(aps, tps, T, TD.Liquid())
 
-        xr = CM2.pdf_rain_parameters(SB2006.pdf_r, q_rai, ρ, N_rai).xr
-        Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr^FT(1 / 3)
+        (; xr_mean) = CM2.pdf_rain_parameters(SB2006.pdf_r, q_rai, ρ, N_rai)
+        Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr_mean^FT(1 / 3)
 
-        t_star = (FT(6) * x_star / xr)^FT(1 / 3)
+        t_star = (FT(6) * x_star / xr_mean)^FT(1 / 3)
         a_vent_0 = av * FT(SF.gamma(-1, t_star)) / FT(6)^FT(-2 / 3)
         b_vent_0 =
             bv * FT(SF.gamma((-1 / 2) + (3 / 2) * β, t_star)) /
@@ -43,11 +43,11 @@ function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T)
         b_vent_1 =
             bv * SF.gamma(FT(5 / 2) + FT(3 / 2) * β) / FT(6)^FT(β / 2 + 1 / 2)
 
-        N_Re = α * xr^β * sqrt(ρ0 / ρ) * Dr / ν_air
+        N_Re = α * xr_mean^β * sqrt(ρ0 / ρ) * Dr / ν_air
         Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
         Fv1 = a_vent_1 + b_vent_1 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re)
 
-        evap_rate_0 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv0 / xr)
+        evap_rate_0 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv0 / xr_mean)
         evap_rate_1 = min(FT(0), FT(2) * FT(π) * G * S * N_rai * Dr * Fv1 / ρ)
     end
 

src/DistributionTools.jl

@@ -11,6 +11,7 @@ Mostly used for the 2-moment microphysics.
 module DistributionTools
 
 import SpecialFunctions as SF
+import LogExpFunctions as LEF
 
 """
     generalized_gamma_quantile(ν, μ, B, Y)
@@ -58,7 +59,7 @@ function generalized_gamma_cdf(ν, μ, B, x)
     x ≤ 0 && return zero(x)
 
     # Compute the regularized incomplete gamma function
-    p, _ = SF.gamma_inc((ν + 1) / μ, (B * x)^μ)
+    p, _ = SF.gamma_inc((ν + 1) / μ, B * x^μ)
     return p
 end
 
@@ -110,7 +111,8 @@ function exponential_cdf(D_mean, D)
     # Handle edge cases
     D < 0 && return zero(D)
     # Calculate CDF: P(X ≤ D) = 1 - exp(-D/D_mean)
-    return 1 - exp(-D / D_mean)
+    logcdf = LEF.log1mexp(-D / D_mean)  # careful calculation in log-space
+    return exp(logcdf)
 end
 
 """
@@ -135,7 +137,8 @@ function exponential_quantile(D_mean, Y)
     (0 ≤ Y ≤ 1) || throw(DomainError(Y, "Probability Y must be in [0,1]"))
     D_mean > 0 || throw(DomainError(D_mean, "Mean parameter must be positive"))
     # Calculate quantile: x = -D_mean * ln(1-Y)
-    return -D_mean * log(1 - Y)
+    logquantile = log(D_mean) + LEF.cloglog(Y)  # careful calculation in log-space
+    return exp(logquantile)
 end
 
 """

src/Microphysics2M.jl

@@ -81,7 +81,7 @@ function pdf_rain_parameters(pdf_r::CMP.RainParticlePDF_SB2006_limited, qᵣ, ρ
     # Sequence of limiting steps in Seifert and Beheng 2006:
     x̃r = clamp(Lᵣ / Nᵣ, xr_min, xr_max)  # Eq. (94)  # TODO: Ill-defined for 0 / 0
     N₀r = clamp(Nᵣ * ∛(π * ρw / x̃r), N0_min, N0_max)  # Eq. (95)
-    λr = clamp((π * ρw * N₀r / Lᵣ)^(1 // 4), λ_min, λ_max)  # Eq. (96)
+    λr = clamp(∜(π * ρw * N₀r / Lᵣ), λ_min, λ_max)  # Eq. (96)
     xr_mean = clamp(Lᵣ * λr / N₀r, xr_min, xr_max)  # Eq. (97)
 
     Dr_mean = 1 / λr  # The inverse of λr is the mean diameter of the raindrops (units: `m`)

test/microphysics2M_tests.jl

@@ -535,32 +535,32 @@ function test_microphysics2M(FT)
             Mⁿ(n, psd) = y -> y^n * psd(y)
 
             # integral bounds computed based on the size distribution
-            p = eps(FT)
+            p = FT(1e-6)
             D_min, D_max = CM2.get_size_distribution_bounds(pdf_r, ρₐ, qᵣ, Nᵣ, p)
             x_min = DT.generalized_gamma_quantile(νr, μr, Br, p)
             x_max = DT.generalized_gamma_quantile(νr, μr, Br, 1 - p)
 
             # Test that these bounds correspond to the correct probability levels
-            TT.@test_broken DT.generalized_gamma_cdf(νr, μr, Br, x_min) ≈ p # it's a bit unstable at very low values
-            TT.@test_broken DT.generalized_gamma_cdf(νr, μr, Br, x_max) ≈ 1 - p  # TODO: fix this -- not clear why this doesn't work, seems like an issue with SpecialFunctions.jl??
+            TT.@test DT.generalized_gamma_cdf(νr, μr, Br, x_min) ≈ p
+            TT.@test DT.generalized_gamma_cdf(νr, μr, Br, x_max) ≈ 1 - p
             TT.@test DT.exponential_cdf(Dr_mean, D_min) ≈ p
             TT.@test DT.exponential_cdf(Dr_mean, D_max) ≈ 1 - p
 
             # Sanity checks for number concentrations for rain
             ND = QGK.quadgk(f_D, D_min, D_max)[1]
             Nx = QGK.quadgk(f_x, x_min, x_max)[1]
             ND_psd = QGK.quadgk(psd, D_min, D_max)[1]
-            TT.@test ND ≈ Nᵣ
-            TT.@test Nx ≈ Nᵣ
-            TT.@test ND_psd ≈ Nᵣ
+            TT.@test ND ≈ Nᵣ rtol = 1e-5
+            TT.@test Nx ≈ Nᵣ rtol = 5e-3
+            TT.@test ND_psd ≈ Nᵣ rtol = 1e-5
 
             # Sanity checks for specific contents for rain
             qD = QGK.quadgk(Mⁿ(3, f_D), D_min, D_max)[1] * k_m / ρₐ
             qx = QGK.quadgk(Mⁿ(1, f_x), x_min, x_max)[1] / ρₐ
             qD_psd = QGK.quadgk(Mⁿ(3, psd), D_min, D_max)[1] * k_m / ρₐ
-            TT.@test qD ≈ qᵣ
-            TT.@test qx ≈ qᵣ
-            TT.@test qD_psd ≈ qᵣ
+            TT.@test qD ≈ qᵣ rtol = 5e-3
+            TT.@test qx ≈ qᵣ rtol = 5e-3
+            TT.@test qD_psd ≈ qᵣ rtol = 5e-3
 
             # Test relationship between exponential moments in diameter space and generalized gamma moments in mass space
             # For raindrops, we expect:
@@ -621,33 +621,33 @@ function test_microphysics2M(FT)
         Mⁿ(n, psd) = y -> y^n * psd(y)
 
         # integral bounds guesstimated in meters for the mass distribution
-        p = eps(FT)
+        p = FT(1e-6)
         D_min, D_max = CM2.get_size_distribution_bounds(pdf_c, ρₐ, qₗ, Nₗ, p)
         x_min = DT.generalized_gamma_quantile(νc, μc, Bc, p)
         x_max = DT.generalized_gamma_quantile(νc, μc, Bc, 1 - p)
 
         # Test that these bounds correspond to the correct probability levels
         TT.@test DT.generalized_gamma_cdf(νc, μc, Bc, x_min) ≈ p
         TT.@test DT.generalized_gamma_cdf(νc, μc, Bc, x_max) ≈ 1 - p
-        TT.@test_broken DT.generalized_gamma_cdf(νcD, μcD, λc, D_min) ≈ p  # TODO: fix this -- not clear why this doesn't work, seems like an issue with SpecialFunctions.jl??
+        TT.@test DT.generalized_gamma_cdf(νcD, μcD, λc, D_min) ≈ p
         TT.@test DT.generalized_gamma_cdf(νcD, μcD, λc, D_max) ≈ 1 - p
 
 
         # Sanity checks of specific content and number concentration with mass distribution
         Nx = QGK.quadgk(Mⁿ(0, f_x), x_min, x_max)[1]
         qx = QGK.quadgk(Mⁿ(1, f_x), x_min, x_max)[1] / ρₐ
-        TT.@test qx ≈ qₗ
-        TT.@test Nx ≈ Nₗ
+        TT.@test qx ≈ qₗ rtol = 1e-5
+        TT.@test Nx ≈ Nₗ rtol = 1e-5
 
         # Sanity checks of specific content and number concentration with diameter distribution
         ND = QGK.quadgk(Mⁿ(0, f_D), D_min, D_max)[1]
         ND_psd = QGK.quadgk(Mⁿ(0, psd), D_min, D_max)[1]
         qD = QGK.quadgk(Mⁿ(3, f_D), D_min, D_max)[1] * k_m / ρₐ
         qD_psd = QGK.quadgk(Mⁿ(3, psd), D_min, D_max)[1] * k_m / ρₐ
-        TT.@test ND ≈ Nₗ
-        TT.@test ND_psd ≈ Nₗ
-        TT.@test qD ≈ qₗ
-        TT.@test qD_psd ≈ qₗ
+        TT.@test ND ≈ Nₗ rtol = 1e-5
+        TT.@test ND_psd ≈ Nₗ rtol = 1e-5
+        TT.@test qD ≈ qₗ rtol = 1e-5
+        TT.@test qD_psd ≈ qₗ rtol = 1e-5
     end
 end