Merge pull request #397 from CliMA/ap/integral

P3 integral vs gamma tests
CliMA · May 29, 2024 · ffcacc4 · ffcacc4
2 parents bf68040 + 6bee931
commit ffcacc4
Show file tree

Hide file tree

Showing 11 changed files with 413 additions and 105 deletions.
diff --git a/Project.toml b/Project.toml
@@ -7,6 +7,7 @@ version = "0.20.2"
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -24,6 +25,7 @@ ClimaParams = "0.10.6"
 DataFrames = "1.6"
 DocStringExtensions = "0.8, 0.9"
 ForwardDiff = "0.10"
+HCubature = "1.6"
 MLJ = "0.20"
 QuadGK = "2.9"
 RootSolvers = "0.3, 0.4"

diff --git a/docs/src/API.md b/docs/src/API.md
@@ -65,6 +65,12 @@ P3Scheme.thresholds
 P3Scheme.distribution_parameter_solver
 ```
 
+# Terminal Velocity
+```@docs 
+TerminalVelocity 
+TerminalVelocity.velocity_chen
+```
+
 # Aerosol model
 
 ```@docs

diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl
@@ -9,99 +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) 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} = P3.α_va_si(p3) * D^p3.β_va
-# for partially rimed ice
-mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} =
-    P3.α_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}
-    D_th = P3.D_th_helper(p3)
-    if D_th > D
-        return mass_s(D, p3.ρ_i)          # small spherical ice
-    elseif F_r == 0
-        return mass_nl(p3, D)             # large nonspherical unrimed ice
-    elseif th.D_gr > D >= D_th
-        return mass_nl(p3, D)             # dense nonspherical ice
-    elseif th.D_cr > D >= th.D_gr
-        return mass_s(D, th.ρ_g)          # graupel
-    elseif D >= th.D_cr
-        return mass_r(p3, D, F_r)         # partially rimed ice
-    end
-end
-
-"""
-    A_(p3, D)
-
- - p3 - a struct with P3 scheme parameters
- - D - maximum particle dimension
-
-Returns particle projected area as a function of size for different particle regimes
-"""
-# for spherical particles
-A_s(D::FT) = FT(π) / 4 * D^2
-# for nonspherical particles
-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)
-
-"""
-    area(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 area(D), used to create figures for the documentation.
-"""
-function area(
-    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}
-    # Area regime:
-    if P3.D_th_helper(p3) > D
-        return A_s(D)                      # small spherical ice
-    elseif F_r == 0
-        return A_ns(p3, D)                 # large nonspherical unrimed ice
-    elseif th.D_gr > D >= P3.D_th_helper(p3)
-        return A_ns(p3, D)                 # dense nonspherical ice
-    elseif th.D_cr > D >= th.D_gr
-        return A_s(D)                      # graupel
-    elseif D >= th.D_cr
-        return A_r(p3, F_r, D)             # partially rimed ice
-    end
-
-end
 
 function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect)
     return CMK.Axis(
@@ -143,13 +51,13 @@ function p3_relations_plot()
     sol4_5 = P3.thresholds(p3, 400.0, 0.5)
     sol4_8 = P3.thresholds(p3, 400.0, 0.8)
     # m(D)
-    fig1_0 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.0        ) for D in D_range], color = cl[1], linewidth = lw)
-    fig1_5 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
-    fig1_8 = CMK.lines!(ax1, D_range * 1e3, [p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.0        ) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
+    fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
     # a(D)
-    fig2_0 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.0        ) for D in D_range], color = cl[1], linewidth = lw)
-    fig2_5 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
-    fig2_8 = CMK.lines!(ax2, D_range * 1e3, [area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0        ) for D in D_range], color = cl[1], linewidth = lw)
+    fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
+    fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
     # plot verical lines
     for ax in [ax1, ax2]
         global d_tha  = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
@@ -164,13 +72,13 @@ function p3_relations_plot()
     sol_4 = P3.thresholds(p3, 400.0, 0.95)
     sol_8 = P3.thresholds(p3, 800.0, 0.95)
     # m(D)
-    fig3_200 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
-    fig3_400 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
-    fig3_800 = CMK.lines!(ax3, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig3_200 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig3_400 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig3_800 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw)
     # a(D)
-    fig3_200 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw)
-    fig3_400 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
-    fig3_800 = CMK.lines!(ax4, D_range * 1e3, [area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
     # plot verical lines
     for ax in [ax3, ax4]
         global d_thb    = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)

diff --git a/src/CloudMicrophysics.jl b/src/CloudMicrophysics.jl
@@ -14,6 +14,7 @@ include("Common.jl")
 include("Microphysics0M.jl")
 include("Microphysics1M.jl")
 include("Microphysics2M.jl")
+include("TerminalVelocity.jl")
 include("P3Scheme.jl")
 include("MicrophysicsNonEq.jl")
 include("AerosolModel.jl")

diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl
@@ -14,10 +14,12 @@ module P3Scheme
 import SpecialFunctions as SF
 import QuadGK as QGK
 import RootSolvers as RS
+import HCubature as HC
 
 import ClimaParams as CP
 import CloudMicrophysics.Parameters as CMP
 import CloudMicrophysics.Common as CO
+import CloudMicrophysics.TerminalVelocity as TV
 
 const PSP3 = CMP.ParametersP3
 
@@ -156,6 +158,102 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT}
     end
 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}
