diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 6bd92cc93..1350e5299 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -87,6 +87,17 @@ @article{Bigg1953 doi = {10.1088/0370-1301/66/8/309} } +@article{Blahak2010, + title = {Efficient approximation of the incomplete gamma function for use in cloud model applications}, + author = {Blahak, U.}, + journal = {Geoscientific Model Development}, + volume = {3}, + year = {2010}, + number = {2}, + pages = {329--336}, + doi = {10.5194/gmd-3-329-2010} +} + @article{BrownFrancis1995, title = {Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe}, author = {Philip R. A. Brown and Peter N. Francis}, diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index 42b9c4186..1f03e4b37 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -134,6 +134,14 @@ As a result ``L_{p3, tot}`` can be expressed as a sum of incomplete gamma functi where ``\Gamma \,(a, z) = \int_{z}^{\infty} \! t^{a - 1} e^{-t} \mathrm{d}D`` and ``\Gamma \,(a) = \Gamma \,(a, 0)`` for simplicity. +We rely on [Blahak2010](@cite) parameterization to compute the incomplete gamma functions. +Below plot roughly reproduce Figure 4 from [Blahak2010](@cite). +The relative error differences are mostly due to too small values returned by + our reference incomplete gamma function from Julia SpecialFunctions package. +```@example +include("plots/P3GammaAprox.jl") +``` +![](P3_GammaInc_error.svg) In our solver, we approximate ``\mu`` from ``L/N`` and keep it constant throughout the solving step. We approximate ``\mu`` by an exponential function given by the ``L/N`` points @@ -174,7 +182,7 @@ include("plots/P3LambdaErrorPlots.jl") ## Terminal Velocity -We use the [Chen2022](@cite) velocity parametrization: +We use the [Chen2022](@cite) velocity parameterization: ```math V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} ``` @@ -291,7 +299,7 @@ Visualizing mass-weighted terminal velocity as a function of ``F_{liq}``, ``F_{r medium, and large particles, we have mostly graupel (``D_{gr} < D < D_{cr}``) for small ``D_m`` and mostly partially rimed ice (``D > D_{cr}``) for medium and large ``D_m``. Thus, we can attribute the non-monotonic behavior of velocity with ``F_liq`` in the medium and large ``D_m`` plots to the variations in ``\phi`` caused by nonspherical particle shape, whereas the small ``D_m`` plot - confirms a mmore monotonic change in ``V_m`` for spherical ice. The ``L`` and ``N`` values used to generate small, medium, and large ``D_{m}`` + confirms a more monotonic change in ``V_m`` for spherical ice. The ``L`` and ``N`` values used to generate small, medium, and large ``D_{m}`` are the same as in the plot above. ```@example @@ -301,7 +309,7 @@ include("plots/P3TerminalVelocity_F_liq_rim.jl") When modifying process rates, we now need to consider whether they are concerned with the ice core or the whole particle, in addition to whether they become sources - and sinks of different prognostic variables with the inclusion of ``F_{liq}``. + and sinks of different prognostic variables with the inclusion of ``F_{liq}``. With the addition of liquid fraction, too, come new process rates. diff --git a/docs/src/plots/P3GammaAprox.jl b/docs/src/plots/P3GammaAprox.jl new file mode 100644 index 000000000..3779776e0 --- /dev/null +++ b/docs/src/plots/P3GammaAprox.jl @@ -0,0 +1,145 @@ +import CairoMakie as PL +PL.activate!(type = "svg") + +import SpecialFunctions as SF + +import Thermodynamics as TD +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.P3Scheme as P3 + +#! format: off + +FT = Float64 +res = 100 # plot resolution + +# Plot Fig 3 from Blahak 2010 +# --------------------------- + +a_range = range(1e-1, stop = 30, length = res) +c1 = [P3.c₁(a) for a in a_range] +c2 = [P3.c₂(a) for a in a_range] +c3 = [P3.c₃(a) for a in a_range] +c4 = [P3.c₄(a) for a in a_range] + +fig = PL.Figure(size = (1500, 1000), fontsize=22, linewidth=3) + +ax1 = PL.Axis(fig[1, 1]) +ax2 = PL.Axis(fig[1, 2]) +ax3 = PL.Axis(fig[2, 1]) +ax4 = PL.Axis(fig[2, 2]) + +ax1.ylabel = "c₁" +ax2.ylabel = "c₂" +ax3.ylabel = "c₃" +ax4.ylabel = "c₄" + +ax3.xlabel = "a" +ax4.xlabel = "a" + +PL.ylims!(ax2, 0, 2) +PL.ylims!(ax4, 0, 10) + +c1l = PL.lines!(ax1, a_range, c1) +c2l = PL.lines!(ax2, a_range, c2) +c3l = PL.lines!(ax3, a_range, c3) +c4l = PL.lines!(ax4, a_range, c4) + +PL.save("P3_GammaInc_c_coeffs.svg", fig) + + +# Plot figure 1 from Blahak 2010 +# ------------------------------ + +fig2 = PL.Figure(size = (1000, 500), fontsize=22) +ax = PL.Axis(fig2[1,1]) +ax.xlabel = "z" +ax.ylabel = "P" + +z_range = range(1e-1, stop = 25, length=res) + +P_a1 = [P3.Γ_lower(FT(1), z) / SF.gamma(FT(1)) for z in z_range] +P_a3 = [P3.Γ_lower(FT(3), z) / SF.gamma(FT(3)) for z in z_range] +P_a5 = [P3.Γ_lower(FT(5), z) / SF.gamma(FT(5)) for z in z_range] +P_a10 = [P3.Γ_lower(FT(10), z) / SF.gamma(FT(10)) for z in z_range] +P_a15 = [P3.Γ_lower(FT(15), z) / SF.gamma(FT(15)) for z in z_range] +P_a30 = [P3.Γ_lower(FT(30), z) / SF.gamma(FT(30)) for z in z_range] + +G_a1 = [SF.gamma_inc(FT(1), z)[1] for z in z_range] +G_a3 = [SF.gamma_inc(FT(3), z)[1] for z in z_range] +G_a5 = [SF.gamma_inc(FT(5), z)[1] for z in z_range] +G_a10 = [SF.gamma_inc(FT(10), z)[1] for z in z_range] +G_a15 = [SF.gamma_inc(FT(15), z)[1] for z in z_range] +G_a30 = [SF.gamma_inc(FT(30), z)[1] for z in z_range] + +l1 = PL.lines!(ax, z_range, P_a1, label="a=1", color="midnightblue", linewidth=3) +l3 = PL.lines!(ax, z_range, P_a3, label="a=3", color="slateblue", linewidth=3) +l5 = PL.lines!(ax, z_range, P_a5, label="a=5", color="steelblue1", linewidth=3) +l10 = PL.lines!(ax, z_range, P_a10, label="a=10", color="orange", linewidth=3) +l15 = PL.lines!(ax, z_range, P_a15, label="a=15", color="orangered2", linewidth=3) +l30 = PL.lines!(ax, z_range, P_a30, label="a=30", color="brown", linewidth=3) + +g1 = PL.lines!(ax, z_range, G_a1, label="a=1", color="midnightblue", linewidth=3, linestyle=:dash) +g3 = PL.lines!(ax, z_range, G_a3, label="a=1", color="slateblue", linewidth=3, linestyle=:dash) +g5 = PL.lines!(ax, z_range, G_a5, label="a=1", color="steelblue1", linewidth=3, linestyle=:dash) +g10 = PL.lines!(ax, z_range, G_a10, label="a=1", color="orange", linewidth=3, linestyle=:dash) +g15 = PL.lines!(ax, z_range, G_a15, label="a=1", color="orangered2", linewidth=3, linestyle=:dash) +g30 = PL.lines!(ax, z_range, G_a30, label="a=1", color="brown", linewidth=3, linestyle=:dash) + +PL.Legend( + fig2[1, 2], + [l1, l3, l5, l10, l15, l30], + ["a=1", "a=3", "a=5", "a=10", "a=15", "a=30"], +) +PL.save("P3_GammaInc_P_check.svg", fig2) + +# Plot figure 4 from Blahak 2010 +# ------------------------------ +# The relative error blows up for cases where the Γ_lower is very small. +# The reference value we obtain from Julia Special Functions is off +# by several orders of magnitude when compared with Wolfram Mathematica + +y_range = range(1e-1, stop = 2.5, length = res) +z_range = y_range .* (a_range .+ 1) + +val = zeros(res, res) +err_rel = zeros(res, res) +err_abs = zeros(res, res) +for it in range(1, stop=res, step=1) # a + for jt in range(1, stop=res, step=1) # z + a = FT(a_range[it]) + z = FT(z_range[jt]) + val[it, jt] = P3.Γ_lower(a, z) / SF.gamma(a) + err_abs[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) - SF.gamma_inc(a, z)[1] + err_rel[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) / SF.gamma_inc(a, z)[1] - 1 + end +end + +fig3 = PL.Figure(size = (1800, 600), fontsize=22) +ax1 = PL.Axis(fig3[1, 1], xlabel = "a", ylabel = "z / (a + 1)", title = "Relative error") +ax2 = PL.Axis(fig3[1, 3], xlabel = "a", title = "Absolute error") +ax3 = PL.Axis(fig3[1, 5], xlabel = "a", title = "Γ_lower(a, z) / Γ(a)") + +PL.ylims!(ax1, 0, 2.5) +PL.ylims!(ax2, 0, 2.5) +PL.ylims!(ax3, 0, 2.5) + +cplot1 = PL.heatmap!( + ax1, a_range, y_range, err_rel, + colormap = PL.cgrad(:viridis, 20, categorical=true), + colorrange = (-0.1, 10), highclip=:red, lowclip=:orange, +) +PL.Colorbar(fig3[1,2], cplot1) + +cplot2 = PL.contourf!( + ax2, a_range, y_range, err_abs, + levels=-0.03:0.005:0.03, + colormap="RdBu" +) +PL.Colorbar(fig3[1,4], cplot2) + +cplot3 = PL.contourf!( + ax3, a_range, y_range, val, +) +PL.Colorbar(fig3[1,6], cplot3) + +PL.save("P3_GammaInc_error.svg", fig3) diff --git a/src/P3_helpers.jl b/src/P3_helpers.jl index 38c9c40b0..f5a560d38 100644 --- a/src/P3_helpers.jl +++ b/src/P3_helpers.jl @@ -1,8 +1,99 @@ -# Some wrappers to cast types from SF.gamma +# Wrapper 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)) +# TODO - replace with parameterization of our own Γ(a::FT) where {FT <: Real} = FT(SF.gamma(a)) +# helper functions for Γ_approx +function c₁(a::FT) where {FT <: Real} + p1 = FT(9.4368392235e-03) + p2 = -FT(1.0782666481e-04) + p3 = -FT(5.8969657295e-06) + p4 = FT(2.8939523781e-07) + p5 = FT(1.0043326298e-01) + p6 = FT(5.5637848465e-01) + p = [p1, p2, p3, p4, p5, p6] + + ret = 1 + p[5] * (exp(-p[6] * a) - 1) + for it in range(1, 4, step = 1) + ret += p[it] * a^it + end + return ret +end +function c₂(a::FT) where {FT <: Real} + q1 = FT(1.1464706419e-01) + q2 = FT(2.6963429121) + q3 = -FT(2.9647038257) + q4 = FT(2.1080724954) + q = [q1, q2, q3, q4] + + ret = q[1] + for it in range(2, 4, step = 1) + ret += q[it] / a^(it - 1) + end + return ret +end +function c₃(a::FT) where {FT <: Real} + r1 = FT(0) + r2 = FT(1.1428716184) + r3 = -FT(6.6981186438e-03) + r4 = FT(1.0480765092e-04) + r = [r1, r2, r3, r4] + + ret = 0 + for it in range(1, 4, step = 1) + ret += r[it] * a^(it - 1) + end + return ret +end +function c₄(a::FT) where {FT <: Real} + s1 = FT(1.0356711153) + s2 = FT(2.3423452308) + s3 = -FT(3.6174503174e-01) + s4 = -FT(3.1376557650) + s5 = FT(2.9092306039) + s = [s1, s2, s3, s4, s5] + + ret = 0 + for it in range(1, 5, step = 1) + ret += s[it] / a^(it - 1) + end + return ret +end + +""" + Γ_lower(a, z) + +An approximated lower incomplete Gamma function based on Blahak 2010: +doi:10.5194/gmd-3-329-2010 +https://gmd.copernicus.org/articles/3/329/2010/gmd-3-329-2010.pdf +""" +function Γ_lower(a::FT, z::FT) where {FT <: Real} + + if isnan(z) || z <= FT(0) + return FT(0) + else + W(a, z) = FT(0.5) + FT(0.5) * tanh(c₂(a) * (z - c₃(a))) + + return exp(-z) * + z^a * + ( + 1 / a + + c₁(a) * z / a / (a + 1) + + (c₁(a) * z)^2 / a / (a + 1) / (a + 2) + ) * + (1 - W(a, z)) + Γ(a) * W(a, z) * (1 - c₄(a)^(-z)) + end +end + +""" + Γ_upper(a, z) + +An approximated upper incomplete Gamma function computed as Γ(a) - Γ_lower(a, z) +""" +function Γ_upper(a::FT, z::FT) where {FT <: Real} + return Γ(a) - Γ_lower(a, z) +end + """ ∫_Γ(x₀, x_end, c1, c2, c3) @@ -16,11 +107,13 @@ 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 == FT(Inf) - return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3)) + return c1 * c3^(-c2 - 1) * (Γ_upper(1 + c2, x₀ * c3)) elseif x₀ == 0 - return c1 * c3^(-c2 - 1) * (Γ(1 + c2) - Γ(1 + c2, x_end * c3)) + return c1 * c3^(-c2 - 1) * (Γ_lower(1 + c2, x_end * c3)) else - return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3) - Γ(1 + c2, x_end * c3)) + return c1 * + c3^(-c2 - 1) * + (Γ_upper(1 + c2, x₀ * c3) - Γ_upper(1 + c2, x_end * c3)) end end diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl index f0d046f99..2cf61001c 100644 --- a/src/P3_particle_properties.jl +++ b/src/P3_particle_properties.jl @@ -312,12 +312,12 @@ end - F_rim - rime mass fraction (L_rim/ L_ice) [-] - th - P3 particle properties thresholds -Returns the aspect ratio (ϕ) for an ice particle with mass, cross-sectional area, -and ice density determined using the size-dependent particle property regimes -following Morrison and Milbrandt (2015). The density of nonspherical -particles is assumed to be equal to the particle mass divided by the volume of a -spherical particle with the same D_max. -Assuming zero liquid fraction. +Returns the aspect ratio (ϕ) for an ice particle with mass, cross-sectional area, +and ice density determined using the size-dependent particle property regimes +following Morrison and Milbrandt (2015). The density of nonspherical +particles is assumed to be equal to the particle mass divided by the volume of a +spherical particle with the same D_max. +Assuming zero liquid fraction. """ function ϕᵢ(p3::PSP3, D::FT, F_rim::FT, th) where {FT} F_liq = FT(0) diff --git a/src/P3_processes.jl b/src/P3_processes.jl index f2d292cb5..07e528770 100644 --- a/src/P3_processes.jl +++ b/src/P3_processes.jl @@ -90,7 +90,7 @@ function ice_melt( distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_) N(D) = N′ice(p3, D, λ, N_0) # ... and D_max for the integral - bound = get_ice_bound(p3, λ, FT(1e-6)) + bound = get_ice_bound(p3, λ, N_ice, FT(1e-6)) # Ice particle terminal velocity v(D) = ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th) diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index caa7e1f0e..848a47c5a 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -50,49 +50,49 @@ end 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) + #@assert D > FT(0) + return ifelse(D > 0, N_0 * D^DSD_μ(p3, λ) * exp(-λ * D), FT(0)) end -# """ -# get_ice_bound(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::PSP3, λ::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 """ - get_ice_bound(p3, λ, tolerance) + get_ice_bound(p3, λ, N, rtol) - p3 - a struct containing p3 parameters - λ - shape parameters of ice distribution - - tolerance - tolerance for how much of the distribution we want to integrate over + - N - ice number concentration + - rtol - relative tolerance for root solver - Ice size distribution upper bound. Set to a fixed number as a temporary fix. - The commented above function is not stable right now. + Ice size distribution upper bound for numerical integrals, such that + total number concentration is equal to the integral over the size distribution + with relative tolerance rtol. """ -function get_ice_bound(p3::PSP3, λ::FT, tolerance::FT) where {FT} - return FT(1e-2) +function get_ice_bound(p3::PSP3, λ::FT, N::FT, rtol::FT) where {FT} + #return FT(1e-2) - From what I'm seeing so far, this is not such a bad guess + #We might want to reconsider using it in the future (after more testing). + μ = DSD_μ(p3, λ) + N₀ = DSD_N₀(p3, N, λ) + + ice_bound_problem(D_max_hat) = + N - N₀ / λ^(μ + 1) * Γ_lower(μ + 1, D_max_hat) + + D_guess_low = FT(1e-6) + D_guess_high = FT(1e-1) + + if ice_bound_problem(D_guess_low * λ) <= rtol + return D_guess_low + elseif ice_bound_problem(D_guess_high * λ) <= rtol + return D_guess_high + else + bound = + RS.find_zero( + ice_bound_problem, + RS.SecantMethod(D_guess_low * λ, D_guess_high * λ), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(rtol), + 5, + ).root + return bound / λ + end end """ @@ -160,7 +160,6 @@ function L_over_N_gamma( D_th = D_th_helper(p3) λ = exp(log_λ) N = Γ(1 + μ) / (λ^(1 + μ)) - return ifelse( F_rim == FT(0), (p3_F_liq_average( diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index 0598eb722..c20569a22 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -181,7 +181,7 @@ function ice_terminal_velocity( # get the size distribution parameters (λ, N₀) = distribution_parameter_solver(p3, L, N, ρ_r, F_rim, F_liq) # get the integral limit - D_max = get_ice_bound(p3, λ, eps(FT)) + D_max = get_ice_bound(p3, λ, N, 1e-8) # ∫N(D) m(D) v(D) dD v_m = QGK.quadgk( diff --git a/test/p3_tests.jl b/test/p3_tests.jl index 16ba42c2e..8d9b02538 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -8,13 +8,35 @@ import Thermodynamics as TD import QuadGK as QGK -@info "P3 Scheme Tests" +@info("P3 Scheme Tests") -function test_p3_thresholds(FT) +function test_incomplete_gamma_function_approximation(FT) + # A smoke tests to see if the reference values are not changed + TT.@testset "Incomplete Γ function approximation" begin + + Γ_lower_reference = [ + [0.6171173028812539, 0.9999647375921944, 0.9999999999999561], + [66.43238525128415, 192195.93434392574, 362553.66f0], + [9.801186821098075e11, 8.309791f14, 1.1845235f17], + ] + + a = [FT(1), FT(10), FT(20)] + z = [FT(1), FT(10), FT(30)] + + for it in [1, 2, 3] + for jt in [1, 2, 3] + TT.@test P3.Γ_lower(a[it], z[jt]) ≈ Γ_lower_reference[it][jt] rtol = + 1e-3 + end + end + end +end + +function test_thresholds_solver(FT) p3 = CMP.ParametersP3(FT) - TT.@testset "thresholds (nonlinear solver function)" begin + TT.@testset "Thresholds - nonlinear solver" begin # initialize test values: ρ_r = FT(400) @@ -89,7 +111,7 @@ function test_p3_thresholds(FT) end end - TT.@testset "mass, area, density and aspect ratio tests" begin + TT.@testset "Thresholds - mass, area, density, aspect ratio" begin # values ρ_r = FT(500) F_rim = FT(0.5) @@ -186,11 +208,11 @@ function test_p3_thresholds(FT) end end -function test_p3_shape_solver(FT) +function test_shape_solver(FT) p3 = CMP.ParametersP3(FT) - TT.@testset "shape parameters (nonlinear solver function)" begin + TT.@testset "Shape parameters - nonlinear solver" begin # initialize test values: ep = 1 #1e4 * eps(FT) @@ -248,7 +270,7 @@ function test_particle_terminal_velocities(FT) Chen2022 = CMP.Chen2022VelType(FT) ρ_a = FT(1.2) - TT.@testset "Chen 2022 - Rain" begin + TT.@testset "Smoke tests for rain particle terminal vel from Chen 2022" begin Ds = range(FT(1e-6), stop = FT(1e-5), length = 5) expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144] for i in axes(Ds, 1) @@ -258,7 +280,7 @@ function test_particle_terminal_velocities(FT) end end - TT.@testset "Chen 2022 - Ice" begin + TT.@testset "Smoke tests for ice particle terminal vel from Chen 2022" begin F_rim = FT(0.5) ρ_r = FT(500) F_liq = FT(0) @@ -302,7 +324,7 @@ function test_particle_terminal_velocities(FT) end end - TT.@testset "Chen 2022 - Mixed-Phase" begin + TT.@testset "Smoke tests for mixed phase particle terminal velocity" begin F_rim = FT(0.5) F_liq = FT(0.5) ρ_r = FT(500) @@ -354,202 +376,136 @@ function test_bulk_terminal_velocities(FT) L = FT(0.22) N = FT(1e6) ρ_a = FT(1.2) - ρ_rs = [FT(200), FT(400), FT(600), FT(800)] - F_rims = [FT(0), FT(0.2), FT(0.4), FT(0.6), FT(0.8)] - F_liqs = [FT(0), FT(0.5), FT(1)] + ρ_r = FT(800) + F_rims = [FT(0), FT(0.6)] + F_liqs = [FT(0.5), FT(1)] TT.@testset "Mass and number weighted terminal velocities" begin - # commented out reference values: - # these are computed with the old get_ice_bound() - # reference_vals_m = [ - # [ - # [7.79, 7.27, 6.66, 5.94, 5.25], - # [7.79, 7.26, 6.62, 5.83, 4.82], - # [7.79, 7.26, 6.62, 5.81, 4.70], - # [7.79, 7.25, 6.61, 5.80, 4.65], - # ], - # [ - # [5.19, 5.17, 5.15, 5.12, 5.09], - # [5.19, 5.17, 5.13, 5.08, 4.97], - # [5.19, 5.17, 5.13, 5.07, 4.92], - # [5.19, 5.16, 5.13, 5.06, 4.89], - # ], - # [ - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # ], - # ] - - # reference_vals_n = [ - # [ - # [3.65, 3.37, 3.05, 2.64, 2.14], - # [3.65, 3.37, 3.04, 2.62, 2.04], - # [3.65, 3.37, 3.04, 2.62, 2.02], - # [3.65, 3.37, 3.04, 2.62, 2.01], - # ], - # [ - # [1.69, 1.68, 1.68, 1.66, 1.64], - # [1.69, 1.68, 1.68, 1.66, 1.62], - # [1.69, 1.68, 1.67, 1.66, 1.61], - # [1.69, 1.68, 1.67, 1.66, 1.61], - # ], - # [ - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # ], - # ] - - # reference_vals_m_ϕ = [ - # [ - # [4.23, 4.65, 4.9, 4.94, 4.89], - # [4.23, 4.65, 4.88, 4.84, 4.44], - # [4.23, 4.65, 4.87, 4.82, 4.33], - # [4.23, 4.65, 4.87, 4.81, 4.28], - # ], - # [ - # [4.36, 4.51, 4.67, 4.82, 4.99], - # [4.36, 4.51, 4.64, 4.76, 4.82], - # [4.36, 4.50, 4.64, 4.74, 4.77], - # [4.36, 4.50, 4.64, 4.74, 4.74], - # ], - # [ - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # [5.42, 5.42, 5.42, 5.42, 5.42], - # ], - # ] - - # reference_vals_n_ϕ = [ - # [ - # [2.21, 2.34, 2.39, 2.31, 2.04], - # [2.21, 2.33, 2.37, 2.27, 1.94], - # [2.21, 2.33, 2.36, 2.25, 1.91], - # [2.21, 2.33, 2.36, 2.24, 1.90], - # ], - # [ - # [1.47, 1.52, 1.56, 1.59, 1.60], - # [1.47, 1.52, 1.56, 1.59, 1.60], - # [1.47, 1.51, 1.55, 1.58, 1.59], - # [1.47, 1.51, 1.55, 1.58, 1.58], - # ], - # [ - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # [1.62, 1.62, 1.62, 1.62, 1.62], - # ], - # ] - #! format: off - reference_vals_m = [6.63 6.63 6.63 6.63; 5.19 5.19 5.19 5.19; 5.42 5.42 5.42 5.42;;; 6.56 6.55 6.55 6.55; 5.17 5.17 5.17 5.16; 5.42 5.42 5.42 5.42;;; 6.32 6.29 6.28 6.28; 5.15 5.13 5.13 5.13; 5.42 5.42 5.42 5.42;;; 5.84 5.74 5.72 5.71; 5.12 5.08 5.07 5.06; 5.42 5.42 5.42 5.42;;; 5.24 4.81 4.69 4.65; 5.09 4.97 4.92 4.89; 5.42 5.42 5.42 5.42] - + # Liquid fraction = 0 - reference_vals_n = [3.59 3.59 3.59 3.59; 1.69 1.69 1.69 1.69; 1.62 1.62 1.62 1.62;;; 3.34 3.34 3.34 3.34; 1.68 1.68 1.68 1.68; 1.62 1.62 1.62 1.62;;; 3.03 3.03 3.03 3.03; 1.68 1.68 1.67 1.67; 1.62 1.62 1.62 1.62;;; 2.64 2.62 2.62 2.62; 1.66 1.66 1.66 1.66; 1.62 1.62 1.62 1.62;;; 2.14 2.04 2.02 2.01; 1.64 1.62 1.61 1.61; 1.62 1.62 1.62 1.62] + ref_v_n = [3.637320375996216, 2.6130243256593313] + ref_v_n_ϕ = [2.2058692424028923, 2.24224820764301] + ref_v_m = [7.7208191246659705, 5.751352472270078] + ref_v_m_ϕ = [4.193799087495265, 4.773592179767702] - reference_vals_m_ϕ = [3.66 3.66 3.66 3.66; 4.36 4.36 4.36 4.36; 5.42 5.42 5.42 5.42;;; 4.23 4.23 4.23 4.23; 4.51 4.51 4.5 4.5; 5.42 5.42 5.42 5.42;;; 4.67 4.64 4.64 4.64; 4.67 4.64 4.64 4.64; 5.42 5.42 5.42 5.42;;; 4.86 4.77 4.75 4.74; 4.82 4.76 4.74 4.74; 5.42 5.42 5.42 5.42;;; 4.88 4.44 4.32 4.28; 4.99 4.82 4.77 4.74; 5.42 5.42 5.42 5.42] + for k in 1:length(F_rims) + F_rim = F_rims[k] + F_liq = FT(0) + vel = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + false, + ) + vel_ϕ = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + true, + ) - reference_vals_n_ϕ = [2.18 2.18 2.18 2.18; 1.47 1.47 1.47 1.47; 1.62 1.62 1.62 1.62;;; 2.32 2.32 2.31 2.31; 1.52 1.52 1.51 1.51; 1.62 1.62 1.62 1.62;;; 2.38 2.36 2.35 2.35; 1.56 1.56 1.55 1.55; 1.62 1.62 1.62 1.62;;; 2.3 2.26 2.25 2.25; 1.59 1.59 1.58 1.58; 1.62 1.62 1.62 1.62;;; 2.04 1.94 1.91 1.9; 1.6 1.6 1.59 1.58; 1.62 1.62 1.62 1.62] - #! format: on + # number weighted + TT.@test vel[1] > 0 + TT.@test vel_ϕ[1] > 0 + TT.@test vel[1] ≈ ref_v_n[k] rtol = 1e-6 + TT.@test vel_ϕ[1] ≈ ref_v_n_ϕ[k] rtol = 1e-6 + + # mass weighted + TT.@test vel[2] > 0 + TT.@test vel_ϕ[2] > 0 + TT.@test vel[2] ≈ ref_v_m[k] rtol = 1e-6 + TT.@test vel_ϕ[2] ≈ ref_v_m_ϕ[k] rtol = 1e-6 + + # slower with aspect ratio + TT.@test vel_ϕ[1] <= vel[1] + TT.@test vel_ϕ[2] <= vel[2] + end - for j in 1:length(ρ_rs) - for k in 1:length(F_rims) - for i in 1:length(F_liqs) - ρ_r = ρ_rs[j] - F_rim = F_rims[k] - F_liq = F_liqs[i] - - calculated_vel = P3.ice_terminal_velocity( - p3, - Chen2022, - L, - N, - ρ_r, - F_rim, - F_liq, - ρ_a, - false, - ) - calculated_vel_ϕ = P3.ice_terminal_velocity( - p3, - Chen2022, - L, - N, - ρ_r, - F_rim, - F_liq, - ρ_a, - true, - ) - - # number weighted - TT.@test calculated_vel[1] > 0 - TT.@test reference_vals_n[i, j, k] ≈ calculated_vel[1] atol = - 0.1 - TT.@test calculated_vel_ϕ[1] > 0 - TT.@test reference_vals_n_ϕ[i, j, k] ≈ calculated_vel_ϕ[1] atol = - 0.1 - - # mass weighted - TT.@test calculated_vel[2] > 0 - TT.@test reference_vals_m[i, j, k] ≈ calculated_vel[2] atol = - 0.1 - TT.@test calculated_vel_ϕ[2] > 0 - TT.@test reference_vals_m_ϕ[i, j, k] ≈ calculated_vel_ϕ[2] atol = - 0.1 - - TT.@test calculated_vel_ϕ[1] <= calculated_vel[1] - TT.@test calculated_vel_ϕ[2] <= calculated_vel[2] - end - end + # Liquid fraction != 0 + ref_v_n = [1.674591925057434, 1.6180970319460353] + ref_v_n_ϕ = [1.549777478756061, 1.6180970319460353] + ref_v_m = [5.126648558302173, 5.416679316254198] + ref_v_m_ϕ = [4.6358422594886495, 5.416679316254198] + + for k in 1:length(F_liqs) + F_liq = F_liqs[k] + F_rim = FT(0.4) + vel = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + false, + ) + vel_ϕ = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + true, + ) + # number weighted + TT.@test vel[1] > 0 + TT.@test vel_ϕ[1] > 0 + TT.@test vel[1] ≈ ref_v_n[k] rtol = 1e-6 + TT.@test vel_ϕ[1] ≈ ref_v_n_ϕ[k] rtol = 1e-6 + + # mass weighted + TT.@test vel[2] > 0 + TT.@test vel_ϕ[2] > 0 + TT.@test vel[2] ≈ ref_v_m[k] rtol = 1e-6 + TT.@test vel_ϕ[2] ≈ ref_v_m_ϕ[k] rtol = 1e-6 + + # slower with aspect ratio + TT.@test vel_ϕ[1] <= vel[1] + TT.@test vel_ϕ[2] <= vel[2] end end - TT.@testset "Mass-weighted mean diameters" begin - F_liq = FT(0) # test against paper - 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_rims) - ρ_r = ρ_rs[i] - F_rim = F_rims[j] - - calculated_dm = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) * 1e3 - - TT.@test calculated_dm > 0 - TT.@test paper_vals[i][j] ≈ calculated_dm atol = 3.14 - - end + F_liq = FT(0) + ref_vals = [0.00536934536778844, 0.003321590762128599] + for j in 1:length(F_rims) + Dₘ = P3.D_m(p3, L, N, ρ_r, F_rims[j], F_liq) + TT.@test Dₘ > 0 + TT.@test Dₘ ≈ ref_vals[j] end # nonzero F_liq - ρ_r = ρ_rs[4] - F_rim = F_rims[4] - F_liqs = [FT(0), FT(0.33), FT(0.67), FT(1)] - expected_vals = - [FT(0.00333689), FT(0.000892223), FT(0.00124423), FT(0.00164873)] + F_rim = F_rims[2] + F_liqs = [FT(0.33), FT(1)] + ref_vals = [FT(0.0021371920600012184), FT(0.0016487352655895715)] for i in eachindex(F_liqs) - calculated_dm = P3.D_m(p3, L, N, ρ_r, F_rim, F_liqs[i]) - expected_vals[i] = calculated_dm - TT.@test expected_vals[i] ≈ calculated_dm atol = 1e-6 + Dₘ = P3.D_m(p3, L, N, ρ_r, F_rim, F_liqs[i]) + TT.@test Dₘ ≈ ref_vals[i] end end end -function test_integrals(FT) +function test_numerical_integrals(FT) p3 = CMP.ParametersP3(FT) Chen2022 = CMP.Chen2022VelType(FT) N = FT(1e8) - Ls = range(0.001, stop = 0.005, length = 5) + Ls = range(FT(0.001), stop = FT(0.005), length = 5) ρ_r = FT(500) F_rims = [FT(0), FT(0.5)] ρ_a = FT(1.2) @@ -557,35 +513,42 @@ function test_integrals(FT) tolerance = eps(FT) F_liq = FT(0) # test F_liq = 0 - TT.@testset "Gamma vs Integral Comparison" begin + TT.@testset "Numerical integrals sanity checks for N, velocity and diameter" begin for F_rim in F_rims for i in axes(Ls, 1) L = Ls[i] - # Velocity comparisons - vel_N, vel_m = P3.ice_terminal_velocity( + # Get shape parameters, thresholds and intergal bounds + λ, N_0 = P3.distribution_parameter_solver( p3, - Chen2022, L, N, ρ_r, F_rim, F_liq, - ρ_a, - use_aspect_ratio, ) + th = P3.thresholds(p3, ρ_r, F_rim) + ice_bound = P3.get_ice_bound(p3, λ, N, tolerance) - λ, N_0 = P3.distribution_parameter_solver( + # Number concentration comparison + f(d) = P3.N′ice(p3, d, λ, N_0) + qgk_N, = QGK.quadgk(d -> f(d), FT(0), ice_bound) + TT.@test N ≈ qgk_N rtol = sqrt(eps(FT)) + + # Bulk velocity comparison + vel_N, vel_m = P3.ice_terminal_velocity( p3, + Chen2022, L, N, ρ_r, F_rim, F_liq, + ρ_a, + use_aspect_ratio, ) - th = P3.thresholds(p3, ρ_r, F_rim) - ice_bound = P3.get_ice_bound(p3, λ, tolerance) - vel(d) = P3.ice_particle_terminal_velocity( + + particle_vel(d) = P3.ice_particle_terminal_velocity( p3, d, Chen2022.snow_ice, @@ -594,27 +557,25 @@ function test_integrals(FT) th, use_aspect_ratio, ) - f(d) = vel(d) * P3.N′ice(p3, d, λ, N_0) + g(d) = particle_vel(d) * P3.N′ice(p3, d, λ, N_0) - qgk_vel_N, = QGK.quadgk(d -> f(d) / N, FT(0), ice_bound) + qgk_vel_N, = QGK.quadgk(d -> g(d) / N, FT(0), ice_bound) qgk_vel_m, = QGK.quadgk( - d -> f(d) * P3.p3_mass(p3, d, F_rim, F_liq, th) / L, + d -> g(d) * P3.p3_mass(p3, d, F_rim, F_liq, th) / L, FT(0), ice_bound, ) - - TT.@test vel_N ≈ qgk_vel_N rtol = 1e-5 - TT.@test vel_m ≈ qgk_vel_m rtol = 1e-5 + TT.@test vel_N ≈ qgk_vel_N rtol = 1e-6 + TT.@test vel_m ≈ qgk_vel_m rtol = 1e-6 # Dₘ comparisons D_m = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) - f_d(d) = + h(d) = d * P3.p3_mass(p3, d, F_rim, F_liq, th) * P3.N′ice(p3, d, λ, N_0) - qgk_D_m, = QGK.quadgk(d -> f_d(d) / L, FT(0), ice_bound) - - TT.@test D_m ≈ qgk_D_m rtol = 1e-8 + qgk_D_m, = QGK.quadgk(d -> h(d) / L, FT(0), ice_bound) + TT.@test D_m ≈ qgk_D_m rtol = 1e-2 end end end @@ -771,7 +732,7 @@ function test_p3_het_freezing(FT) expected_freeze_L = [1.544e-22, 1.068e-6, 0.0001428, 0.0001428, 0.0001428, 0.0001428] expected_freeze_N = [1.082e-10, 747647.5, N_liq, N_liq, N_liq, N_liq] - qᵥ_range = range(0.5e-3, stop = 1.5e-3, length = 6) + qᵥ_range = range(FT(0.5e-3), stop = FT(1.5e-3), length = 6) for it in range(1, 6) qₚ = TD.PhasePartition(FT(qᵥ_range[it]), FT(2e-4), FT(0)) @@ -845,8 +806,8 @@ function test_p3_melting(FT) TT.@test rate.dNdt >= 0 TT.@test rate.dLdt >= 0 - TT.@test rate.dNdt ≈ 171236.02291689217 rtol = 1e-6 - TT.@test rate.dLdt ≈ 8.561801145844608e-5 rtol = 1e-6 + TT.@test rate.dNdt ≈ 171326.17992979713 rtol = 1e-6 + TT.@test rate.dLdt ≈ 8.566308996489856e-5 rtol = 1e-6 T_vwarm = FT(273.15 + 0.1) rate = P3.ice_melt( @@ -868,22 +829,31 @@ function test_p3_melting(FT) end end +for FT in [Float32, Float64] + @info("Testing " * string(FT)) + + # numerics + test_incomplete_gamma_function_approximation(FT) + test_thresholds_solver(FT) + + # velocity + test_particle_terminal_velocities(FT) + + # processes + test_p3_het_freezing(FT) +end -println("Testing Float32") -test_p3_thresholds(Float32) -test_particle_terminal_velocities(Float32) -# TODO - fix instability in get_ice_bound for 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) -# test_bulk_terminal_velocities(Float32) - -println("Testing Float64") -test_p3_thresholds(Float64) -test_p3_shape_solver(Float64) -test_particle_terminal_velocities(Float64) -test_bulk_terminal_velocities(Float64) -test_integrals(Float64) -test_p3_het_freezing(Float64) -test_p3_melting(Float64) -#test_tendencies(Float64) +# TODO - make work for Float32 +for FT in [Float64] + @info("Tests that only work for Float64") + + # numerics + test_shape_solver(FT) + test_numerical_integrals(FT) + + # velocity + test_bulk_terminal_velocities(FT) + + # processes + test_p3_melting(FT) +end