Merge pull request #448 from CliMA/aj/actual_bound_fix

Add lower incomplete gamma function approximation

docs/bibliography.bib

@@ -87,6 +87,17 @@ @article{Bigg1953
   doi = {10.1088/0370-1301/66/8/309}
 }
 
+@article{Blahak2010,
+  title = {Efficient approximation of the incomplete gamma function for use in cloud model applications},
+  author = {Blahak, U.},
+  journal = {Geoscientific Model Development},
+  volume = {3},
+  year = {2010},
+  number = {2},
+  pages = {329--336},
+  doi = {10.5194/gmd-3-329-2010}
+}
+
 @article{BrownFrancis1995,
   title = {Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe},
   author = {Philip R. A. Brown and Peter N. Francis},

docs/src/P3Scheme.md

@@ -134,6 +134,14 @@ As a result ``L_{p3, tot}`` can be expressed as a sum of incomplete gamma functi
 
 where ``\Gamma \,(a, z) = \int_{z}^{\infty} \! t^{a - 1} e^{-t} \mathrm{d}D``
   and ``\Gamma \,(a) = \Gamma \,(a, 0)`` for simplicity.
+We rely on [Blahak2010](@cite) parameterization to compute the incomplete gamma functions.
+Below plot roughly reproduce Figure 4 from [Blahak2010](@cite).
+The relative error differences are mostly due to too small values returned by
+  our reference incomplete gamma function from Julia SpecialFunctions package.
+```@example
+include("plots/P3GammaAprox.jl")
+```
+![](P3_GammaInc_error.svg)
 
 In our solver, we approximate ``\mu`` from ``L/N`` and keep it constant throughout the solving step.
 We approximate ``\mu`` by an exponential function given by the ``L/N`` points
@@ -174,7 +182,7 @@ include("plots/P3LambdaErrorPlots.jl")
 
 ## Terminal Velocity
 
-We use the [Chen2022](@cite) velocity parametrization:
+We use the [Chen2022](@cite) velocity parameterization:
 ```math
 V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D}
 ```
@@ -291,7 +299,7 @@ Visualizing mass-weighted terminal velocity as a function of ``F_{liq}``, ``F_{r
   medium, and large particles, we have mostly graupel (``D_{gr} < D < D_{cr}``) for small ``D_m`` and mostly partially rimed ice
   (``D > D_{cr}``) for medium and large ``D_m``. Thus, we can attribute the non-monotonic behavior of velocity with ``F_liq`` in the
   medium and large ``D_m`` plots to the variations in ``\phi`` caused by nonspherical particle shape, whereas the small ``D_m`` plot
-  confirms a mmore monotonic change in ``V_m`` for spherical ice. The ``L`` and ``N`` values used to generate small, medium, and large ``D_{m}``
+  confirms a more monotonic change in ``V_m`` for spherical ice. The ``L`` and ``N`` values used to generate small, medium, and large ``D_{m}``
   are the same as in the plot above.
 
 ```@example
@@ -301,7 +309,7 @@ include("plots/P3TerminalVelocity_F_liq_rim.jl")
 
 When modifying process rates, we now need to consider whether they are concerned with
   the ice core or the whole particle, in addition to whether they become sources
-  and sinks of different prognostic variables with the inclusion of ``F_{liq}``. 
+  and sinks of different prognostic variables with the inclusion of ``F_{liq}``.
   With the addition of liquid fraction, too,
   come new process rates.
 

docs/src/plots/P3GammaAprox.jl

@@ -0,0 +1,145 @@
+import CairoMakie as PL
+PL.activate!(type = "svg")
+
+import SpecialFunctions as SF
+
+import Thermodynamics as TD
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.P3Scheme as P3
+
+#! format: off
+
+FT = Float64
+res = 100 # plot resolution
+
+# Plot Fig 3 from Blahak 2010
+# ---------------------------
+
+a_range = range(1e-1, stop = 30, length = res)
+c1 = [P3.c₁(a) for a in a_range]
+c2 = [P3.c₂(a) for a in a_range]
+c3 = [P3.c₃(a) for a in a_range]
+c4 = [P3.c₄(a) for a in a_range]
+
+fig = PL.Figure(size = (1500, 1000), fontsize=22, linewidth=3)
+
+ax1 = PL.Axis(fig[1, 1])
+ax2 = PL.Axis(fig[1, 2])
+ax3 = PL.Axis(fig[2, 1])
+ax4 = PL.Axis(fig[2, 2])
+
+ax1.ylabel = "c₁"
+ax2.ylabel = "c₂"
+ax3.ylabel = "c₃"
+ax4.ylabel = "c₄"
+
+ax3.xlabel = "a"
+ax4.xlabel = "a"
+
+PL.ylims!(ax2, 0, 2)
+PL.ylims!(ax4, 0, 10)
+
+c1l = PL.lines!(ax1, a_range,  c1)
+c2l = PL.lines!(ax2, a_range,  c2)
+c3l = PL.lines!(ax3, a_range,  c3)
+c4l = PL.lines!(ax4, a_range,  c4)
+
+PL.save("P3_GammaInc_c_coeffs.svg", fig)
+
+
+# Plot figure 1 from Blahak 2010
+# ------------------------------
+
+fig2 = PL.Figure(size = (1000, 500), fontsize=22)
+ax = PL.Axis(fig2[1,1])
+ax.xlabel = "z"
+ax.ylabel = "P"
+
+z_range = range(1e-1, stop = 25, length=res)
+
+P_a1  = [P3.Γ_lower(FT(1), z)  / SF.gamma(FT(1))  for z in z_range]
+P_a3  = [P3.Γ_lower(FT(3), z)  / SF.gamma(FT(3))  for z in z_range]
+P_a5  = [P3.Γ_lower(FT(5), z)  / SF.gamma(FT(5))  for z in z_range]
+P_a10 = [P3.Γ_lower(FT(10), z) / SF.gamma(FT(10)) for z in z_range]
+P_a15 = [P3.Γ_lower(FT(15), z) / SF.gamma(FT(15)) for z in z_range]
+P_a30 = [P3.Γ_lower(FT(30), z) / SF.gamma(FT(30)) for z in z_range]
+
+G_a1   = [SF.gamma_inc(FT(1), z)[1]  for z in z_range]
+G_a3   = [SF.gamma_inc(FT(3), z)[1]  for z in z_range]
+G_a5   = [SF.gamma_inc(FT(5), z)[1]  for z in z_range]
+G_a10  = [SF.gamma_inc(FT(10), z)[1] for z in z_range]
+G_a15  = [SF.gamma_inc(FT(15), z)[1] for z in z_range]
+G_a30  = [SF.gamma_inc(FT(30), z)[1] for z in z_range]
+
+l1  = PL.lines!(ax, z_range, P_a1,  label="a=1",  color="midnightblue", linewidth=3)
+l3  = PL.lines!(ax, z_range, P_a3,  label="a=3",  color="slateblue",    linewidth=3)
+l5  = PL.lines!(ax, z_range, P_a5,  label="a=5",  color="steelblue1",   linewidth=3)
+l10 = PL.lines!(ax, z_range, P_a10, label="a=10", color="orange",       linewidth=3)
+l15 = PL.lines!(ax, z_range, P_a15, label="a=15", color="orangered2",   linewidth=3)
+l30 = PL.lines!(ax, z_range, P_a30, label="a=30", color="brown",        linewidth=3)
+
+g1  = PL.lines!(ax, z_range, G_a1,  label="a=1", color="midnightblue", linewidth=3, linestyle=:dash)
+g3  = PL.lines!(ax, z_range, G_a3,  label="a=1", color="slateblue",    linewidth=3, linestyle=:dash)
+g5  = PL.lines!(ax, z_range, G_a5,  label="a=1", color="steelblue1",   linewidth=3, linestyle=:dash)
+g10 = PL.lines!(ax, z_range, G_a10, label="a=1", color="orange",       linewidth=3, linestyle=:dash)
+g15 = PL.lines!(ax, z_range, G_a15, label="a=1", color="orangered2",   linewidth=3, linestyle=:dash)
+g30 = PL.lines!(ax, z_range, G_a30, label="a=1", color="brown",        linewidth=3, linestyle=:dash)
+
+PL.Legend(
+  fig2[1, 2],
+  [l1, l3, l5, l10, l15, l30],
+  ["a=1", "a=3", "a=5", "a=10", "a=15", "a=30"],
+)
+PL.save("P3_GammaInc_P_check.svg", fig2)
+
+# Plot figure 4 from Blahak 2010
+# ------------------------------
+# The relative error blows up for cases where the Γ_lower is very small.
+# The reference value we obtain from Julia Special Functions is off
+# by several orders of magnitude when compared with Wolfram Mathematica
+
+y_range = range(1e-1, stop = 2.5, length = res)
+z_range = y_range .* (a_range .+ 1)
+
+val     = zeros(res, res)
+err_rel = zeros(res, res)
+err_abs = zeros(res, res)
+for it in range(1, stop=res, step=1)      # a
+    for jt in range(1, stop=res, step=1)  # z
+        a = FT(a_range[it])
+        z = FT(z_range[jt])
+        val[it, jt] = P3.Γ_lower(a, z) / SF.gamma(a)
+        err_abs[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) - SF.gamma_inc(a, z)[1]
+        err_rel[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) / SF.gamma_inc(a, z)[1] - 1
+    end
+end
+
+fig3 = PL.Figure(size = (1800, 600), fontsize=22)
+ax1 = PL.Axis(fig3[1, 1], xlabel = "a", ylabel = "z  / (a + 1)", title = "Relative error")
+ax2 = PL.Axis(fig3[1, 3], xlabel = "a", title = "Absolute error")
+ax3 = PL.Axis(fig3[1, 5], xlabel = "a", title = "Γ_lower(a, z) / Γ(a)")
+
+PL.ylims!(ax1, 0, 2.5)
+PL.ylims!(ax2, 0, 2.5)
+PL.ylims!(ax3, 0, 2.5)
+
+cplot1 = PL.heatmap!(
+    ax1, a_range, y_range, err_rel,
+    colormap = PL.cgrad(:viridis, 20, categorical=true),
+    colorrange = (-0.1, 10), highclip=:red, lowclip=:orange,
+)
+PL.Colorbar(fig3[1,2], cplot1)
+
+cplot2 = PL.contourf!(
+    ax2, a_range, y_range, err_abs,
+    levels=-0.03:0.005:0.03,
+    colormap="RdBu"
+)
+PL.Colorbar(fig3[1,4], cplot2)
+
+cplot3 = PL.contourf!(
+    ax3, a_range, y_range, val,
+)
+PL.Colorbar(fig3[1,6], cplot3)
+
+PL.save("P3_GammaInc_error.svg", fig3)

src/P3_helpers.jl

@@ -1,8 +1,99 @@
-# Some wrappers to cast types from SF.gamma
+# Wrapper to cast types from SF.gamma
 # (which returns Float64 even when the input is Float32)
-Γ(a::FT, z::FT) where {FT <: Real} = FT(SF.gamma(a, z))
+# TODO - replace with parameterization of our own
 Γ(a::FT) where {FT <: Real} = FT(SF.gamma(a))
 
+# helper functions for Γ_approx
+function c₁(a::FT) where {FT <: Real}
+    p1 = FT(9.4368392235e-03)
+    p2 = -FT(1.0782666481e-04)
+    p3 = -FT(5.8969657295e-06)
+    p4 = FT(2.8939523781e-07)
+    p5 = FT(1.0043326298e-01)
+    p6 = FT(5.5637848465e-01)
+    p = [p1, p2, p3, p4, p5, p6]
+
+    ret = 1 + p[5] * (exp(-p[6] * a) - 1)
+    for it in range(1, 4, step = 1)
+        ret += p[it] * a^it
+    end
+    return ret
+end
+function c₂(a::FT) where {FT <: Real}
+    q1 = FT(1.1464706419e-01)
+    q2 = FT(2.6963429121)
+    q3 = -FT(2.9647038257)
+    q4 = FT(2.1080724954)
+    q = [q1, q2, q3, q4]
+
+    ret = q[1]
+    for it in range(2, 4, step = 1)
+        ret += q[it] / a^(it - 1)
+    end
+    return ret
+end
+function c₃(a::FT) where {FT <: Real}
+    r1 = FT(0)
+    r2 = FT(1.1428716184)
+    r3 = -FT(6.6981186438e-03)
+    r4 = FT(1.0480765092e-04)
+    r = [r1, r2, r3, r4]
+
+    ret = 0
+    for it in range(1, 4, step = 1)
+        ret += r[it] * a^(it - 1)
+    end
+    return ret
+end
+function c₄(a::FT) where {FT <: Real}
+    s1 = FT(1.0356711153)
+    s2 = FT(2.3423452308)
+    s3 = -FT(3.6174503174e-01)
+    s4 = -FT(3.1376557650)
+    s5 = FT(2.9092306039)
+    s = [s1, s2, s3, s4, s5]
+
+    ret = 0
+    for it in range(1, 5, step = 1)
+        ret += s[it] / a^(it - 1)
+    end
+    return ret
+end
+
+"""
+    Γ_lower(a, z)
+
+An approximated lower incomplete Gamma function based on Blahak 2010:
+doi:10.5194/gmd-3-329-2010
+https://gmd.copernicus.org/articles/3/329/2010/gmd-3-329-2010.pdf
+"""
+function Γ_lower(a::FT, z::FT) where {FT <: Real}
+
+    if isnan(z) || z <= FT(0)
+        return FT(0)
+    else
+        W(a, z) = FT(0.5) + FT(0.5) * tanh(c₂(a) * (z - c₃(a)))
+
+        return exp(-z) *
+               z^a *
+               (
+                   1 / a +
+                   c₁(a) * z / a / (a + 1) +
+                   (c₁(a) * z)^2 / a / (a + 1) / (a + 2)
+               ) *
+               (1 - W(a, z)) + Γ(a) * W(a, z) * (1 - c₄(a)^(-z))
+    end
+end
+
+"""
+    Γ_upper(a, z)
+
+An approximated upper incomplete Gamma function computed as Γ(a) - Γ_lower(a, z)
+"""
+function Γ_upper(a::FT, z::FT) where {FT <: Real}
+    return Γ(a) - Γ_lower(a, z)
+end
+
 """
     ∫_Γ(x₀, x_end, c1, c2, c3)
 
@@ -16,11 +107,13 @@ Integrates f(D, c1, c2, c3) dD from x₀ to x_end
 """
 function ∫_Γ(x₀::FT, x_end::FT, c1::FT, c2::FT, c3::FT) where {FT}
     if x_end == FT(Inf)
-        return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3))
+        return c1 * c3^(-c2 - 1) * (Γ_upper(1 + c2, x₀ * c3))
     elseif x₀ == 0
-        return c1 * c3^(-c2 - 1) * (Γ(1 + c2) - Γ(1 + c2, x_end * c3))
+        return c1 * c3^(-c2 - 1) * (Γ_lower(1 + c2, x_end * c3))
     else
-        return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3) - Γ(1 + c2, x_end * c3))
+        return c1 *
+               c3^(-c2 - 1) *
+               (Γ_upper(1 + c2, x₀ * c3) - Γ_upper(1 + c2, x_end * c3))
     end
 end
 

src/P3_particle_properties.jl

@@ -312,12 +312,12 @@ end
  - F_rim - rime mass fraction (L_rim/ L_ice) [-]
  - th - P3 particle properties thresholds
 
-Returns the aspect ratio (ϕ) for an ice particle with mass, cross-sectional area,  
-and ice density determined using the size-dependent particle property regimes  
-following Morrison and Milbrandt (2015). The density of nonspherical  
-particles is assumed to be equal to the particle mass divided by the volume of a  
-spherical particle with the same D_max.  
-Assuming zero liquid fraction. 
+Returns the aspect ratio (ϕ) for an ice particle with mass, cross-sectional area,
+and ice density determined using the size-dependent particle property regimes
+following Morrison and Milbrandt (2015). The density of nonspherical
+particles is assumed to be equal to the particle mass divided by the volume of a
+spherical particle with the same D_max.
+Assuming zero liquid fraction.
 """
 function ϕᵢ(p3::PSP3, D::FT, F_rim::FT, th) where {FT}
     F_liq = FT(0)

src/P3_processes.jl

@@ -90,7 +90,7 @@ function ice_melt(
         distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_)
     N(D) = N′ice(p3, D, λ, N_0)
     # ... and D_max for the integral
-    bound = get_ice_bound(p3, λ, FT(1e-6))
+    bound = get_ice_bound(p3, λ, N_ice, FT(1e-6))
 
     # Ice particle terminal velocity
     v(D) = ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th)