From ca4357b04c11e519cc91c518c0b31fffeccfd247 Mon Sep 17 00:00:00 2001 From: Anastasia Popova Date: Wed, 21 Feb 2024 15:02:57 -0800 Subject: [PATCH] Add terminal velocities to P3 scheme --- Project.toml | 4 +- docs/Project.toml | 1 + docs/src/API.md | 1 - docs/src/P3Scheme.md | 29 ++ docs/src/plots/P3SchemePlots.jl | 72 ++++- docs/src/plots/P3TerminalVelocityPlots.jl | 140 +++++++++ src/Common.jl | 1 + src/P3Scheme.jl | 358 +++++++++++++++++----- test/Project.toml | 2 + test/p3_tests.jl | 138 +++++---- test/performance_tests.jl | 16 + 11 files changed, 605 insertions(+), 157 deletions(-) create mode 100644 docs/src/plots/P3TerminalVelocityPlots.jl diff --git a/Project.toml b/Project.toml index 7e74a80cd..4571258b3 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.18.0" ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" @@ -22,13 +23,14 @@ CloudyExt = "Cloudy" EmulatorModelsExt = ["DataFrames", "MLJ"] [compat] -ClimaParams = "0.10" +ClimaParams = "0.10.5" Cloudy = "0.4" DataFrames = "1.6" DocStringExtensions = "0.8, 0.9" EnsembleKalmanProcesses = "1.1.5" ForwardDiff = "0.10" MLJ = "0.20" +QuadGK = "2.9.4" RootSolvers = "0.3, 0.4" SpecialFunctions = "1, 2" Thermodynamics = "0.12.4" diff --git a/docs/Project.toml b/docs/Project.toml index 76d874ac0..10b38097f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,6 +9,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/docs/src/API.md b/docs/src/API.md index c74418078..f7b922397 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -70,7 +70,6 @@ MicrophysicsFlexible.weighted_vt ```@docs P3Scheme P3Scheme.thresholds -P3Scheme.p3_mass P3Scheme.distribution_parameter_solver ``` diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index aa6850448..a197f080f 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -159,3 +159,32 @@ Using this approach we get the following relative errors for ``\lambda`` include("plots/P3LambdaErrorPlots.jl") ``` ![](P3LambdaHeatmap.png) + +## Terminal Velocity + +We use the [Chen2022](@cite) velocity parametrization: +```math +V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} +``` +where ``\phi = (16 \rho_{ice}^2 A(D)^3) / (9 \pi m(D)^2)`` is the aspect ratio, +and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters. + +The mass-weighted fall speed (``V_m``) and the number-weighted fall speed (``V_n``) are calculated as +```math +V_m = \frac{\int_{0}^{\infty} \! V(D) m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D} +``` +```math +V_n = \frac{\int_{0}^{\infty} \! V(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! N'(D) \mathrm{d}D} +``` + +We also plot the mass-weighted mean particle size ``D_m`` which is given by: +```math +D_m = \frac{\int_{0}^{\infty} \! D m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D} +``` + +Below w show these relationships for small, medium, and large ``D_m`` +They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite). +```@example +include("plots/P3TerminalVelocityPlots.jl") +``` +![](MorrisonandMilbrandtFig2.svg) diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index 695c182f7..c106b8f99 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -9,6 +9,54 @@ 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) @@ -43,17 +91,13 @@ function area( # Area regime: if P3.D_th_helper(p3) > D return A_s(D) # small spherical ice - end - if F_r == 0 + elseif F_r == 0 return A_ns(p3, D) # large nonspherical unrimed ice - end - if th.D_gr > D >= P3.D_th_helper(p3) + elseif th.D_gr > D >= P3.D_th_helper(p3) return A_ns(p3, D) # dense nonspherical ice - end - if th.D_cr > D >= th.D_gr + elseif th.D_cr > D >= th.D_gr return A_s(D) # graupel - end - if D >= th.D_cr + elseif D >= th.D_cr return A_r(p3, F_r, D) # partially rimed ice end @@ -99,9 +143,9 @@ 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.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) + 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) # 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) @@ -120,9 +164,9 @@ 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.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) + 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) # 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) diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl new file mode 100644 index 000000000..ce3f2c3ce --- /dev/null +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -0,0 +1,140 @@ +import ClimaParams +import CloudMicrophysics as CM +import CloudMicrophysics.P3Scheme as P3 +import CloudMicrophysics.Parameters as CMP +import CairoMakie as Plt + +const PSP3 = CMP.ParametersP3 + +FT = Float64 + +p3 = CMP.ParametersP3(FT) + +function get_values( + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + q::FT, + N::FT, + ρ_a::FT, + x_resolution::Int, + y_resolution::Int, +) where {FT} + F_rs = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) + ρ_rs = range(FT(25), stop = FT(975), length = y_resolution) + + V_m = zeros(x_resolution, y_resolution) + D_m = zeros(x_resolution, y_resolution) + + for i in 1:x_resolution + for j in 1:y_resolution + F_r = F_rs[i] + ρ_r = ρ_rs[j] + + V_m[i, j] = + P3.terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a)[2] + # get D_m in mm for plots + D_m[i, j] = 1e3 * P3.D_m(p3, q, N, ρ_r, F_r) + end + end + return (; F_rs, ρ_rs, V_m, D_m) +end + +function make_axis_top(fig, col, title) + return Plt.Axis( + fig[1, col], + xlabel = "F_r", + ylabel = "ρ_r", + title = title, + width = 350, + height = 350, + limits = (0, 1.0, 25, 975), + xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], + yticks = [200, 400, 600, 800], + ) +end +function make_axis_bottom(fig, col, title) + return Plt.Axis( + fig[3, col], + height = 350, + width = 350, + xlabel = "F_r", + ylabel = "ρ_r", + title = title, + ) +end + +#! format: off +function figure_2() + Chen2022 = CMP.Chen2022VelType(FT) + # density of air in kg/m^3 + ρ_a = FT(1.2) #FT(1.293) + # small D_m + q_s = FT(0.0008) + N_s = FT(1e6) + # medium D_m + q_m = FT(0.22) + N_m = FT(1e6) + # large D_m + q_l = FT(0.7) + N_l = FT(1e6) + # get V_m and D_m + xres = 100 + yres = 100 + (F_rs, ρ_rs, V_ms, D_ms) = get_values(p3, Chen2022.snow_ice, q_s, N_s, ρ_a, xres, yres) + (F_rm, ρ_rm, V_mm, D_mm) = get_values(p3, Chen2022.snow_ice, q_m, N_m, ρ_a, xres, yres) + (F_rl, ρ_rl, V_ml, D_ml) = get_values(p3, Chen2022.snow_ice, q_l, N_l, ρ_a, xres, yres) + + fig = Plt.Figure() + + # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015 + + args = (color = :black, labels = true, levels = 3, linewidth = 1.5, labelsize = 18) + + ax1 = make_axis_top(fig, 1, "Small Dₘ") + hm = Plt.contourf!(ax1, F_rs, ρ_rs, V_ms, levels = 0.402:0.001:0.413) + Plt.contour!(ax1, F_rs, ρ_rs, D_ms; args...) + Plt.Colorbar(fig[2, 1], hm, vertical = false) + + ax2 = make_axis_top(fig, 2, "Medium Dₘ") + hm = Plt.contourf!(ax2, F_rm, ρ_rm, V_mm, levels = 2:0.1:4) + Plt.contour!(ax2, F_rm, ρ_rm, D_mm; args...) + Plt.Colorbar(fig[2, 2], hm, vertical = false) + + ax3 = make_axis_top(fig, 3, "Large Dₘ") + hm = Plt.contourf!(ax3, F_rl, ρ_rl, V_ml, levels = 2:0.2:5) + Plt.contour!(ax3, F_rl, ρ_rl, D_ml; args...) + Plt.Colorbar(fig[2, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax2) + Plt.linkxaxes!(ax2, ax3) + Plt.linkyaxes!(ax1, ax2) + Plt.linkyaxes!(ax2, ax3) + + ## Plot D_m as second row of comparisons + + ax4 = make_axis_bottom(fig, 1, "Small Dₘ vs F_r and ρ_r") + hm = Plt.contourf!(ax4, F_rs, ρ_rs, D_ms) + Plt.Colorbar(fig[4, 1], hm, vertical = false) + + ax5 = make_axis_bottom(fig, 2, "Medium Dₘ vs F_r and ρ_r") + hm = Plt.contourf!(ax5, F_rm, ρ_rm, D_mm) + Plt.Colorbar(fig[4, 2], hm, vertical = false) + + ax6 = make_axis_bottom(fig, 3, "Large Dₘ vs F_r and ρ_r") + hm = Plt.contourf!(ax6, F_rl, ρ_rl, D_ml) + Plt.Colorbar(fig[4, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax4) + Plt.linkxaxes!(ax4, ax5) + Plt.linkxaxes!(ax5, ax6) + Plt.linkyaxes!(ax1, ax4) + Plt.linkyaxes!(ax4, ax5) + Plt.linkyaxes!(ax5, ax6) + + Plt.resize_to_layout!(fig) + Plt.save("MorrisonandMilbrandtFig2.svg", fig) +end +#! format: on + +# Terminal Velocity figure +figure_2() diff --git a/src/Common.jl b/src/Common.jl index 15f09e861..856d5acf5 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -18,6 +18,7 @@ export a_w_eT export a_w_ice export Chen2022_vel_add export Chen2022_vel_coeffs_small +export Chen2022_vel_coeffs_large """ G_func(air_props, tps, T, Liquid()) diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl index 440baa37e..6992b0b82 100644 --- a/src/P3Scheme.jl +++ b/src/P3Scheme.jl @@ -12,14 +12,16 @@ Note: Particle size is defined as its maximum length (i.e. max dimesion). module P3Scheme import SpecialFunctions as SF - +import QuadGK as QGK import RootSolvers as RS + import ClimaParams as CP import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.Common as CO const PSP3 = CMP.ParametersP3 -export thresholds, p3_mass, distribution_parameter_solver +export thresholds, distribution_parameter_solver """ α_va_si(p3) @@ -154,65 +156,6 @@ 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} - 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 - - # TODO - would something like this be better? - #return ifelse(D_th_helper(p3) > D, mass_s(D, p3.ρ_i), - # ifelse(F_r == 0, mass_nl(p3, D), - # ifelse(th.D_gr > D >= D_th_helper(p3), mass_nl(p3, D), - # ifelse(th.D_cr > D >= th.D_gr, mass_s(D, th.ρ_g), - # mass_r(p3, D, F_r))))) - -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)) @@ -255,7 +198,7 @@ function DSD_μ_approx(p3::PSP3, q::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} # Return approximation between them μ = (p3.μ_max / (q_over_N_max - q_over_N_min)) * (log(q / N) - q_over_N_min) - # Clip approximation between 0 and 6 + # Clip approximation between 0 and 6 return min(p3.μ_max, max(FT(0), μ)) end @@ -288,6 +231,54 @@ function DSD_N₀(p3::PSP3, N::FT, λ::FT) where {FT} return N / Γ(1 + μ) * λ^(1 + μ) end +""" + ∫_Γ(x₀, x_end, c1, c2, c3) + + - x₀ - lower bound + - x_end - upper bound + - c1, c2, c3 - respective constants + +f(D, c1, c2, c3) = c1 * D ^ (c2) * exp(-c3 * D) + +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 == Inf + return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3)) + elseif x₀ == 0 + return c1 * c3^(-c2 - 1) * (Γ(1 + c2) - Γ(1 + c2, x_end * c3)) + else + return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3) - Γ(1 + c2, x_end * c3)) + end +end + +""" + ∫_Γ(x₀, xₘ, x_end, c1, c2, c3, c4, c5, c6) + + - x₀ - lower bound + - xₘ - switch point + - x_end - upper bound + - c1, c2, c3 - respective constants for the first part of the integral + - c4, c5, c6 - respective constants for the second part of the integral + +f(D, c1, c2, c3) = c1 * D ^ (c2) * exp(-c3 * D) + +Integrates f(D, c1, c2, c3) dD from x₀ to xₘ and f(D, c4, c5, c6) dD from xₘ to x_end +""" +function ∫_Γ( + x₀::FT, + xₘ::FT, + x_end::FT, + c1::FT, + c2::FT, + c3::FT, + c4::FT, + c5::FT, + c6::FT, +) where {FT} + return ∫_Γ(x₀, xₘ, c1, c2, c3) + ∫_Γ(xₘ, x_end, c4, c5, c6) +end + """ q_(p3, ρ, F_r, λ, μ, D_min, D_max) @@ -306,25 +297,20 @@ end # or # q_rim > 0 and D_min = D_gr, D_max = D_cr, ρ = ρ_g function q_s(p3::PSP3, ρ::FT, μ::FT, λ::FT, D_min::FT, D_max::FT) where {FT} - x = μ + 4 - return FT(π) / 6 * ρ / λ^x * (Γ(x, λ * D_min) - Γ(x, λ * D_max)) + return ∫_Γ(D_min, D_max, FT(π) / 6 * ρ, μ + 3, λ) end # q_rim = 0 and D_min = D_th, D_max = inf function q_rz(p3::PSP3, μ::FT, λ::FT, D_min::FT) where {FT} - x = μ + p3.β_va + 1 - return α_va_si(p3) / λ^x * (Γ(x) + Γ(x, λ * D_min) - (x - 1) * Γ(x - 1)) + return ∫_Γ(D_min, FT(Inf), α_va_si(p3), μ + p3.β_va, λ) end # q_rim > 0 and D_min = D_th and D_max = D_gr function q_n(p3::PSP3, μ::FT, λ::FT, D_min::FT, D_max::FT) where {FT} - x = μ + p3.β_va + 1 - return α_va_si(p3) / λ^x * (Γ(x, λ * D_min) - Γ(x, λ * D_max)) + return ∫_Γ(D_min, D_max, α_va_si(p3), μ + p3.β_va, λ) end # partially rimed ice or large unrimed ice (upper bound on D is infinity) # q_rim > 0 and D_min = D_cr, D_max = inf function q_r(p3::PSP3, F_r::FT, μ::FT, λ::FT, D_min::FT) where {FT} - x = μ + p3.β_va + 1 - return α_va_si(p3) / (1 - F_r) / λ^x * - (Γ(x) + Γ(x, λ * D_min) - (x - 1) * Γ(x - 1)) + return ∫_Γ(D_min, FT(Inf), α_va_si(p3) / (1 - F_r), μ + p3.β_va, λ) end """ @@ -435,13 +421,13 @@ function distribution_parameter_solver( # Get the thresholds for different particles regimes th = thresholds(p3, ρ_r, F_r) - # Get μ given q and N + # Get μ given q and N μ = DSD_μ_approx(p3, q, N, ρ_r, F_r) # To ensure that λ is positive solve for x such that λ = exp(x) shape_problem(x) = q / N - q_over_N_gamma(p3, F_r, x, μ, th) - # Get intial guess for solver + # Get intial guess for solver (; min, max) = get_bounds(N, q, μ, F_r, p3, th) # Find slope parameter @@ -451,10 +437,234 @@ function distribution_parameter_solver( RS.SecantMethod(min, max), RS.CompactSolution(), RS.RelativeSolutionTolerance(eps(FT)), - 10, + 5, ).root return (; λ = exp(x), N_0 = DSD_N₀(p3, N, exp(x))) end +""" + terminal_velocity(p3, Chen2022, q, N, ρ_r, F_r, ρ_a) + + - p3 - a struct with P3 scheme parameters + - Chen2022 - a struch with terminal velocity parameters as in Chen(2022) + - q - mass mixing ratio + - N - number mixing ratio + - ρ_r - rime density (q_rim/B_rim) [kg/m^3] + - F_r - rime mass fraction (q_rim/q_i) + - ρ_a - density of air + + Returns the mass and number weighted fall speeds + Eq C10 of Morrison and Milbrandt (2015) +""" +function terminal_velocity( + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + q::FT, + N::FT, + ρ_r::FT, + F_r::FT, + ρ_a::FT, +) where {FT} + # Get the pree parameters for terminal velocities of small + # and large particles + small = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + large = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + get_p(prs, it) = (prs[1][it], prs[2][it], prs[3][it]) + + # 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 + + # Get the shape parameters of the particle size distribution + (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) + μ = DSD_μ(p3, λ) + + # TODO: Change when each value used depending on type of particle + # TODO: or keep fixed and add to ClimaParams...? + κ = FT(-1 / 6) #FT(1/3) + # Redefine α_va to be in si units + α_va = α_va_si(p3) + + aₛ(a) = a * N_0 + bₛ(b) = b + μ + cₛ(c) = c + λ + + aₛ_m(a) = aₛ(a) * FT(π) / 6 * p3.ρ_i + bₛ_m(b) = bₛ(b) + 3 + + spheres_n(a, b, c) = (aₛ(a), bₛ(b), cₛ(c)) + spheres_m(a, b, c) = (aₛ_m(a), bₛ_m(b), cₛ(c)) + + aₙₛ(a) = aₛ(a) * (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ + bₙₛ(b) = bₛ(b) + κ * (3 * p3.σ - 2 * p3.β_va) + + aₙₛ_m(a) = aₙₛ(a) * α_va + bₙₛ_m(b) = bₙₛ(b) + p3.β_va + + non_spheres_n(a, b, c) = (aₙₛ(a), bₙₛ(b), cₛ(c)) + non_spheres_m(a, b, c) = (aₙₛ_m(a), bₙₛ_m(b), cₛ(c)) + + aᵣₛ(a) = aₛ(a) * (p3.ρ_i / ρ_g)^(2 * κ) + aᵣₛ_m(a) = aᵣₛ(a) * FT(π) / 6 * ρ_g + + rimed_n(a, b, c) = (aᵣₛ(a), bₛ(b), cₛ(c)) + rimed_m(a, b, c) = (aᵣₛ_m(a), bₛ_m(b), cₛ(c)) + + v_n_D_cr(D, a, b, c) = + a * + N_0 * + D^(b + μ) * + exp((-c - λ) * D) * + ( + 16 * p3.ρ_i^2 * (F_r * π / 4 * D^2 + (1 - F_r) * p3.γ * D^p3.σ)^3 / + (9 * π * (α_va / (1 - F_r) * D^p3.β_va)^2) + )^κ + v_m_D_cr(D, a, b, c) = v_n_D_cr(D, a, b, c) * (α_va / (1 - F_r) * D^p3.β_va) + + v_m = 0 + v_n = 0 + for i in 1:2 + if F_r == 0 + v_m += ∫_Γ(FT(0), D_th, spheres_m(get_p(small, i)...)...) + v_n += ∫_Γ(FT(0), D_th, spheres_n(get_p(small, i)...)...) + + v_m += ∫_Γ( + D_th, + D_ct, + Inf, + non_spheres_m(get_p(small, i)...)..., + non_spheres_m(get_p(large, i)...)..., + ) + v_n += ∫_Γ( + D_th, + D_ct, + Inf, + non_spheres_n(get_p(small, i)...)..., + non_spheres_n(get_p(large, i)...)..., + ) + else + # Velocity coefficients for small particles + v_m += ∫_Γ(FT(0), D_th, spheres_m(get_p(small, i)...)...) + v_n += ∫_Γ(FT(0), D_th, spheres_n(get_p(small, i)...)...) + is_large = false + + # D_th to D_gr + if !is_large && D_gr > D_ct + v_m += ∫_Γ( + D_th, + D_ct, + D_gr, + non_spheres_m(get_p(small, i)...)..., + non_spheres_m(get_p(large, i)...)..., + ) + v_n += ∫_Γ( + D_th, + D_ct, + D_gr, + non_spheres_n(get_p(small, i)...)..., + non_spheres_n(get_p(large, i)...)..., + ) + # Switch to large particles + is_large = true + else + v_m += ∫_Γ(D_th, D_gr, non_spheres_m(get_p(small, i)...)...) + v_n += ∫_Γ(D_th, D_gr, non_spheres_n(get_p(small, i)...)...) + end + + # D_gr to D_cr + if !is_large && D_cr > D_ct + v_m += ∫_Γ( + D_gr, + D_ct, + D_cr, + rimed_m(get_p(small, i)...)..., + rimed_m(get_p(large, i)...)..., + ) + v_n += ∫_Γ( + D_gr, + D_ct, + D_cr, + rimed_n(get_p(small, i)...)..., + rimed_n(get_p(large, i)...)..., + ) + # Switch to large particles + is_large = true + elseif is_large + v_m += ∫_Γ(D_gr, D_cr, rimed_m(get_p(large, i)...)...) + v_n += ∫_Γ(D_gr, D_cr, rimed_n(get_p(large, i)...)...) + else + v_m += ∫_Γ(D_gr, D_cr, rimed_m(get_p(small, i)...)...) + v_n += ∫_Γ(D_gr, D_cr, rimed_n(get_p(small, i)...)...) + end + + # D_cr to Infinity + if !is_large + (Im, em) = + QGK.quadgk(D -> v_m_D_cr(D, get_p(small, i)...), D_cr, D_ct) + (In, en) = + QGK.quadgk(D -> v_n_D_cr(D, get_p(small, i)...), D_cr, D_ct) + v_m += Im + v_n += In + + # Switch to large particles + (Im, em) = + QGK.quadgk(D -> v_m_D_cr(D, get_p(large, i)...), D_ct, Inf) + (In, en) = + QGK.quadgk(D -> v_n_D_cr(D, get_p(large, i)...), D_ct, Inf) + v_m += Im + v_n += In + else + # TODO - check if it should be large or small + (Im, em) = + QGK.quadgk(D -> v_m_D_cr(D, get_p(large, i)...), D_cr, Inf) + (In, en) = + QGK.quadgk(D -> v_n_D_cr(D, get_p(large, i)...), D_cr, Inf) + v_m += Im + v_n += In + end + end + end + return (v_n / N, v_m / q) +end + +""" + D_m (p3, q, N, ρ_r, F_r) + + - p3 - a struct with P3 scheme parameters + - q - mass mixing ratio + - N - number mixing ratio + - ρ_r - rime density (q_rim/B_rim) [kg/m^3] + - F_r - rime mass fraction (q_rim/q_i) + + Return the mass weighted mean particle size [m] +""" +function D_m(p3::PSP3, q::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} + # Get the thresholds for different particles regimes + th = thresholds(p3, ρ_r, F_r) + D_th = D_th_helper(p3) + + # Get the shape parameters + (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) + μ = DSD_μ(p3, λ) + + # Redefine α_va to be in si units + α_va = α_va_si(p3) + + # Calculate numerator + n = 0 + if F_r == 0 + n += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n += ∫_Γ(D_th, Inf, α_va * N_0, μ + p3.β_va + 1, λ) + else + n += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n += ∫_Γ(D_th, th.D_gr, α_va * N_0, μ + p3.β_va + 1, λ) + n += ∫_Γ(th.D_gr, th.D_cr, π / 6 * th.ρ_g * N_0, μ + 4, λ) + n += ∫_Γ(th.D_cr, Inf, α_va / (1 - F_r) * N_0, μ + p3.β_va + 1, λ) + end + # Normalize by q + return n / q +end + end diff --git a/test/Project.toml b/test/Project.toml index 30b163471..5e6ba8d83 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -23,6 +23,7 @@ MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -30,5 +31,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] +ClimaParams = "0.10.5" KernelAbstractions = "0.9" Optim = "<1.9.3" diff --git a/test/p3_tests.jl b/test/p3_tests.jl index 14c38aec4..e2dab5a4f 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -85,71 +85,6 @@ 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) #TODO - this is very large for eps - - for ρ in ρs - for F_r in F_rs - D_th = P3.D_th_helper(p3) - D1 = D_th / 2 - th = P3.thresholds(p3, ρ, F_r) - - 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, th) == - P3.mass_s(D1, p3.ρ_i) - TT.@test P3.p3_mass(p3, D2, F_r, th) == P3.mass_nl(p3, D2) - end - - end - end - end - -end - function test_p3_shape_solver(FT) p3 = CMP.ParametersP3(FT) @@ -200,14 +135,83 @@ function test_p3_shape_solver(FT) end end +function test_velocities(FT) + Chen2022 = CMP.Chen2022VelType(FT) + p3 = CMP.ParametersP3(FT) + q = FT(0.22) + N = FT(1e6) + ρ_a = FT(1.2) + ρ_rs = [FT(200), FT(400), FT(600), FT(800)] + F_rs = [FT(0), FT(0.2), FT(0.4), FT(0.6), FT(0.8)] + + TT.@testset "Mass and number weighted terminal velocities" begin + paper_vals = [ + [1.5, 1.5, 1.5, 1.5, 1.5], + [1.5, 1.5, 2.5, 2.5, 2.5], + [1.5, 2.5, 2.5, 2.5, 2.5], + [1.5, 2.5, 3.5, 3.5, 3.5], + ] + expected_vals = [ + [1.52, 1.46, 1.41, 1.36, 1.24], + [1.52, 1.47, 1.44, 1.42, 1.35], + [1.52, 1.47, 1.45, 1.44, 1.42], + [1.52, 1.47, 1.45, 1.45, 1.45], + ] + for i in 1:length(ρ_rs) + for j in 1:length(F_rs) + ρ_r = ρ_rs[i] + F_r = F_rs[j] + + calculated_vel = P3.terminal_velocity( + p3, + Chen2022.snow_ice, + q, + N, + ρ_r, + F_r, + ρ_a, + ) + + # number weighted + TT.@test calculated_vel[2] > 0 + TT.@test expected_vals[i][j] ≈ calculated_vel[1] atol = 0.1 + + # mass weighted + TT.@test calculated_vel[1] > 0 + TT.@test paper_vals[i][j] ≈ calculated_vel[2] atol = 3.14 + end + end + end + + TT.@testset "Mass-weighted mean diameters" begin + paper_vals = [ + [5, 5, 5, 5, 5], + [4.5, 4.5, 4.5, 4.5, 4.5], + [3.5, 3.5, 3.5, 3.5, 3.5], + [3.5, 3.5, 2.5, 2.5, 2.5], + ] + for i in 1:length(ρ_rs) + for j in 1:length(F_rs) + ρ_r = ρ_rs[i] + F_r = F_rs[j] + + calculated_dm = P3.D_m(p3, q, N, ρ_r, F_r) * 1e3 + + TT.@test calculated_dm > 0 + TT.@test paper_vals[i][j] ≈ calculated_dm atol = 3.14 + + end + end + end +end + println("Testing Float32") test_p3_thresholds(Float32) -test_p3_mass(Float32) #TODO - only works for Float64 now. We should switch the units inside the solver # from SI base to something more managable #test_p3_shape_solver(Float32) println("Testing Float64") test_p3_thresholds(Float64) -test_p3_mass(Float64) test_p3_shape_solver(Float64) +test_velocities(Float64) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 2b78dea7a..6090889ad 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -91,6 +91,7 @@ function benchmark_test(FT) ρ_r = FT(400.0) F_r = FT(0.95) + N = FT(1e8) T_air_2 = FT(250) T_air_cold = FT(230) @@ -126,6 +127,21 @@ function benchmark_test(FT) # P3 scheme bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80) + if FT == Float64 + bench_press( + P3.distribution_parameter_solver, + (p3, q_ice, N, ρ_r, F_r), + 1e5, + ) + bench_press( + P3.terminal_velocity, + (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air), + 2e5, + 3e4, + 2e3, + ) + bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5) + end # aerosol activation bench_press(