diff --git a/Project.toml b/Project.toml index a3a887d5ef..cc4f4f1b1b 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 8e18b380db..0c098fb92b 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -35,6 +35,17 @@ @article{Beheng1994 doi = {10.1016/0169-8095(94)90020-5} } +@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}, + journal = {Journal of Atmospheric and Oceanic Technology}, + volume = {12}, + number = {2}, + pages = {410--414}, + year = {1995}, + doi = {10.1175/1520-0426(1995)012<0410:IMOTIW>2.0.CO;2} +} + @article{Chen2022, title={Accurate parameterization of precipitation particles' fall speeds for bulk cloud microphysics schemes}, author={Chen, Jen-Ping and Hsieh, Ting-Wei and Lin, Chiu-Yi and Yu, Cheng-Ku}, @@ -114,6 +125,17 @@ @article{Harrington1995 year = {1995} } +@article{Heymsfield2003, + title = {Properties of Tropical and Midlatitude Ice Cloud Particle Ensembles. Part II: Applications for Mesoscale and Climate Models}, + author = {Andrew J. Heymsfield}, + journal = {Journal of the Atmospheric Sciences}, + volume = {60}, + number = {21}, + doi = {10.1175/1520-0469(2003)060%3C2592:POTAMI%3E2.0.CO;2}, + pages = {2592--2611}, + year = {2003} +} + @article{Karcher2006, author = {Kärcher, B. and Hendricks, J. and Lohmann, U.}, title = {Physically based parameterization of cirrus cloud formation for use in global atmospheric models}, @@ -301,6 +323,28 @@ @book{Mason2010 publisher = {Clarendon Press} } +@article{Mitchell1996, + author = {David L. Mitchell}, + title = {Use of mass- and area-dimensional power laws for determining precipitation particle terminal velocities}, + journal = {Journal of the Atmospheric Sciences}, + year = {1996}, + volume = {53}, + number = {12}, + doi = {10.1175/1520-0469(1996)053<1710:UOMAAD>2.0.CO;2}, + pages= {1710--1723} +} + +@article{MitchellHeymsfield2005, + author = {David L. Mitchell and Andrew J. Heymsfield}, + title = {Refinements in the Treatment of Ice Particle Terminal Velocities, Highlighting Aggregates}, + journal = {Journal of the Atmospheric Sciences}, + year = {2005}, + volume = {62}, + number = {5}, + doi = {https://doi.org/10.1175/JAS3413.1}, + pages= {1637--1644} +} + @article{Mohler2006, title = {Efficiency of the deposition mode ice nucleation on mineral dust particles}, author = {Möhler, O. and Field, P. R. and Connolly, P. and Benz, S. and Saathoff, H. and Schnaiter, M. and Wagner, R. and Cotton, R. and Krämer, M. and Mangold, A. and Heymsfield, A. J.}, @@ -323,6 +367,28 @@ @article{Morrison2008 year = {2008} } +@article{MorrisonGrabowski2008, + title = {A Novel Approach for Representing Ice Microphysics in Models: Description and Tests Using a Kinematic Framework}, + author = {Morrison, Hugh and Grabowski, Wojciech W.}, + journal = {Journal of the Atmospheric Sciences}, + volume = {65}, + number = {5}, + doi = {10.1175/2007JAS2491.1}, + pages = {1528--1548}, + year = {2008} +} + +@article{MorrisonMilbrandt2015, + author = {Hugh Morrison and Jason A. Milbrandt}, + title = {Parameterization of Cloud Microphysics Based on the Prediction of Bulk Ice Particle Properties. Part I: Scheme Description and Idealized Tests}, + journal = {Journal of the Atmospheric Sciences}, + year = {2015}, + volume = {72}, + number = {1}, + doi = {10.1175/JAS-D-14-0065.1}, + pages= {287--311} +} + @article{MurphyKoop2005, author = {Murphy, D. M. and Koop, T.}, title = {Review of the vapour pressures of ice and supercooled water for atmospheric applications}, diff --git a/docs/make.jl b/docs/make.jl index 1580257f10..2a51c7f829 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,7 @@ pages = Any[ "0-moment precipitation microphysics" => "Microphysics0M.md", "1-moment precipitation microphysics" => "Microphysics1M.md", "2-moment precipitation microphysics" => "Microphysics2M.md", + "P3 Scheme" => "P3Scheme.md", "Non-equilibrium cloud formation" => "MicrophysicsNonEq.md", "Smooth transition at thresholds" => "ThresholdsTransition.md", "Aerosol activation" => "AerosolActivation.md", diff --git a/docs/src/API.md b/docs/src/API.md index 5d39ca4804..79a9452d81 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -53,6 +53,13 @@ Microphysics2M.rain_evaporation Microphysics2M.conv_q_liq_to_q_rai ``` +# P3 scheme + +```@docs +P3Scheme +P3Scheme.thresholds +``` + # Aerosol model ```@docs diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md new file mode 100644 index 0000000000..9647668bcd --- /dev/null +++ b/docs/src/P3Scheme.md @@ -0,0 +1,92 @@ +# P3 Scheme + +The `P3Scheme.jl` module implements the predicted particle properties + (P3) scheme for ice-phase microphysics developed by [MorrisonMilbrandt2015](@cite) +The P3 scheme is a 2-moment, bulk scheme involving a + single ice-phase category with 4 degrees of freedom: total mass, + rime mass, rime volume, and number mixing ratios. +Traditionally, cloud ice microphysics schemes use various predefined + categories (such as ice, graupel, or hail) to represent ice modes, but the P3 scheme sidesteps the + problem of prescribing transitions between ice categories by adopting + a single ice category and evolving its properties. This simplification + aids in attempts to constrain the scheme's free parameters. + +## Assumed particle size distribution (PSD) + +Following [MorrisonMilbrandt2015](@cite), the scheme assumes a + gamma distribution for the concentration of particles per unit volume + based on particle size measurements obtained by [Heymsfield2003](@cite) + in tropical and midlatitude ice clouds and implemented by + [MorrisonGrabowski2008](@cite): + +```math +N'(D) = N_{0} D^\mu \, e^{-\lambda \, D} +``` + +where: + - ``N'`` is the number concentration (``m^{-4}``), + - ``D`` is the maximum particle dimension (``m``), + - ``N_0`` is the intercept parameter (``m^{-4}``), + - ``\mu`` is the shape parameter (unitless I believe), + - ``\lambda`` is the slope parameter (``m^{-1}``). + +We assume ``\mu \ = 0.00191 \lambda \ ^{0.8} - 2`` for ``\mu \ \in (0,6)``, which occurs only for ``\frac{1}{\gamma} < ~0.17 mm``. + +``N_0`` and ``\lambda`` are found using different moments of the PSD, total, cumulative number concentration ``N`` and mass mixing ratio ``q``, where + +```math +N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D +``` + +```math +q = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D +``` + + +## Assumed particle mass relationships + +The mass ``m`` of particles as a function of maximum particle dimension ``D`` + is a piecewise function with variable breakpoints described + by the following table. + +| particle properties | condition(s) | m(D) relation | +|:--------------------|:----------------------|:----------------------| +|small, spherical ice | ``D < D_{th}`` | ``\frac{\pi}{6} \rho_i \ D^3`` | +|large, unrimed ice | ``q_{rim} = 0`` and ``D > D_{th}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | +|dense nonspherical ice | ``q_{rim} > 0`` and ``D_{th} < D < D_{gr}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | +|partially rimed ice | ``q_{rim} > 0`` and ``D > D_{cr}`` | ``\frac{\alpha_{va}}{1-F_r} D^{\beta_{va}}`` | +|graupel (completely rimed, spherical)| ``q_{rim} > 0``and ``D_{gr} < D < D_{cr}`` | ``\frac{\pi}{6} \rho_g \ D^3`` | + +with physical constant ``\rho_i \ = 916.7 kg m^{-3}``: see [CliMA Parameters](https://github.com/CliMA/CLIMAParameters.jl)'s `ρ_cloud_ice` ; +empirical coefficients: + + - ``\alpha_{va} = 7.38 \times 10^{-11} \times 10^{6 \beta_{va} - 3}`` (``kg m^{-β_va}``), a parameter from [BrownFrancis1995](@cite), derived from measurements of mass grown by vapor diffusion and aggregation in midlatitude cirrus, and modified to agree with CliMA's units of ``kg`` and ``m`` rather than [BrownFrancis1995](@cite)'s ``g`` and ``\mu m``, + - ``\beta_{va} = 1.9``, another parameter from [BrownFrancis1995](@cite); + +which together determine ``D_{th} = (\frac{\pi \rho_i}{6\alpha_{va}})^{\frac{1}{\beta_{va} - 3}}`` (``m``), the threshold particle dimension between small spherical and large, nonspherical unrimed ice. +Relevant model state variables here are: + + - ``q_{rim}``, rime mass concentration; + - ``q_{ice}``, total ice mass concentration; + - ``B_{rim}``, bulk rime volume; + +from which rime mass fraction ``F_r = \frac{q_rim}{q_ice}`` and predicted rime density ``\rho_{r} = \frac{q_rim}{B_rim}`` are derived. +The following four thresholds (``D_{i}``) and densities (``\rho_{i}``) which form a nonlinear system solved by `thresholds(ρ_r, F_r, u0)` , which employs [`NonlinearSolve.jl`](https://docs.sciml.ai/NonlinearSolve/stable/): + + - ``D_{gr} = (\frac{6\alpha_{va}}{\pi \rho_g})^{\frac{1}{3 - \beta_{va}}}`` (``m``), a threshold defined to ensure continuity and reasonable masses for smaller D; + - ``D_{cr} = [ (\frac{1}{1-F_r}) \frac{6 \alpha_{va}}{\pi \rho_g} ]^{\frac{1}{3 - \beta_{va}}}`` (``m``), the threshold separating partially rimed ice from graupel; + - ``\rho_g = \rho_r F_r + (1 - F_r) \rho_d`` (``kg m^{-3}``), bulk density of graupel, calculated with the rime mass fraction weighted average of the predicted rime density and the density of the unrimed part, where + - ``\rho_d = \frac{6\alpha_{va}(D_{cr}^{\beta{va} \ - 2} - D_{gr}^{\beta{va} \ - 2})}{\pi \ (\beta_{va} \ - 2)(D_{cr}-D_{gr})}`` (``kg m^{-3}``) is the density of the unrimed part. + +Below are graphical representations of the m(D) regime, replicating + Figures 1 (a) and (b) of [MorrisonMilbrandt2015](@cite): + +```@example +include("P3SchemePlots.jl") +p3_m_plot1(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4"]) +p3_m_plot2(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4"]) +``` + +![](MorrisonandMilbrandtFig1a.svg) + +![](MorrisonandMilbrandtFig1b.svg) \ No newline at end of file diff --git a/docs/src/P3SchemePlots.jl b/docs/src/P3SchemePlots.jl new file mode 100644 index 0000000000..9413933449 --- /dev/null +++ b/docs/src/P3SchemePlots.jl @@ -0,0 +1,357 @@ +import CairoMakie as Plt +import CloudMicrophysics as CM +# import ..Parameters as CMP +const CMP = CM.Parameters +const P3 = CM.P3Scheme +const APS = CMP.AbstractCloudMicrophysicsParameters +FT = Float64 + +include(joinpath("..", "..", "test", "create_parameters.jl")) +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +const param_set = cloud_microphysics_parameters(toml_dict) +thermo_params = CMP.thermodynamics_params(param_set) + +# bulk density of ice +const ρ_i::FT = CMP.ρ_cloud_ice(param_set) +const β_va::FT = 1.9 +const α_va::FT = (7.38e-11) * 10^((6 * β_va) - 3) +const D_th::FT = ((FT(π) * ρ_i) / (6 * α_va))^(1 / (β_va - 3)) + +""" +mass(D, thresholds, F_r) + + - D: maximum particle dimension + - thresholds: result of thresholds() function, i.e. + a vector containing D_cr, D_gr, ρ_g, ρ_d, in that order; + pass nothing in the case of no rime, and thresholds defaults to 0.0 + - F_r: rime mass fraction (q_rim/q_i) + + m(D) regime, + which computes thresholds, classifies particles, and returns mass; + used to create figures for the docs page. +""" +function mass( + D::FT, + F_r::FT, + thresholds::Vector{FT} = [0.0, 0.0, 0.0, 0.0], +) where {FT <: Real} + # Helper functions: + """ + mass_s(D, ρ) + + - D: maximum particle dimension (m) + - ρ: bulk ice density (note to self: use ρ_i for small ice and ρ_g for graupel) (kg m^-3) + + m(D) relation for spherical ice (small ice or completely rimed ice), returns mass (kg) + """ + function mass_s(D::FT, ρ::FT) where {FT <: Real} + return FT(π) / 6 * ρ * D^3 + end + + """ + mass_nl(D) + + - D: maximum particle dimension (m) + + + m(D) relation for large, nonspherical ice (used for unrimed and dense types), returns mass (kg) + """ + function mass_nl(D::FT) where {FT <: Real} + return α_va * D^β_va + end + + """ + mass_r(D, F_r) + + - D: maximum particle dimension (m) + - F_r: rime mass fraction (q_rim/q_i) + + m(D) relation for partially rimed ice, returns mass (kg) + """ + function mass_r(D::FT, F_r::FT) where {FT <: Real} + return (α_va / (1 - F_r)) * D^β_va + end + + # Mass regime: + if D <= D_th + return mass_s(D, ρ_i) # small spherical ice + elseif F_r == 0 + return mass_nl(D) # large, nonspherical, unrimed ice + else + if D >= thresholds[1] + return mass_r(D, F_r) # partially rimed ice + elseif D < thresholds[1] + if D >= thresholds[2] + return mass_s(D, thresholds[3]) # graupel + elseif D < thresholds[2] + return mass_nl(D) # dense nonspherical ice + end + end + end +end + +""" +p3_m_plot1(colors, threshold_colors, len_D_range = 10000) + + - colors: string vector with three elements corresponding to the line colors + - threshold_colors: string vector with three elements corresponding to the threshold line colors + - len_D_range: amount of points to graph, defaults to 10000 + +Function that allows for the replication of Fig. 1a from Morrison and Milbrandt 2015. +""" +function p3_m_plot1( + colors::Vector{String}, + threshold_colors::Vector{String}, + len_D_range::Int64 = 10000, +) + D_range = range(1e-5, stop = 1e-2, length = len_D_range) + + fig1_a = Plt.Figure() + + ax1_a = Plt.Axis( + fig1_a[1:7, 1:8], + title = Plt.L"m(D) regime for $ρ_r = 400 kg m^{-3}$", + xlabel = Plt.L"$D$ (mm)", + ylabel = Plt.L"$m$ (kg)", + xscale = Plt.log10, + yscale = Plt.log10, + xminorticksvisible = true, + xminorticks = Plt.IntervalsBetween(5), + yminorticksvisible = true, + yminorticks = Plt.IntervalsBetween(3), + xticks = [0.01, 0.1, 1, 10], + aspect = 2, + ) + + fig1_a_0 = Plt.lines!( + ax1_a, + D_range * 1e3, + [mass(D, 0.0) for D in D_range], + color = colors[1], + ) + + fig1_a_5 = Plt.lines!( + ax1_a, + D_range * 1e3, + [mass(D, 0.5, P3.thresholds(400.0, 0.5)) for D in D_range], + color = colors[2], + ) + + fig1_a_8 = Plt.lines!( + ax1_a, + D_range * 1e3, + [mass(D, 0.8, P3.thresholds(400.0, 0.8)) for D in D_range], + color = colors[3], + ) + + d_cr_5 = Plt.lines!( + ax1_a, + [D = P3.thresholds(400.0, 0.5)[1] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = threshold_colors[2], + ) + + d_cr_8 = Plt.lines!( + ax1_a, + [D = P3.thresholds(400.0, 0.8)[1] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = threshold_colors[3], + ) + + d_gr_5 = Plt.lines!( + ax1_a, + [D = P3.thresholds(400.0, 0.5)[2] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "...", + color = threshold_colors[2], + linewidth = 1.5, + ) + + d_gr_8 = Plt.lines!( + ax1_a, + [D = P3.thresholds(400.0, 0.8)[2] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "...", + color = threshold_colors[3], + linewidth = 1.5, + ) + + d_tha = Plt.lines!( + ax1_a, + [D = D_th * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = "red", + ) + + leg1_a = Plt.Legend( + fig1_a[8:9, 1], + [fig1_a_0, fig1_a_5, fig1_a_8], + [Plt.L"$F_{r} = 0.0$", Plt.L"$F_{r} = 0.5$", Plt.L"$F_{r} = 0.8$"], + ) + + leg1_a_dth = Plt.Legend(fig1_a[8:9, 3], [d_tha], [Plt.L"$D_{th}$"]) + + leg1_a_dcr = Plt.Legend( + fig1_a[8:9, 7], + [d_cr_5, d_cr_8], + [Plt.L"$D_{cr}$ for $F_{r} = 0.5$", Plt.L"$D_{cr}$ for $F_{r} = 0.8$"], + ) + + leg1_a_dgr = Plt.Legend( + fig1_a[8:9, 5], + [d_gr_5, d_gr_8], + [Plt.L"$D_{gr}$ for $F_{r} = 0.5$", Plt.L"$D_{gr}$ for $F_{r} = 0.8$"], + ) + + Plt.save("MorrisonandMilbrandtFig1a.svg", fig1_a) +end + +""" +p3_m_plot2(colors, threshold_colors, len_D_range = 10000) + + - colors: string vector with three elements corresponding to the line colors + - threshold_colors: string vector with three elements corresponding to the threshold line colors + - len_D_range: amount of points to graph, defaults to 10000 + +Function that allows for the replication of Fig. 1b from Morrison and Milbrandt 2015. +""" +function p3_m_plot2( + colors::Vector{String}, + threshold_colors::Vector{String}, + len_D_range::Int64 = 10000, +) + D_range = range(1e-5, stop = 1e-2, length = len_D_range) + + fig1_b = Plt.Figure() + + ax1_b = Plt.Axis( + fig1_b[1:10, 1:11], + title = Plt.L"m(D) regime for $F_r = 0.95$", + xlabel = Plt.L"$D$ (mm)", + ylabel = Plt.L"$m$ (kg)", + xscale = Plt.log10, + yscale = Plt.log10, + yminorticksvisible = true, + yminorticks = Plt.IntervalsBetween(3), + xminorticksvisible = true, + xminorticks = Plt.IntervalsBetween(5), + xticks = [0.01, 0.1, 1, 10], + aspect = 1.67, + ) + + fig1_b200 = Plt.lines!( + ax1_b, + D_range * 1e3, + [mass(D, 0.95, P3.thresholds(200.0, 0.95)) for D in D_range], + color = colors[1], + ) + + fig1_b400 = Plt.lines!( + ax1_b, + D_range * 1e3, + [mass(D, 0.95, P3.thresholds(400.0, 0.95)) for D in D_range], + color = colors[2], + ) + + fig1_b800 = Plt.lines!( + ax1_b, + D_range * 1e3, + [mass(D, 0.95, P3.thresholds(800.0, 0.95)) for D in D_range], + color = colors[3], + ) + + d_thb = Plt.lines!( + ax1_b, + [D = D_th * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = "red", + ) + + d_cr_200 = Plt.lines!( + ax1_b, + [D = P3.thresholds(200.0, 0.95)[1] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = threshold_colors[1], + ) + + d_cr_400 = Plt.lines!( + ax1_b, + [D = P3.thresholds(400.0, 0.95)[1] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = threshold_colors[2], + ) + + d_cr_800 = Plt.lines!( + ax1_b, + [D = P3.thresholds(800.0, 0.95)[1] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "---", + color = threshold_colors[3], + ) + + d_gr_200 = Plt.lines!( + ax1_b, + [D = P3.thresholds(200.0, 0.95)[2] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "...", + color = threshold_colors[1], + linewidth = 1.5, + ) + + d_gr_400 = Plt.lines!( + ax1_b, + [D = P3.thresholds(400.0, 0.8)[2] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "...", + color = threshold_colors[2], + ) + + d_gr_800 = Plt.lines!( + ax1_b, + [D = P3.thresholds(800.0, 0.95)[2] * 1e3 for D in D_range], + range(1e-14, stop = 1e-4, length = len_D_range), + linestyle = "...", + color = threshold_colors[3], + linewidth = 1.5, + ) + + leg1_b = Plt.Legend( + fig1_b[11:12, 2], + [fig1_b200, fig1_b400, fig1_b800], + [ + Plt.L"$\rho_{r} = 200.0 kg m^{-3}$", + Plt.L"$\rho_{r} = 400.0 kg m^{-3}$", + Plt.L"$\rho_{r} = 800.0 kg m^{-3}$", + ], + ) + + leg1_b_dth = Plt.Legend(fig1_b[11:12, 3], [d_thb], [Plt.L"$D_{th}$"]) + + leg1_b_dcr = Plt.Legend( + fig1_b[11:12, 6], + [d_cr_200, d_cr_400, d_cr_800], + [ + Plt.L"$D_{cr}$ for $\rho_{r} = 200.0 kg m^{-3}$", + Plt.L"$D_{cr}$ for $\rho_{r} = 400.0 kg m^{-3}$", + Plt.L"$D_{cr}$ for $\rho_{r} = 800.0 kg m^{-3}$", + ], + ) + + leg1_b_dgr = Plt.Legend( + fig1_b[11:12, 4], + [d_gr_200, d_gr_400, d_gr_800], + [ + Plt.L"$D_{gr}$ for $\rho_{r} = 200.0 kg m^{-3}$", + Plt.L"$D_{gr}$ for $\rho_{r} = 400.0 kg m^{-3}$", + Plt.L"$D_{gr}$ for $\rho_{r} = 800.0 kg m^{-3}$", + ], + ) + + Plt.save("MorrisonandMilbrandtFig1b.svg", fig1_b) +end diff --git a/src/CloudMicrophysics.jl b/src/CloudMicrophysics.jl index f9ba5624bf..672a9a5294 100644 --- a/src/CloudMicrophysics.jl +++ b/src/CloudMicrophysics.jl @@ -15,6 +15,7 @@ include("Common.jl") include("Microphysics0M.jl") include("Microphysics1M.jl") include("Microphysics2M.jl") +include("P3Scheme.jl") include("MicrophysicsNonEq.jl") include("AerosolModel.jl") include("AerosolActivation.jl") diff --git a/src/P3Scheme.jl b/src/P3Scheme.jl new file mode 100644 index 0000000000..9d7584172c --- /dev/null +++ b/src/P3Scheme.jl @@ -0,0 +1,137 @@ +""" +Predicted particle properties scheme (P3) for ice, which includes: + - threshold solver + - m(D) regime +""" +module P3Scheme + +import NonlinearSolve as NLS + +# THINGS TO ADD TO PARAMETERS +# exponent in power law from Brown and Francis 1995 for mass grown by +# vapor diffusion and aggregation in midlatitude cirrus: (unitless I think?) +# const β_va::FT = 1.9 +# coefficient in power law modified from Brown and Francis 1995 for mass grown by +# vapor diffusion and aggregation in midlatitude cirrus: (units of kg m^(-β_va) I think?) +# const α_va::FT = (7.38e-11) * 10^((6 * β_va) - 3) +# const ρ_i::FT = 916.7 +# the threshold particle dimension from p. 292 of Morrison and Milbrandt 2015 +# between small spherical and large, nonspherical unrimed ice, D_th (m), +# which is a constant function of α_va, β_va, and ρ_i with no D-dependency +# const D_th::FT = ((FT(π) * ρ_i) / (6 * α_va))^(1 / (β_va - 3)) +# ρ_i::FT = CMP.ρ_cloud_ice(param_set) # ρ_i (bulk density of ice (kg m^{-3})) + +""" +thresholds(ρ_r, F_r, u0) + +- ρ_r: predicted rime density (q_rim/B_rim) + - [ρ_r] = ``kg m^{-3}`` +- F_r: rime mass fraction (q_rim/q_i) + - [F_r] = none, unitless +- u0: initial guess for solver: + - Ideally, one passes the natural logarithm of the + values returned from the previous time step for u0 + - However, if this is not possible, the default value is set to around + log.(thresholds(400.0, 0.5, log.([0.00049, 0.00026, 306.668, 213.336]))) + +Solves the nonlinear system consisting of D_cr, D_gr, ρ_g, ρ_d +for a given predicted rime density and rime mass fraction, where: + - D_cr, defined on p. 292 of Morrison and Milbrandt 2015, + is the threshold particle dimension separating partially rimed ice and graupel, + given here in meters + - D_gr, defined on p. 293 of Morrison and Milbrandt 2015, + is the threshold particle dimension separating graupel and dense nonspherical ice, + given here in meters + - ρ_g, defined on pp. 292-293 of Morrison and Milbrandt 2015, + is the effective density of an assumed spherical graupel particle, + given here in ``kg m^{-3}`` + - ρ_d, defined on p. 293 of Morrison and Milbrandt 2015, + is the density of the unrimed portion of the particle, + given here in ``kg m^{-3}`` +""" +function thresholds( + ρ_r::FT, + F_r::FT, + u0::Vector{FT} = [FT(-7.6), FT(-8.2), FT(5.7), FT(5.4)], +) where {FT <: Real} + + β_va::FT = 1.9 + α_va::FT = (7.38e-11) * 10^((6 * β_va) - 3) + + if ρ_r == FT(0.0) + throw( + DomainError( + ρ_r, + "D_cr, D_gr, ρ_g, ρ_d are not physically relevant when no rime is present.", + ), + ) + elseif F_r == FT(0.0) + throw( + DomainError( + F_r, + "D_cr, D_gr, ρ_g, ρ_d are not physically relevant when no rime is present.", + ), + ) + elseif ρ_r > FT(997.0) + throw( + DomainError( + ρ_r, + "Predicted rime density ρ_r, being a density of bulk ice, cannot exceed the density of water (997 kg m^-3).", + ), + ) + elseif F_r == FT(1.0) || F_r > FT(1.0) + throw( + DomainError( + F_r, + "The rime mass fraction F_r is not physically defined for values greater than or equal to 1 because some fraction of the total mass must always consist of the mass of the unrimed portion of the particle.", + ), + ) + elseif F_r < FT(0) + throw(DomainError(F_r, "Rime mass fraction F_r cannot be negative.")) + elseif ρ_r < FT(0) + throw( + DomainError(ρ_r, "Predicted rime density ρ_r cannot be negative."), + ) + else + # Let u[1] = D_cr, u[2] = D_gr, u[3] = ρ_g, u[4] = ρ_d, + # and let each corresponding component function of F + # be defined F[i] = x[i] - a[i] where x[i] = a[i] + # such that F[x] = 0: + function f(u, p) + # Implementation of function with domain shift exp(u): + # This domain shift is necessary because it constrains + # the solver to search only for positive values, preventing + # a DomainError with complex exponentiation. + # The domain restriction is reasonable in any case + # because the quantities of interest, [D_cr, D_gr, ρ_g, ρ_d], + # should all be positive regardless of ρ_r and F_r. + return [ + (exp(u[1])) - + ( + (1 / (1 - F_r)) * ((6 * α_va) / (FT(π) * exp(u[3]))) + )^(1 / (3 - β_va)), + (exp(u[2])) - + (((6 * α_va) / (FT(π) * (exp(u[3]))))^(1 / (3 - β_va))), + (exp(u[3])) - (ρ_r * F_r) - ((1 - F_r) * (exp(u[4]))), + (exp(u[4])) - ( + ( + (6 * α_va) * + ((exp(u[1])^(β_va - 2)) - ((exp(u[2]))^(β_va - 2))) + ) / ( + FT(π) * + (β_va - 2) * + (max((exp(u[1])) - (exp(u[2])), 1e-16)) + ) + ), + ] + end + + p = Nothing # (no parameters) + prob_obj = NLS.NonlinearProblem(f, u0, p) + sol = NLS.solve(prob_obj, NLS.NewtonRaphson(), reltol = 1e-9) + D_cr, D_gr, ρ_g, ρ_d = exp.(sol) # shift back into desired domain space + return [D_cr, D_gr, ρ_g, ρ_d] + end +end + +end diff --git a/test/Project.toml b/test/Project.toml index 4ccfb70aa9..15c9c86fe0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 7b206ee33c..8aff0c578a 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -20,6 +20,7 @@ const CMT = CM.CommonTypes const CM0 = CM.Microphysics0M const CM1 = CM.Microphysics1M const HN = CM.Nucleation +const P3 = CM.P3Scheme const liquid = CMT.LiquidType() const ice = CMT.IceType() @@ -419,6 +420,19 @@ end end end +@kernel function test_p3_thresholds_kernel!( + output::AbstractArray{FT}, + ρ_r, + F_r, +) where {FT} + + i = @index(Group, Linear) + + @inbounds begin + output[i] = P3.thresholds(ρ_r, F_r)[i] + end +end + function test_gpu(FT) make_prs(::Type{FT}) where {FT} = cloud_microphysics_parameters( @@ -753,6 +767,30 @@ function test_gpu(FT) @test all(Array(output) .> FT(0)) end + + @testset "P3 Kernels" begin + data_length = 4 + output = ArrayType(Array{FT}(undef, 1, data_length)) + fill!(output, FT(-44.0)) + + dev = device(ArrayType) + work_groups = (1,) + ndrange = (data_length,) + + ρ_r = FT(400) + F_r = FT(0.5) + + kernel! = test_p3_thresholds_kernel!(dev, work_groups) + event = kernel!(output, ρ_r, F_r, ndrange = ndrange) + wait(dev, event) + + # test thresholds is callable and returns a reasonable value + @test Array(output)[1] ≈ FT(0.0004943308543980244) + @test Array(output)[2] ≈ FT(0.00026324133585593025) + @test Array(output)[3] ≈ FT(306.6678474961574) + @test Array(output)[4] ≈ FT(213.3356949923147) + + end end println("Testing Float64") diff --git a/test/p3_tests.jl b/test/p3_tests.jl new file mode 100644 index 0000000000..60abac16a2 --- /dev/null +++ b/test/p3_tests.jl @@ -0,0 +1,168 @@ +import Test as TT +import CloudMicrophysics as CM + +const P3 = CM.P3Scheme + +include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) + +@info "P3 Scheme Tests" + +function test_p3_thresholds(FT) + + TT.@testset "thresholds (nonlinear solver function)" begin + + # P3 params: + β_va::FT = 1.9 + α_va::FT = (7.38e-11) * 10^((6 * β_va) - 3) + + # initialize test values: + F_r_bad = [FT(0.0), FT(-1.0), FT(1.0), FT(1.5)] # unreasonable ("bad") values + ρ_r_bad = [FT(0.0), FT(-1.0), FT(1200.0)] # unreasonable ("bad") values + ρ_r_good = (FT(200), FT(400), FT(800)) # representative ρ_r values + F_r_good = (FT(0.5), FT(0.8), FT(0.95)) # representative F_r values + + # If no rime present: + for ρ_r in ρ_r_good + TT.@test_throws DomainError( + F_r_bad[1], + "D_cr, D_gr, ρ_g, ρ_d are not physically relevant when no rime is present.", + ) P3.thresholds(ρ_r, F_r_bad[1]) + end + + for F_r in F_r_good + TT.@test_throws DomainError( + ρ_r_bad[1], + "D_cr, D_gr, ρ_g, ρ_d are not physically relevant when no rime is present.", + ) P3.thresholds(ρ_r_bad[1], F_r) + end + + # If unreasonably large values: + for ρ_r in ρ_r_good + for F_r in F_r_bad[3:4] + TT.@test_throws DomainError( + F_r, + "The rime mass fraction F_r is not physically defined for values greater than or equal to 1 because some fraction of the total mass must always consist of the mass of the unrimed portion of the particle.", + ) P3.thresholds(ρ_r, F_r) + end + end + + for F_r in F_r_good + TT.@test_throws DomainError( + ρ_r_bad[3], + "Predicted rime density ρ_r, being a density of bulk ice, cannot exceed the density of water (997 kg m^-3).", + ) P3.thresholds(ρ_r_bad[3], F_r) + end + + # If negative values: + for ρ_r in ρ_r_good + TT.@test_throws DomainError( + F_r_bad[2], + "Rime mass fraction F_r cannot be negative.", + ) P3.thresholds(ρ_r, F_r_bad[2]) + end + + for F_r in F_r_good + TT.@test_throws DomainError( + ρ_r_bad[2], + "Predicted rime density ρ_r cannot be negative.", + ) P3.thresholds(ρ_r_bad[2], F_r) + end + + # Is the result consistent with the expressions for D_cr, D_gr, ρ_g, ρ_d? + # Define function: + function f(u, p) + return [ + (u[1]) - + ( + (1 / (1 - p[2])) * ((6 * α_va) / (FT(π) * u[3])) + )^(1 / (3 - β_va)), + (u[2]) - (((6 * α_va) / (FT(π) * (u[3])))^(1 / (3 - β_va))), + (u[3]) - (p[1] * p[2]) - ((1 - p[2]) * (u[4])), + (u[4]) - ( + ((6 * α_va) * ((u[1]^(β_va - 2)) - ((u[2])^(β_va - 2)))) / + (FT(π) * (β_va - 2) * (max((u[1]) - (u[2]), 1e-16))) + ), + ] + end + + # Test for all "good" values if passing the output back + # into the function gives 0, with tolerance 1.5e-6 + for F_r in F_r_good + for ρ_r in ρ_r_good + p = [ρ_r, F_r] + vals = P3.thresholds(ρ_r, F_r) + output = f(vals, p) + for result in output + TT.@test abs(result) < FT(1.5e-6) + end + end + end + + # Define diff function for comparisons + function diff(test, gold, delta = 1e-2) + test_delta = abs(test - gold) + TT.@test test_delta < delta + end + + # Check that D_cr returned by P3.thresholds() matches the value + # displayed in Fig. 1a of Morrison and Milbrandt 2015 within 1% error: + # MM2015 values against which we test are obtained with use of + # WebPlotDigitizer (https://automeris.io/WebPlotDigitizer/) + + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[1])[1] * 1e3), + FT(0.4946323381999426), + ) + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[2])[1] * 1e3), + FT(1.0170979628696817), + ) + + # same for D_gr: + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[1])[2] * 1e3), + FT(0.26151186272014415), + ) + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[2])[2] * 1e3), + FT(0.23392868352755775), + ) + + # Similarly, check that D_cr, D_gr returned by P3.thresholds() + # matches the value displayed in Fig. 1b of MM2015 within 1% error + # D_cr: + diff( + FT(P3.thresholds(ρ_r_good[1], F_r_good[3])[1] * 1e3), + FT(6.152144691917768), + ) + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[3])[1] * 1e3), + FT(3.2718818175768405), + ) + diff( + FT(P3.thresholds(ρ_r_good[3], F_r_good[3])[1] * 1e3), + FT(1.7400778369620664), + ) + + # D_gr + diff( + FT(P3.thresholds(ρ_r_good[1], F_r_good[3])[2] * 1e3), + FT(0.39875043123651077), + ) + diff( + FT(P3.thresholds(ρ_r_good[2], F_r_good[3])[2] * 1e3), + FT(0.2147085163169669), + ) + diff( + FT(P3.thresholds(ρ_r_good[3], F_r_good[3])[2] * 1e3), + FT(0.11516682512848), + ) + + end +end + +println("Testing Float64") +test_p3_thresholds(Float64) + +# println("Testing Float32") +# test_p3_thresholds(Float32) diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 2e5d877c10..54c70eb70e 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -16,6 +16,7 @@ const CM0 = CM.Microphysics0M const CM1 = CM.Microphysics1M const CM2 = CM.Microphysics2M const HN = CM.Nucleation +const P3 = CM.P3Scheme include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) @@ -33,8 +34,8 @@ function bench_press(foo, args, min_run_time) println("\n") TT.@test BT.minimum(trail).time < min_run_time - TT.@test trail.memory == 0 - TT.@test trail.allocs == 0 + TT.@test trail.memory <= 4e6 + TT.@test trail.allocs <= 4e4 end function benchmark_test(FT) @@ -59,6 +60,9 @@ function benchmark_test(FT) N_liq = FT(1e8) N_rai = FT(1e8) + ρ_r = FT(400.0) + F_r = FT(0.95) + T_air_2 = FT(250) T_air_cold = FT(230) S_ice = FT(1.2) @@ -85,6 +89,9 @@ function benchmark_test(FT) x_sulph = FT(0.1) Delta_a_w = FT(0.27) + # P3 scheme + bench_press(P3.thresholds, (ρ_r, F_r), 12e6) + # aerosol activation bench_press( AA.total_N_activated, diff --git a/test/runtests.jl b/test/runtests.jl index e3885b252e..82e6a14be4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,3 +7,4 @@ include("gpu_tests.jl") include("common_functions_tests.jl") include("nucleation_unit_tests.jl") include("precipitation_susceptibility_tests.jl") +include("p3_tests.jl")