Merge pull request #293 from CliMA/apaj/p3_docs

Move P3 mass functions to src, add gamma functions to docs

docs/src/API.md

@@ -58,6 +58,7 @@ Microphysics2M.conv_q_liq_to_q_rai
 ```@docs
 P3Scheme
 P3Scheme.thresholds
+P3Scheme.p3_mass
 ```
 
 # Aerosol model

docs/src/P3Scheme.md

@@ -20,43 +20,6 @@ The prognostic variables are:
 !!! note
     TODO - At some point we should switch to specific humidities...
 
-## Assumed particle size distribution (PSD)
-
-Following [MorrisonMilbrandt2015](@cite), the scheme assumes a
- gamma distribution for the concentration of particles per unit volume
- based on particle size measurements obtained by [Heymsfield2003](@cite)
- in tropical and midlatitude ice clouds and implemented by
- [MorrisonGrabowski2008](@cite):
-
-```math
-N'(D) = N_{0} D^\mu \, e^{-\lambda \, D}
-```
-where:
- - ``N'`` is the number concentration in ``m^{-4}``
- - ``D`` is the maximum particle dimension in ``m``,
- - ``N_0`` is the intercept parameter in ``m^{-4}``,
- - ``\mu`` is the shape parameter (dimensionless),
- - ``\lambda`` is the slope parameter in ``m^{-1}``.
-
-We assume ``\mu \ = 0.00191 \; \lambda \ ^{0.8} - 2``.
-Following [MorrisonGrabowski2008](@cite) we limit ``\mu \ \in (0,6)``.
-A negative ``\mu`` can occur only for very small mean particle sizes``\frac{1}{\lambda} < ~0.17 mm``.
-``N_0`` and ``\lambda`` can be found using different moments of the PSD,
-namely the total number concentration ``N`` and mass mixing ratio ``q``, where
-
-```math
-N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D
-```
-
-```math
-q = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D
-```
-
-For liquid droplets, these equations are solved without issue, but for ice, the third moment of the size distribution listed above (i.e. ``\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D``) varies as the mass relation varies across the PSD (see below for the mass regime documentation). On the other hand, the first moment of the PSD, the number concentration, does not vary across the PSD and yields ``N = \frac{N_0}{\lambda}``.
-
-!!! note
-    TODO - The scheme uses a mean particle size value ``D_m`` for each time step to determine which mass relation to employ. In other words, ``N_0`` and ``\lambda`` must be calculated for the five different mass relations below to accommodate ranges of ``D_m`` corresponding to each mass relation. For mean particle sizes that employ the mass relations characterized by graupel and by partially rimed ice, the mass relations are time-dependent due to the presence of ``\rho_g`` and ``F_r``. This complicates the scheme's use of the PSD, and as a result, writing analytical formulas for the PSD parameters is challenging.
-
 ## Assumed particle mass relationships
 
 The mass ``m`` of particles as a function of maximum particle dimension ``D``
@@ -91,10 +54,6 @@ The remaining thresholds: ``D_{gr}``, ``D_{cr}``, as well as the
 where
  - ``F_r = \frac{q_{rim}}{q_{ice}}`` is the rime mass fraction,
  - ``\rho_{r} = \frac{q_{rim}}{B_{rim}}`` is the predicted rime density.
