Lambda testing

Project.toml

@@ -7,6 +7,7 @@ version = "0.15.2"
 CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Integrals = "de52edbc-65ea-441a-8357-d3a637375a31"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"

docs/src/P3Scheme.md

@@ -144,3 +144,29 @@ include("plots/P3SchemePlots.jl")
 ![](MorrisonandMilbrandtFig1b.svg)
 ![](P3Scheme_Area_1.svg)
 ![](P3Scheme_Area_2.svg)
+
+
+## Calculating shape parameters 
+
+Solving  for the shape parameters of the ``N'`` distribution proves to be difficult for ice particles. We solve the non-linear system from ``N`` and ``q`` for ``N_0`` and ``\lambda`` through the use of incomplete gamma functions. 
+
+Solving the function for ``N`` is relatively straightforward.
+
+```math 
+N = \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 \,))
+```
+
+Calculating ``q`` in terms of incomplete gamma functions proves to be increasingly difficult due to the variation in the mass relation across the PSD (as stated above). We solve ``q`` in a piece-wise fashion defined by the same thresholds as ``m(D)``:
+
+|      condition(s)                            |    ``q = \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_{cr}}^{D_{gr}} \! \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.
+
+Therefore, for any ``q_{rim}`` we are able to piece together a total ``q`` value across the PSD using these incomplete gamma functions. 
+

src/P3Scheme.jl

@@ -228,7 +228,7 @@ Returns the shape parameter μ for a given λ value
 Eq. 3 in Morrison and Milbrandt (2015).
 """
 function DSD_μ(p3::PSP3, λ::FT) where {FT}
-    @assert λ > FT(0)
+    #@assert λ > FT(0)
     return min(p3.μ_max, max(FT(0), p3.a * λ^p3.b - p3.c))
 end
 
@@ -317,6 +317,46 @@ function q_gamma(
     )
 end
 
+"""
+    get_bounds(N, q, F_r, p3, th)
+
+ - N - ice number concentration [1/m3]
+ - q - mass mixing ratio
+ - F_r -rime mass fraction [q_rim/q_i]
+ - p3 - a struct with P3 scheme parameters
+ - th -  thresholds() nonlinear solve output tuple (D_cr, D_gr, ρ_g, ρ_d)
+
+ Returns estimated guess for λ from q to be used in distribution_parameter_solver()
+"""
+function get_bounds(
+    N::FT,
+    N̂::FT,
+    q::FT,
+    F_r::FT,
+    p3::PSP3,
+    th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)),
+) where {FT}
+
+    leftpt = FT(1)
+    rightpt = FT(1e6)
+    radius = FT(0.8)
+
+    q_left = q_gamma(p3, F_r, N/N̂, log(leftpt), th)
+    q_right = q_gamma(p3, F_r, N/N̂, log(rightpt), th)
+
+    guess =
+        (q / q_left)^((log(rightpt) - log(leftpt)) / (log(q_right) - log(q_left)))
+
+    max = log(guess * exp(radius))
+    min = log(guess)
+
+    if guess < FT(2.5 * 1e4) #|| isequal(guess, NaN)
+        min = log(FT(20000))
+        max = log(FT(50000))
+    end
+    return (min = FT(min), max = FT(max))
+end
+
 """
     distrbution_parameter_solver()
 
