diff --git a/Project.toml b/Project.toml index 799f65811..a3c15ece7 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" @@ -29,6 +30,7 @@ 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 0711e8470..56125d67b 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -67,7 +67,6 @@ MicrophysicsFlexible.weighted_vt ```@docs P3Scheme P3Scheme.thresholds -P3Scheme.p3_mass P3Scheme.distribution_parameter_solver ``` diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index 695c182f7..f80093d6d 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 @@ -76,6 +120,69 @@ function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect) aspect = aspect, limits = ((0.02, 10.0), nothing), ) + sol_5 = P3.thresholds(p3, 400.0, 0.5) + sol_8 = P3.thresholds(p3, 400.0, 0.8) + + #! format: off + fig1_a_0 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) + fig1_a_5 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.5, sol_5) for D in D_range], color = cl[2], linewidth = lw) + fig1_a_8 = CMK.lines!(ax1_a, D_range * 1e3, [p3_mass(p3, D, 0.8, sol_8) for D in D_range], color = cl[3], linewidth = lw) + + d_tha = CMK.vlines!(ax1_a, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) + d_cr_5 = CMK.vlines!(ax1_a, sol_5[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) + d_cr_8 = CMK.vlines!(ax1_a, sol_8[1] * 1e3, linestyle = :dot, color = cl[3], linewidth = lw) + d_gr_5 = CMK.vlines!(ax1_a, sol_5[2] * 1e3, linestyle = :dash, color = cl[2], linewidth = lw) + d_gr_8 = CMK.vlines!(ax1_a, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) + + leg1_a = CMK.Legend(fig1_a[8:9, 1], [fig1_a_0, fig1_a_5, fig1_a_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) + leg1_a_dth = CMK.Legend(fig1_a[8:9, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) + leg1_a_dcr = CMK.Legend(fig1_a[8:9, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) + leg1_a_dgr = CMK.Legend(fig1_a[8:9, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) + + #! format: on + CMK.save("MorrisonandMilbrandtFig1a.svg", fig1_a) + + fig1_b = CMK.Figure() + ax1_b = CMK.Axis( + fig1_b[1:10, 1:11], + title = CMK.L"m(D) regime for $F_r = 0.95$", + xlabel = CMK.L"$D$ (mm)", + ylabel = CMK.L"$m$ (kg)", + xscale = CMK.log10, + yscale = CMK.log10, + yminorticksvisible = true, + yminorticks = CMK.IntervalsBetween(3), + xminorticksvisible = true, + xminorticks = CMK.IntervalsBetween(5), + xticks = [0.01, 0.1, 1, 10], + aspect = 1.67, + limits = ((0.02, 10.0), nothing), + ) + + sol_2 = P3.thresholds(p3, 200.0, 0.95) + sol_4 = P3.thresholds(p3, 400.0, 0.95) + sol_8 = P3.thresholds(p3, 800.0, 0.95) + + #! format: off + fig1_b200 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig1_b400 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig1_b800 = CMK.lines!(ax1_b, D_range * 1e3, [p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw) + + d_thb = CMK.vlines!(ax1_b, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw) + d_cr_200 = CMK.vlines!(ax1_b, sol_2[1] * 1e3, linestyle = :dot, color = cl[1], linewidth = lw) + d_cr_400 = CMK.vlines!(ax1_b, sol_4[1] * 1e3, linestyle = :dot, color = cl[2], linewidth = lw) + d_cr_800 = CMK.vlines!(ax1_b, sol_8[1] * 1e3, linestyle = :dot, color = cl[3], linewidth = lw) + d_gr_200 = CMK.vlines!(ax1_b, sol_2[2] * 1e3, linestyle = :dash, color = cl[1], linewidth = lw) + d_gr_400 = CMK.vlines!(ax1_b, sol_4[2] * 1e3, linestyle = :dash, color = cl[2], linewidth = lw) + d_gr_800 = CMK.vlines!(ax1_b, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) + + leg1_b = CMK.Legend(fig1_b[11:12, 4], [fig1_b200, fig1_b400, fig1_b800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + leg1_b_dth = CMK.Legend(fig1_b[11:12, 5], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) + leg1_b_dcr = CMK.Legend(fig1_b[11:12, 8], [d_cr_200, d_cr_400, d_cr_800], [CMK.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + leg1_b_dgr = CMK.Legend(fig1_b[11:12, 6], [d_gr_200, d_gr_400, d_gr_800], [CMK.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) + + #! format: on + CMK.save("MorrisonandMilbrandtFig1b.svg", fig1_b) end #! format: off @@ -99,9 +206,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 +227,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..71fcabf15 --- /dev/null +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -0,0 +1,255 @@ +import CLIMAParameters +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_mass(p3, Chen2022, q, N, ρ_r, F_r, ρ_a) + # 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 figure_2() + Chen2022 = CMP.Chen2022VelType(FT) + # density of air in kg/m^3 + ρ_a = FT(1.2) #FT(1.293) + + f = Plt.Figure() + xres = 100 + yres = 100 + min = FT(0) + max = FT(10) + + # 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 + (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) + + ax1 = Plt.Axis( + f[1, 1], + xlabel = "F_r", + ylabel = "ρ_r", + title = "Small Dₘ", + 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], + ) + + Plt.contourf!( + F_rs, + ρ_rs, + V_ms, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.contour!(F_rs, ρ_rs, D_ms, color = :black, labels = true, levels = 3) + Plt.Colorbar( + f[2, 1], + limits = ( + minimum(v for v in V_ms if !isnan(v)), + maximum(v for v in V_ms if !isnan(v)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) + + ax2 = Plt.Axis( + f[1, 2], + xlabel = "F_r", + ylabel = "ρ_r", + title = "Medium Dₘ", + 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], + ) + + Plt.contourf!( + F_rm, + ρ_rm, + V_mm, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.contour!(F_rm, ρ_rm, D_mm, color = :black, labels = true, levels = 3) + Plt.Colorbar( + f[2, 2], + limits = ( + minimum(v for v in V_mm if !isnan(v)), + maximum(v for v in V_mm if !isnan(v)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) + + ax3 = Plt.Axis( + f[1, 3], + xlabel = "F_r", + ylabel = "ρ_r", + title = "Large Dₘ", + 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], + ) + + Plt.contourf!( + F_rl, + ρ_rl, + V_ml, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.contour!(F_rl, ρ_rl, D_ml, color = :black, labels = true, levels = 3) + Plt.Colorbar( + f[2, 3], + limits = ( + minimum(v for v in V_ml if !isnan(v)), + maximum(v for v in V_ml if !isnan(v)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + 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 = Plt.Axis( + f[3, 1], + height = 350, + width = 350, + xlabel = "F_r", + ylabel = "ρ_r", + title = "Small Dₘ vs F_r and ρ_r", + ) + Plt.contourf!( + F_rs, + ρ_rs, + D_ms, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.Colorbar( + f[4, 1], + limits = ( + minimum(d for d in D_ms if !isnan(d)), + maximum(d for d in D_ms if !isnan(d)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) + ax5 = Plt.Axis( + f[3, 2], + height = 350, + width = 350, + xlabel = "F_r", + ylabel = "ρ_r", + title = "Medium Dₘ vs F_r and ρ_r", + ) + Plt.contourf!( + F_rm, + ρ_rm, + D_mm, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.Colorbar( + f[4, 2], + limits = ( + minimum(d for d in D_mm if !isnan(d)), + maximum(d for d in D_mm if !isnan(d)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + vertical = false, + ) + ax6 = Plt.Axis( + f[3, 3], + height = 350, + width = 350, + xlabel = "F_r", + ylabel = "ρ_r", + title = "Large Dₘ vs F_r and ρ_r", + ) + Plt.contourf!( + F_rl, + ρ_rl, + D_ml, + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + ) + Plt.Colorbar( + f[4, 3], + limits = ( + minimum(d for d in D_ml if !isnan(d)), + maximum(d for d in D_ml if !isnan(d)), + ), + colormap = Plt.reverse(Plt.cgrad(:Spectral)), + flipaxis = true, + 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!(f) + Plt.save("MorrisonandMilbrandtFig2.svg", f) +end + +# Terminal Velocity figure +figure_2() diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl index 440baa37e..a42de1647 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)) @@ -288,6 +231,26 @@ function DSD_N₀(p3::PSP3, N::FT, λ::FT) where {FT} return N / Γ(1 + μ) * λ^(1 + μ) end +""" + integrate_to_gamma(a, b, c1, c2, c3) + + - a - lower bound + - b - upper bound + - c1, c2, c3 - respective constants + + Integrates the function c1 * D ^ (c2) * exp(-c3 * D) dD from a to b + Returns the result +""" +function integrate_to_gamma(a::FT, b::FT, c1::FT, c2::FT, c3::FT) where {FT} + if b == Inf + return c1 * c3^(-c2 - 1) * (Γ(1 + c2, a * c3)) + elseif a == 0 + return c1 * c3^(-c2 - 1) * (Γ(1 + c2) - Γ(1 + c2, b * c3)) + else + return c1 * c3^(-c2 - 1) * (Γ(1 + c2, a * c3) - Γ(1 + c2, b * c3)) + end +end + """ q_(p3, ρ, F_r, λ, μ, D_min, D_max) @@ -306,25 +269,26 @@ 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 integrate_to_gamma(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 integrate_to_gamma(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 integrate_to_gamma(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 integrate_to_gamma( + D_min, + FT(Inf), + α_va_si(p3) / (1 - F_r), + μ + p3.β_va, + λ, + ) end """ @@ -451,10 +415,447 @@ 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_mass(p3, Chen2022, q, N, ρ_r, F_r) + + - p3 - a struct with P3 scheme parameters + - Chen 2022 - 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 (total)-weighted fall speed + Eq C10 of Morrison and Milbrandt (2015) +""" +function terminal_velocity_mass( + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + q::FT, + N::FT, + ρ_r::FT, + F_r::FT, + ρ_a::FT, +) where {FT} + + # Get the thresholds for different particles regimes + th = thresholds(p3, ρ_r, F_r) + D_th = D_th_helper(p3) + cutoff = FT(0.000625) # TODO add to the struct + + # Get the shape parameters + (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) + μ = DSD_μ(p3, λ) + + # TO DO: Change when each value used depending on type of particle + κ = FT(-1 / 6) #FT(1/3) + + # Redefine α_va to be in si units + α_va = α_va_si(p3) + + v = 0 + for i in 1:2 + if F_r == 0 + # Velocity coefficients for small particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + v += integrate_to_gamma( + FT(0), + D_th, + FT(π) / 6 * p3.ρ_i * ai[i] * N_0, + bi[i] + μ + 3, + ci[i] + λ, + ) + v += integrate_to_gamma( + D_th, + cutoff, + α_va * + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + + # Get velocity coefficients for large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + v += integrate_to_gamma( + cutoff, + Inf, + α_va * + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + else + # Velocity coefficients for small particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + large = false + + v += integrate_to_gamma( + FT(0), + D_th, + FT(π) / 6 * p3.ρ_i * ai[i] * N_0, + bi[i] + μ + 3, + ci[i] + λ, + ) + + # D_th to D_gr + if !large && th.D_gr > cutoff + v += integrate_to_gamma( + D_th, + cutoff, + α_va * + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + v += integrate_to_gamma( + cutoff, + th.D_gr, + α_va * + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + else + v += integrate_to_gamma( + D_th, + th.D_gr, + α_va * + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + p3.β_va + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + end + + # D_gr to D_cr + if !large && th.D_cr > cutoff + v += integrate_to_gamma( + th.D_gr, + cutoff, + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ + 3, + ci[i] + λ, + ) + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + v += integrate_to_gamma( + cutoff, + th.D_cr, + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ + 3, + ci[i] + λ, + ) + else + v += integrate_to_gamma( + th.D_gr, + th.D_cr, + FT(π) / 6 * + th.ρ_g * + ai[i] * + N_0 * + (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ + 3, + ci[i] + λ, + ) + end + + # D_cr to Infinity + v_m_D_cr(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) + )^κ * + ai[i] * + D^(bi[i]) * + exp(-ci[i] * D) * + (α_va / (1 - F_r) * D^p3.β_va) * + N_0 * + D^μ * + exp(-λ * D) + if !large + (I, e) = QGK.quadgk(D -> v_m_D_cr(D), th.D_cr, cutoff) + v += I + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + (I, e) = QGK.quadgk(D -> v_m_D_cr(D), cutoff, Inf) + v += I + else + (I, e) = QGK.quadgk(D -> v_m_D_cr(D), th.D_cr, Inf) + v += I + end + end + end + + return v / q +end + +""" + terminal_velocity_number(p3, Chen2022, q, N, ρ_r, F_r) + + - p3 - a struct with P3 scheme parameters + - Chen 2022 - 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 number (total)-weighted fall speed + Eq C11 of Morrison and Milbrandt (2015) +""" +function terminal_velocity_number( + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + q::FT, + N::FT, + ρ_r::FT, + F_r::FT, + ρ_a::FT, +) where {FT} + + # Get the thresholds for different particles regimes + th = thresholds(p3, ρ_r, F_r) + D_th = D_th_helper(p3) + cutoff = FT(0.000625) # TODO add to the struct + + # Get the shape parameters + (λ, N_0) = distribution_parameter_solver(p3, q, N, ρ_r, F_r) + μ = DSD_μ(p3, λ) + + # TO DO: Change when each value used depending on type of particle + κ = FT(-1 / 6) #FT(1/3) + + # Redefine α_va to be in si units + α_va = α_va_si(p3) + + v = 0 + for i in 1:2 + if F_r == 0 + # Velocity coefficients for small particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + v += integrate_to_gamma( + FT(0), + D_th, + ai[i] * N_0, + bi[i] + μ, + ci[i] + λ, + ) + v += integrate_to_gamma( + D_th, + cutoff, + ai[i] * N_0 * (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + + # Get velocity coefficients for large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + v += integrate_to_gamma( + cutoff, + Inf, + ai[i] * N_0 * (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + else + # Velocity coefficients for small particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_small(Chen2022, ρ_a) + large = false + + v += integrate_to_gamma( + FT(0), + D_th, + ai[i] * N_0, + bi[i] + μ, + ci[i] + λ, + ) + + # D_th to D_gr + if !large && th.D_gr > cutoff + v += integrate_to_gamma( + D_th, + cutoff, + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + v += integrate_to_gamma( + cutoff, + th.D_gr, + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + else + v += integrate_to_gamma( + D_th, + th.D_gr, + ai[i] * + N_0 * + (16 * p3.ρ_i^2 * p3.γ^3 / (9 * FT(π) * α_va^2))^κ, + bi[i] + μ + κ * (3 * p3.σ - 2 * p3.β_va), + ci[i] + λ, + ) + end + + # D_gr to D_cr + if !large && th.D_cr > cutoff + v += integrate_to_gamma( + th.D_gr, + cutoff, + ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ, + ci[i] + λ, + ) + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + v += integrate_to_gamma( + cutoff, + th.D_cr, + ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ, + ci[i] + λ, + ) + else + v += integrate_to_gamma( + th.D_gr, + th.D_cr, + ai[i] * N_0 * (p3.ρ_i / th.ρ_g)^(2 * κ), + bi[i] + μ, + ci[i] + λ, + ) + end + + # D_cr to Infinity + v_n_D_cr(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) + )^κ * + ai[i] * + D^(bi[i]) * + exp(-ci[i] * D) * + N_0 * + D^μ * + exp(-λ * D) + if !large + (I, e) = QGK.quadgk(D -> v_n_D_cr(D), th.D_cr, cutoff) + v += I + + # Switch to large particles + (ai, bi, ci) = CO.Chen2022_vel_coeffs_large(Chen2022, ρ_a) + large = true + + (I, e) = QGK.quadgk(D -> v_n_D_cr(D), cutoff, Inf) + v += I + else + (I, e) = QGK.quadgk(D -> v_n_D_cr(D), th.D_cr, Inf) + v += I + end + end + end + + return v / N +end + +# TODO: add accurate number weighted velocity function that mimics the mass weighted one + +""" + 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 += integrate_to_gamma(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n += integrate_to_gamma(D_th, Inf, α_va * N_0, μ + p3.β_va + 1, λ) + else + n += integrate_to_gamma(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n += integrate_to_gamma(D_th, th.D_gr, α_va * N_0, μ + p3.β_va + 1, λ) + n += + integrate_to_gamma(th.D_gr, th.D_cr, π / 6 * th.ρ_g * N_0, μ + 4, λ) + n += integrate_to_gamma( + 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 4f5e4da24..3d60df363 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -19,6 +19,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7" MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845" MLJModels = "d491faf4-2d78-11e9-2867-c94bc002c0b7" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/test/p3_tests.jl b/test/p3_tests.jl index 14c38aec4..72953dfdb 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,102 @@ 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-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], + ] + 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_mass( + p3, + Chen2022.snow_ice, + q, + N, + ρ_r, + F_r, + ρ_a, + ) + + TT.@test calculated_vel > 0 + TT.@test paper_vals[i][j] ≈ calculated_vel atol = 3.14 + + end + end + end + + TT.@testset "Number-weighted terminal velocities" begin + 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_number( + p3, + Chen2022.snow_ice, + q, + N, + ρ_r, + F_r, + ρ_a, + ) + + TT.@test calculated_vel > 0 + TT.@test expected_vals[i][j] ≈ calculated_vel atol = 0.1 + + 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 cbe34b7e3..a8f8c43a1 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -86,6 +86,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) @@ -121,6 +122,28 @@ 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_mass, + (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air), + 2e6, + 4e4, + 2e3, + ) + bench_press( + P3.terminal_velocity_number, + (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(