+    D_th = D_th_helper(p3)
+    if D_th > D
+        return mass_s(D, p3.ρ_i)          # small spherical ice
+    elseif F_r == 0
+        return mass_nl(p3, D)             # large nonspherical unrimed ice
+    elseif th.D_gr > D >= D_th
+        return mass_nl(p3, D)             # dense nonspherical ice
+    elseif th.D_cr > D >= th.D_gr
+        return mass_s(D, th.ρ_g)          # graupel
+    elseif D >= th.D_cr
+        return mass_r(p3, D, F_r)         # partially rimed ice
+    end
+end
+
+"""
+    A_(p3, D)
+
+ - p3 - a struct with P3 scheme parameters
+ - D - maximum particle dimension
+
+Returns particle projected area as a function of size for different particle regimes
+"""
+# for spherical particles
+A_s(D::FT) where {FT <: Real} = FT(π) / 4 * D^2
+# for nonspherical particles
+A_ns(p3::PSP3, D::FT) where {FT <: Real} = p3.γ * D^p3.σ
+# partially rimed ice
+A_r(p3::PSP3, F_r::FT, D::FT) where {FT <: Real} =
+    F_r * A_s(D) + (1 - F_r) * A_ns(p3, D)
+
+"""
+    p3_area(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 area(D), used to create figures for the documentation.
+"""
+function p3_area(
+    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}
+    # Area regime:
+    if D_th_helper(p3) > D
+        return A_s(D)                      # small spherical ice
+    elseif F_r == 0
+        return A_ns(p3, D)                 # large nonspherical unrimed ice
+    elseif th.D_gr > D >= D_th_helper(p3)
+        return A_ns(p3, D)                 # dense nonspherical ice
+    elseif th.D_cr > D >= th.D_gr
+        return A_s(D)                      # graupel
+    elseif D >= th.D_cr
+        return A_r(p3, F_r, D)             # partially rimed ice
+    else
+        throw("D not in range")
+    end
+end
+
 # Some wrappers 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))
@@ -475,7 +573,7 @@ function terminal_velocity(
     # Get the thresholds for different particles regimes
     (; D_cr, D_gr, ρ_g, ρ_d) = thresholds(p3, ρ_r, F_r)
     D_th = D_th_helper(p3)
-    D_ct = FT(0.000625) # TODO add to ClimaParams
+    D_ct = Chen2022.cutoff
 
     # Get the shape parameters of the particle size distribution
     (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r)
@@ -667,4 +765,49 @@ function D_m(p3::PSP3, q::FT, N::FT, ρ_r::FT, F_r::FT) where {FT}
     return n / q
 end
 
+"""
+    N′ice(p3, D, λ, N0)
+
+ - p3 - a struct containing P3 scheme parameters
+ - D - diameter of particle 
+ - λ - shape parameter of distribution 
+ - N0 - shape parameter of distribution
+
+ Returns the distribution of ice particles (assumed to be of the form 
+ N'(D) = N0 * D ^ μ * exp(-λD)) at given D 
+"""
+function N′ice(p3::PSP3, D::FT, λ::FT, N_0::FT) where {FT}
+    return N_0 * D^DSD_μ(p3, λ) * exp(-λ * D)
+end
+
+"""
+    get_ice_bounds(p3, λ, tolerance)
+
+ - p3 - a struct containing p3 parameters 
+ - λ - shape parameters of ice distribution 
+ - tolerance - tolerance for how much of the distribution we want to integrate over
+
+ Returns the bound on the distribution that would guarantee that 1-tolerance 
+ of the ice distribution is integrated over. This is calculated by setting 
+ N_0(1 - tolerance) = ∫ N'(D) dD from 0 to bound and solving for bound. 
+ This was further simplified to cancel out the N_0 from both sides. 
+ The guess was calculated through a linear approximation extrapolated from 
+ numerical solutions. 
+"""
+function get_ice_bound(p3, λ::FT, tolerance::FT) where {FT}
+    ice_problem(x) =
+        tolerance - Γ(1 + DSD_μ(p3, λ), FT(exp(x)) * λ) / Γ(1 + DSD_μ(p3, λ))
+    guess = log(19 / 6 * (DSD_μ(p3, λ) - 1) + 39) - log(λ)
+    log_ice_x =
+        RS.find_zero(
+            ice_problem,
+            RS.SecantMethod(guess - 1, guess),
+            RS.CompactSolution(),
+            RS.RelativeSolutionTolerance(eps(FT)),
+            5,
+        ).root
+
+    return exp(log_ice_x)
+end
+
 end