@@ -344,20 +384,33 @@ function distribution_parameter_solver(
 
     # To ensure that λ is positive solve for x such that λ = exp(x)
     # We divide by N̂ to deal with large N₀ values for Float32
-    N̂ = FT(1e30)
+    N̂ = FT(1e20)
     shape_problem(x) = q / N̂ - q_gamma(p3, F_r, N / N̂, x, th)
 
+    # Get intial guess for solver 
+    (min, max) = get_bounds(N, N̂, q, F_r, p3, th)
+
     # Find slope parameter
+    sol = RS.find_zero(
+        shape_problem,
+        RS.SecantMethod(FT(min), FT(max)),
+        RS.CompactSolution(),
+        RS.RelativeSolutionTolerance(eps(FT)),                                                                                                                                                                                                                                                                                                                                                                  
+        10,
+    )
+    x = sol.root
+    converge = sol.converged
+    #= TO DO: switch back to this once testing is done and converged is not needed as a return 
     x =
         RS.find_zero(
             shape_problem,
-            RS.SecantMethod(FT(log(20000)), FT(log(50000))),
+            RS.SecantMethod(FT(log(4 * 1e4)), FT(log(5 * 1e4))),
             RS.CompactSolution(),
             RS.RelativeSolutionTolerance(eps(FT)),
             10,
         ).root
-
-    return (; λ = exp(x), N_0 = DSD_N₀(p3, N, exp(x)))
+    =#
+    return (; λ = exp(x), N_0 = DSD_N₀(p3, N, exp(x)), converged = converge)
 end
 
 end

test/P3SchemeSandbox.jl

@@ -0,0 +1,151 @@
+import CloudMicrophysics as CM
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.Parameters as CMP
+import CLIMAParameters as CP
+import Integrals as IN
+import SpecialFunctions as SF
+import RootSolvers as RS
+
+FT = Float64
+
+const PSP3 = CMP.ParametersP3
+p3 = CMP.ParametersP3(FT)
+
+"""
+    N′(D, p)
+
+ - D - maximum particle dimension
+ - p - a tuple containing N_0, λ, μ (intrcept, slope, and shape parameters for N′ respectively)
+
+ Returns the value of N′
+ Eq. 2 in Morrison and Milbrandt (2015).
+"""
+N′(D, p) = p.N_0 * exp(-p.λ * D)
+
+
+"""
+    q_helper(N_0, λ)
+
+ - p3 - a struct with P3 scheme parameters
+ - N_0 - intercept parameter of N′ gamma distribution
+ - λ - slope parameter of N′ gamma distribution
+ - F_r - rime mass fraction (q_rim/q_i)
+ - ρ_r - rime density (q_rim/B_rim) [kg/m^3]
+
+Returns the prognostic mass mixing ratio
+Eq. 5 in Morrison and Milbrandt (2015).
+"""
+function q_helper(p3::PSP3{FT}, N_0::FT, λ::FT, F_r::FT, ρ_r::FT) where {FT}
+    th = P3.thresholds(p3, ρ_r, F_r)
+    q′(D, p) = P3.p3_mass(p3, D, F_r, th) * N′(D, p)
+    problem = IN.IntegralProblem(q′, 0, Inf, (N_0 = N_0, λ = λ))
+    sol = IN.solve(problem, IN.HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
+    return FT(sol.u)
+end
+
+"""
+    lambda_expected()
+
+Returns the expected lambda to be compared to the non linear solver
+"""
+function λ_expected(q_i::FT, n_0::FT, ρ_i ::FT) where{FT}
+    m_e = FT(3)
+    r_0 = FT(1e-5)
+    ρ = FT(1)       # density of air
+    m_0 = FT(4/3 * π * ρ_i * r_0^3)
+    X_m = FT(1)
+    Δ_m = FT(0)
+
+    return ((SF.gamma(m_e + Δ_m + 1) * X_m * m_0 * n_0)/(q_i * ρ * r_0 ^ (m_e + Δ_m)))^(1/(m_e + Δ_m + 1))
+end
+
+function test_1M(q_i, n_0, ρ_i)
+
+    m_e = FT(3)
+    r_0 = FT(1e-5)
+    ρ = FT(1)       # density of air
+    m_0 = FT(4/3 * π * ρ_i * r_0^3)
+    X_m = FT(1)
+    Δ_m = FT(0)
+
+
+    shape_problem(λ) = q_i - 4/3 * π * ρ_i * n_0 * (6* λ^(-4))
+    #shape_problem(λ) = q_i - SF.gamma(m_e + Δ_m + 1) * X_m * m_0 * n_0  / (λ^(m_e+Δ_m+1)) / ρ / r_0^(m_e + Δ_m)
+
+    λ = RS.find_zero(
+        shape_problem,
+        RS.SecantMethod(FT(10), FT(10000)),
+        RS.CompactSolution(),
+        RS.RelativeSolutionTolerance(eps(FT)),
+        #100,
+    ).root
+
+    println(shape_problem(λ))
+
+    println("True λ = ", λ_expected(q_i, n_0, ρ_i))
+    println("Modeled λ = ", λ)
+end
+
+function test_solver(p3, q, N, ρ_r, F_r)
+    (λ, N_0) = P3.distribution_parameter_solver(p3, q, N, ρ_r, F_r)
+    println("λ solved = ", λ)
+    println("N_0 solved = ", N_0)
+end
+
+q_i = FT(1e-6)
+n_0 = FT(16 * 1e6)
+ρ_w = FT(1e3)
+ρ_i = FT(917)   # kg/m^3
+
+N = FT(1e8)
+q = FT(1e-3)
+q_r = FT(0.5 * 1e-6)
+B_r = FT(1/200 * 1e-4)
+n_0 = FT(16 * 1e6)
+ρ_w = FT(1e3)
+
+F_r = q_r/q_i
+ρ_r = ifelse(B_r == 0, FT(0), q_r/B_r)
+
+th = P3.thresholds(p3, ρ_r, F_r)
+
+
+#λ_ex = λ_expected(q_i, n_0, ρ_i)
+
+λ_ex = FT(25000)
+μ_ex = P3.μ_calc(λ_ex)
+n_0_ex = P3.N_0_helper(N, λ_ex, μ_ex)
+
+println(F_r, n_0_ex, log(λ_ex), μ_ex, th)
+q_calculated1 = P3.q_gamma(p3, F_r, n_0_ex, log(λ_ex), μ_ex, th)
+
+println(" ")
+println("λ_ex = ", λ_ex)
+println("n_0 expected = ", n_0_ex)
+println("q_calc = ", q_calculated1)
+
+#N = n_0/λ_ex
+#println("N_expected = ", N)
+    # N = n_0 * λ^(-0.00191*λ_ex^(0.8) + 2 - 1) * SF.gamma(1 + 0.00191*λ_ex^0.8 - 2)
+
+#q_calculated2 = P3.q_helper(p3, n_0, λ_ex, F_r, ρ_r)
+#N_calculated = P3.N_helper(n_0, λ_ex)
+
+#println("with gamma functions, q = ", q_calculated1)
+#println("with integral, q = ", q_calculated2)
+#println()
+
+#println("n_0 entered = ", n_0)
+#println("λ expected = ", λ_ex)
+
+println("testing 1M")
+test_1M(q_i, n_0, ρ_i)
+
+println(" ")
+println("Tetsing P3")
+test_solver(p3, q_calculated1, N, ρ_r, F_r)
+
+
+println("λ_ex = ", λ_ex)
+println("n_0 expected = ", n_0_ex)
+println("q_calc = ", q_calculated1)