-The system is solved using [`NonlinearSolve.jl`](https://docs.sciml.ai/NonlinearSolve/stable/).
-
-!!! note
-    TODO - The use of NonlinearSolve.jl is not ideal because of its runtime and memory allocation requirements. Currently, there is also a look-up table NetCDF file which could be used to look up values of the quantities which form the nonlinear system. However, the look-up table is not GPU-compatible and would require too much memory in an Earth System Model. The current approach may be of use for testing and for visualization of the system, but other options, such as using RootSolvers.jl or using a simpler fit that approximates the solver output, are more suitable long term solutions which do not require outside packages which employ auto-differentiation or use memory, both of which do not suit the needs of CliMA.
 
 ## Assumed particle projected area relationships
 
@@ -116,7 +75,63 @@ where all variables from the m(D) regime are as defined above, and:
  - ``\sigma = 1.88`` (dimensionless), both from the aggregates of side planes, columns, bullets, and planar polycrystals in [Mitchell1996](@cite).
 
 !!! note
-    TODO - As mentioned in [issue #151](https://github.com/CliMA/CloudMicrophysics.jl/issues/151), the units of ``\gamma`` and ``\sigma`` are not immediately clear from [Mitchell1996](@cite) and [MorrisonMilbrandt2015](@cite). To resolve this issue, it may be useful to contact the authors of the paper, or, examine the representative figures below to check units. It has occured to me that the units of ``D`` are probably m and that the units of ``A`` are probably m2. I have assumed these dimensions for the time being. Another likely scenario would be if ``D`` had units of mm, in which case we would have ``\gamma = 0.2285 \; 10^{3 \sigma}`` to correct for units. However, the plots of area versus particle dimension look outlandish in this case.
+    TODO - Check units, see in [issue #151](https://github.com/CliMA/CloudMicrophysics.jl/issues/151)
+
+## Assumed particle size distribution (PSD)
+
+Following [MorrisonMilbrandt2015](@cite), the scheme assumes a
+ gamma distribution for the concentration of ice particles per unit volume
+ based on particle size measurements obtained by [Heymsfield2003](@cite)
+ in tropical and midlatitude ice clouds and implemented by
+ [MorrisonGrabowski2008](@cite):
+
+```math
+N'(D) = N_{0} D^\mu \, e^{-\lambda \, D}
+```
+where:
+ - ``N'`` is the number concentration in ``m^{-4}``
+ - ``D`` is the maximum particle dimension in ``m``,
+ - ``N_0`` is the intercept parameter in ``m^{-4}``,
+ - ``\mu`` is the shape parameter (dimensionless),
+ - ``\lambda`` is the slope parameter in ``m^{-1}``.
+
+We assume ``\mu \ = 0.00191 \; \lambda \ ^{0.8} - 2``.
+Following [MorrisonGrabowski2008](@cite) we limit ``\mu \ \in (0,6)``.
+A negative ``\mu`` can occur only for very small mean particle sizes``\frac{1}{\lambda} < ~0.17 mm``.
+
+The model predicted ice number concentration and ice mass density are defined as
+```math
+N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D
+```
+```math
+q_{ice} = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D
+```
+
+## Calculating shape parameters
+
+As a next step we need to find the mapping between
+  predicted moments of the size distribution ``N_{ice}`` and ``q_{ice}``
+  and the shape parameters ``\lambda`` and ``N_0``.
+Solving for ``N_{ice}`` is relatively straightforward:
+```math
+N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D = \int_{0}^{\infty} \! N_{0} D^\mu \, e^{-\lambda \, D} \mathrm{d}D = N_{0} (\lambda \,^{-(\mu \, + 1)} \Gamma \,(1 + \mu \,))
+```
+
+``q_{ice}`` depends on the variable mass-size relation ``m(D)`` defined above.
+We solve for ``q_{ice}`` in a piece-wise fashion defined by the same thresholds as ``m(D)``.
+As a result ``q\_{ice}`` can be expressed as a sum of inclomplete gamma functions.
+  and the shape parameters are found using iterative solver.
+
+|      condition(s)                            |    ``q_{ice} = \int \! m(D) N'(D) \mathrm{d}D``                                          |         gamma representation          |
+|:---------------------------------------------|:-----------------------------------------------------------------------------------------|:---------------------------------------------|
+| ``D < D_{th}``                               | ``\int_{0}^{D_{th}} \! \frac{\pi}{6} \rho_i \ D^3 N'(D) \mathrm{d}D``                    | ``\frac{\pi}{6} \rho_i N_0 \lambda \,^{-(\mu \, + 4)} (\Gamma \,(\mu \, + 4) - \Gamma \,(\mu \, + 4, \lambda \,D_{th}))``|
+| ``q_{rim} = 0`` and ``D > D_{th}``           | ``\int_{D_{th}}^{\infty} \! \alpha_{va} \ D^{\beta_{va}} N'(D) \mathrm{d}D``             | ``\alpha_{va} \ N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1) + \Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{th}) - (\mu \, + \beta_{va} \,)\Gamma \,(\mu \, + \beta_{va} \,))`` |
+| ``q_{rim} > 0`` and ``D_{gr} > D > D_{th}``  | ``\int_{D_{th}}^{D_{gr}} \! \alpha_{va} \ D^{\beta_{va}} N'(D) \mathrm{d}D``             | ``\alpha_{va} \ N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{th}) - \Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{gr}))`` |
+| ``q_{rim} > 0`` and ``D_{cr} > D > D_{gr}``  | ``\int_{D_{gr}}^{D_{cr}} \! \frac{\pi}{6} \rho_g \ D^3 N'(D) \mathrm{d}D``               | ``\frac{\pi}{6} \rho_g N_0 \lambda \,^{-(\mu \, + 4)} (\Gamma \,(\mu \, + 4, \lambda \,D_{gr}) - \Gamma \,(\mu \, + 4, \lambda \,D_{cr}))`` |
+| ``q_{rim} > 0`` and ``D > D_{cr}``           | ``\int_{D_{cr}}^{\infty} \! \frac{\alpha_{va}}{1-F_r} D^{\beta_{va}} N'(D) \mathrm{d}D`` | ``\frac{\alpha_{va}}{1-F_r} N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1) + \Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{cr}) - (\mu \, + \beta_{va} \,)\Gamma \,(\mu \, + \beta_{va} \,))``  |
+
+where ``\Gamma \,(a, z) = \int_{z}^{\infty} \! t^{a - 1} e^{-t} \mathrm{d}D``
+  and ``\Gamma \,(a) = \Gamma \,(a, 0)`` for simplicity.
 
 ## Example figures
 

docs/src/plots/P3SchemePlots.jl

@@ -9,22 +9,7 @@ FT = Float64
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
 
-"""
-    mass_(p3, D, ρ, F_r)
-
- - p3 - a struct with P3 scheme parameters
- - D - maximum particle dimension [m]
- - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3]
- - F_r - rime mass fraction [q_rim/q_i]
 
-Returns mass as a function of size for differen particle regimes [kg]
-"""
-# for spherical ice (small ice or completely rimed ice)
-mass_s(D::FT, ρ::FT) = FT(π) / 6 * ρ * D^3
-# for large nonspherical ice (used for unrimed and dense types)
-mass_nl(p3::PSP3, D::FT) = P3.α_va_si(p3) * D^p3.β_va
-# for partially rimed ice
-mass_r(p3::PSP3, D::FT, F_r::FT) = P3.α_va_si(p3) / (1 - F_r) * D^p3.β_va
 
 """
     A_(p3, D)
@@ -40,38 +25,6 @@ A_s(D::FT) = FT(π) / 4 * D^2
 A_ns(p3::PSP3, D::FT) = p3.γ * D^p3.σ
 # partially rimed ice
 A_r(p3::PSP3, F_r::FT, D::FT) = F_r * A_s(D) + (1 - F_r) * A_ns(p3, D)
-"""
-    mass(p3, D, F_r, th)
-
- - p3 - a struct with P3 scheme parameters
- - D - maximum particle dimension
- - F_r - rime mass fraction (q_rim/q_i)
- - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
-
-Returns mass(D) regime, used to create figures for the docs page.
-"""
-function mass(
-    p3::PSP3,
-    D::FT,
-    F_r::FT,
-    th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
-) where {FT <: Real}
-    if P3.D_th_helper(p3) > D
-        return mass_s(D, p3.ρ_i)          # small spherical ice
-    end
-    if F_r == 0
-        return mass_nl(p3, D)             # large nonspherical unrimed ice
-    end
-    if th.D_gr > D >= P3.D_th_helper(p3)
-        return mass_nl(p3, D)             # dense nonspherical ice
-    end
-    if th.D_cr > D >= th.D_gr
-        return mass_s(D, th.ρ_g)          # graupel
-    end
-    if D >= th.D_cr
-        return mass_r(p3, D, F_r)         # partially rimed ice
-    end
-end
 
 """
     area(p3, D, F_r, th)
@@ -137,9 +90,9 @@ function p3_mass_plot()
     sol_8 = P3.thresholds(p3, 400.0, 0.8)
 
     #! format: off
-    fig1_a_0 = Plt.lines!(ax1_a, D_range * 1e3, [mass(p3, D, 0.0       ) for D in D_range], color = cl[1], linewidth = lw)
-    fig1_a_5 = Plt.lines!(ax1_a, D_range * 1e3, [mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw)
-    fig1_a_8 = Plt.lines!(ax1_a, D_range * 1e3, [mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig1_a_0 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.0       ) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_a_5 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_a_8 = Plt.lines!(ax1_a, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw)
 
     d_tha  = Plt.vlines!(ax1_a, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
     d_cr_5 = Plt.vlines!(ax1_a, sol_5[1]           * 1e3, linestyle = :dot,  color = cl[2], linewidth = lw)
@@ -177,9 +130,9 @@ function p3_mass_plot()
     sol_8 = P3.thresholds(p3, 800.0, 0.95)
 
     #! format: off
-    fig1_b200 = Plt.lines!(ax1_b, D_range * 1e3, [mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
-    fig1_b400 = Plt.lines!(ax1_b, D_range * 1e3, [mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
-    fig1_b800 = Plt.lines!(ax1_b, D_range * 1e3, [mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig1_b200 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_b400 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_b800 = Plt.lines!(ax1_b, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
 
     d_thb    = Plt.vlines!(ax1_b, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
     d_cr_200 = Plt.vlines!(ax1_b, sol_2[1] * 1e3,           linestyle = :dot,  color = cl[1], linewidth = lw)

src/P3Scheme.jl

@@ -11,12 +11,12 @@ Note: Particle size is defined as its maximum length (i.e. max dimesion).
 module P3Scheme
 
 import RootSolvers as RS
-import CLIMAParameters as CP
-import ..Parameters as CMP
+import SpecialFunctions as SF
 
+import CloudMicrophysics.Parameters as CMP
 const PSP3 = CMP.ParametersP3
 
-export thresholds
+export thresholds, p3_mass
 
 """
     α_va_si(p3)
@@ -148,4 +148,55 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT}
     )
 end
 
+"""
+    mass_(p3, D, ρ, F_r)
+
+ - p3 - a struct with P3 scheme parameters
+ - D - maximum particle dimension [m]
+ - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3]
+ - F_r - rime mass fraction [q_rim/q_i]
+
+Returns mass as a function of size for differen particle regimes [kg]
+"""
+# for spherical ice (small ice or completely rimed ice)
+mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3
+# for large nonspherical ice (used for unrimed and dense types)
+mass_nl(p3::PSP3, D::FT) where {FT <: Real} = α_va_si(p3) * D^p3.β_va
+# for partially rimed ice
+mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} =
+    α_va_si(p3) / (1 - F_r) * D^p3.β_va
+
+"""
+    p3_mass(p3, D, F_r, th)
+
+ - p3 - a struct with P3 scheme parameters
+ - D - maximum particle dimension
+ - F_r - rime mass fraction (q_rim/q_i)
+ - th - P3Scheme nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
+
+Returns mass(D) regime, used to create figures for the docs page.
+"""
+function p3_mass(
+    p3::PSP3,
+    D::FT,
+    F_r::FT,
+    th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
+) where {FT <: Real}
+    if D_th_helper(p3) > D
+        return mass_s(D, p3.ρ_i)          # small spherical ice
+    end
+    if F_r == 0
+        return mass_nl(p3, D)             # large nonspherical unrimed ice
+    end
+    if th.D_gr > D >= D_th_helper(p3)
+        return mass_nl(p3, D)             # dense nonspherical ice
+    end
+    if th.D_cr > D >= th.D_gr
+        return mass_s(D, th.ρ_g)          # graupel
+    end
+    if D >= th.D_cr
+        return mass_r(p3, D, F_r)         # partially rimed ice
+    end
+end
+
 end

test/p3_tests.jl

@@ -84,8 +84,74 @@ function test_p3_thresholds(FT)
     end
 end
 
+function test_p3_mass(FT)
+
+    p3 = CMP.ParametersP3(FT)
+
+    TT.@testset "p3 mass_ functions" begin
+
+        # Initialize test values
+        Ds = (FT(1e-5), FT(1e-4), FT(1e-3))
+        ρs = (FT(400), FT(600), FT(800))
+        F_rs = (FT(0.0), FT(0.5), FT(0.8))
+
+        # Test to see that the mass_ functions give positive, not NaN values
+        for D in Ds
+            for ρ in ρs
+                for F_r in F_rs
+                    TT.@test P3.mass_s(D, ρ) >= 0
+                    TT.@test P3.mass_nl(p3, D) >= 0
+                    TT.@test P3.mass_r(p3, D, F_r) >= 0
+                end
+            end
+        end
+    end
+
+    # Test to see that p3_mass gives correct mass 
+    TT.@testset "p3_mass() accurate values" begin
+
+        # Initialize test values
+        Ds = (FT(1e-5), FT(1e-4), FT(1e-3))
+        ρs = (FT(400), FT(600), FT(800))
+        F_rs = (FT(0.0), FT(0.5), FT(0.8))
+        eps = FT(1e-3)
+
+        for ρ in ρs
+            for F_r in F_rs
+                D_th = P3.D_th_helper(p3)
+                D1 = D_th / 2
+
+                if (F_r > 0)
+                    th = P3.thresholds(p3, ρ, F_r)
+
+                    D2 = (th.D_gr + D_th) / 2
+                    D3 = (th.D_cr + th.D_gr) / 2
+                    D4 = th.D_cr + eps
+
+                    TT.@test P3.p3_mass(p3, D1, F_r, th) ==
+                             P3.mass_s(D1, p3.ρ_i)
+                    TT.@test P3.p3_mass(p3, D2, F_r, th) == P3.mass_nl(p3, D2)
+                    TT.@test P3.p3_mass(p3, D3, F_r, th) ==
+                             P3.mass_s(D3, th.ρ_g)
+                    TT.@test P3.p3_mass(p3, D4, F_r, th) ==
+                             P3.mass_r(p3, D4, F_r)
+                else
+                    D2 = D1 + eps
+                    TT.@test P3.p3_mass(p3, D1, F_r, ()) ==
+                             P3.mass_s(D1, p3.ρ_i)
+                    TT.@test P3.p3_mass(p3, D2, F_r, ()) == P3.mass_nl(p3, D2)
+                end
+
+            end
+        end
+    end
+
+end
+
 println("Testing Float32")
 test_p3_thresholds(Float32)
+test_p3_mass(Float32)
 
 println("Testing Float64")
 test_p3_thresholds(Float64)
+test_p3_mass(Float64)