diff --git a/Project.toml b/Project.toml index 2405ae423..74f8f4368 100644 --- a/Project.toml +++ b/Project.toml @@ -13,7 +13,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" [compat] -CLIMAParameters = ">=0.7.12" +CLIMAParameters = ">=0.7.13" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" KernelAbstractions = "0.8, 0.9" diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index fa2bf26a2..cacbad509 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -107,3 +107,37 @@ p3_m_plot2(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4" ![](MorrisonandMilbrandtFig1a.svg) ![](MorrisonandMilbrandtFig1b.svg) + +## Assumed particle projected area relationships + +The projected area ``A`` of particles as a function of maximum particle dimension ``D`` + is another piecewise function with variable breakpoints described + by the following table. The mean particle dimension ``D_m``, a predicted property, + determines which portion of the piecewise function to use for each time step. + +| particle properties | condition(s) | A(D) relation | +|:------------------------------------|:--------------------------------------------|:-----------------------------------------------------------| +|small, spherical ice | ``D < D_{th}`` | ``\frac{\pi}{4} D^2`` | +|large, unrimed ice | ``q_{rim} = 0`` and ``D > D_{th}`` | ``\gamma \ D^{\sigma}`` | +|dense nonspherical ice | ``q_{rim} > 0`` and ``D_{th} < D < D_{gr}`` | ``\gamma \ D^{\sigma}`` | +|graupel (completely rimed, spherical)| ``q_{rim} > 0``and ``D_{gr} < D < D_{cr}`` | ``\frac{\pi}{4} D^2`` | +|partially rimed ice | ``q_{rim} > 0`` and ``D > D_{cr}`` | ``F_{r} \frac{\pi}{4} D^2 + (1-F_{r})\gamma \ D^{\sigma}`` | + +where all variables from the m(D) regime are as defined above, and: + - ``\gamma = 0.2285`` (``m^{2 - \sigma}``) where + - ``\sigma = 1.88`` (dimensionless), both from the aggregates of side planes, columns, bullets, and planar polycrystals in [Mitchell1996](@cite). + +!!! note + TODO - As mentioned in [issue #151](https://github.com/CliMA/CloudMicrophysics.jl/issues/151), the units of ``\gamma`` and ``\sigma`` are not immediately clear from [Mitchell1996](@cite) and [MorrisonMilbrandt2015](@cite). To resolve this issue, it may be useful to contact the authors of the paper, or, examine the representative figures below to check units. It has occured to me that the units of ``D`` are probably m and that the units of ``A`` are probably m2. I have assumed these dimensions for the time being. Another likely scenario would be if ``D`` had units of mm, in which case we would have ``\gamma = 0.2285 \; 10^{3 \sigma}`` to correct for units. However, the plots of area versus particle dimension look outlandish in this case. + +The figures below, imitating the above mass regime figures, illustrate the particle projected area regime: + +```@example +include("P3SchemePlots.jl") +p3_a_plot1(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4"]) +p3_a_plot2(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4"]) +``` + +![](P3Scheme_Area_1.svg) + +![](P3Scheme_Area_2.svg) diff --git a/docs/src/P3SchemePlots.jl b/docs/src/P3SchemePlots.jl index d5c1e4338..d32072836 100644 --- a/docs/src/P3SchemePlots.jl +++ b/docs/src/P3SchemePlots.jl @@ -17,6 +17,10 @@ const α_va::FT = P3.α_va(prs) const β_va = CMP.β_va_BF1995(prs) # threshold particle dimension const D_th::FT = P3.D_th(prs, FT) +# area parameters +const γ::FT = CMP.γ_M1996(prs) +const σ::FT = CMP.σ_M1996(prs) + """ mass(D, thresholds, F_r) @@ -28,7 +32,7 @@ mass(D, thresholds, F_r) - F_r: rime mass fraction (q_rim/q_i) m(D) regime, - which computes thresholds, classifies particles, and returns mass; + which uses thresholds, classifies particles, and returns mass; used to create figures for the docs page. """ function mass( @@ -91,6 +95,65 @@ function mass( end end +""" +area(D, F_r, thresholds) + +- D: maximum particle dimension +- F_r: rime mass fraction (q_rim/q_i) +- 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 + +Area regime, which uses thresholds, classifies particles, and returns particle projected area; +used to create figures for the documentation. +""" +function area( + D::FT, + F_r::FT, + thresholds::Vector{FT} = [0.0, 0.0, 0.0, 0.0], +) where {FT <: Real} + """ + A_s(D) + + + - D: maximum particle dimension + + Particle projected area relation for assumed spherical particles + """ + function A_s(D::FT) where {FT <: Real} + return (FT(π) / 4) * D^2 + end + + """ + A_ns(D) + + - D: maximum particle dimension + + Particle projected area relation for assumed nonspherical particles, + from Mitchell 1996 + """ + function A_ns(D::FT) where {FT <: Real} + return γ * D^σ + end + + # Area regime: + if D <= D_th + return A_s(D) # small spherical ice + elseif F_r == 0 + return A_ns(D) # large, nonspherical, unrimed ice + else + if D >= thresholds[1] + return (F_r * A_s(D)) + (1 - F_r) * A_ns(D) # partially rimed ice + elseif D < thresholds[1] + if D >= thresholds[2] + return A_s(D) # graupel + elseif D < thresholds[2] + return A_ns(D) # dense nonspherical ice + end + end + end +end + """ p3_m_plot1(colors, threshold_colors, len_D_range = 10000) @@ -374,3 +437,286 @@ function p3_m_plot2( Plt.save("MorrisonandMilbrandtFig1b.svg", fig1_b) end + +""" +p3_a_plot1(colors, threshold_colors, len_D_range = 5000) + + - 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 5000 + + Function that creates a visual representation of the area regime for the documentation. +""" +function p3_a_plot1( + colors::Vector{String}, + threshold_colors::Vector{String}, + len_D_range::Int64 = 5000, +) + D_range = range(0.00005, stop = 0.003, length = len_D_range) + + lw = 3 + + fig = Plt.Figure() + + ax = Plt.Axis( + fig[1:7, 1:8], + title = Plt.L"A(D) regime for $ρ_r = 400 kg m^{-3}$", + xlabel = Plt.L"$D$ (mm)", + ylabel = Plt.L"$A$ ($m^2$)", + xscale = Plt.log10, + yscale = Plt.log10, + yminorticksvisible = true, + yminorticks = Plt.IntervalsBetween(2), + xticks = [0.1, 1], + xminorticksvisible = true, + xminorticks = Plt.IntervalsBetween(4), + aspect = 2, + ) + + fig_0 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.0) for D in D_range], + color = colors[1], + linewidth = lw, + ) + + fig_5 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.5, P3.thresholds(400.0, 0.5)) for D in D_range], + color = colors[2], + linewidth = lw, + ) + + fig_8 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.8, P3.thresholds(400.0, 0.8)) for D in D_range], + color = colors[3], + linewidth = lw, + ) + + d_cr_5 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.5)[1] * 1e3 for D in D_range], + range(1e-10, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = threshold_colors[2], + linewidth = lw, + ) + + d_cr_8 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.8)[1] * 1e3 for D in D_range], + range(1e-10, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = threshold_colors[3], + linewidth = lw, + ) + + d_gr_5 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.5)[2] * 1e3 for D in D_range], + range(1e-10, stop = 1e-5, length = len_D_range), + linestyle = "...", + color = threshold_colors[2], + linewidth = lw, + ) + + d_gr_8 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.8)[2] * 1e3 for D in D_range], + range(1e-10, stop = 1e-5, length = len_D_range), + linestyle = "...", + color = threshold_colors[3], + linewidth = lw, + ) + + d_th = Plt.lines!( + ax, + [D = D_th * 1e3 for D in D_range], + range(1e-10, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = "red", + linewidth = lw, + ) + + leg = Plt.Legend( + fig[8:9, 1], + [fig_0, fig_5, fig_8], + [Plt.L"$F_{r} = 0.0$", Plt.L"$F_{r} = 0.5$", Plt.L"$F_{r} = 0.8$"], + ) + + leg_dth = Plt.Legend(fig[8:9, 3], [d_th], [Plt.L"$D_{th}$"]) + + leg_dcr = Plt.Legend( + fig[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$"], + ) + + leg_dgr = Plt.Legend( + fig[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("P3Scheme_Area_1.svg", fig) +end + +""" +p3_a_plot2(colors, threshold_colors, len_D_range = 11000) + + - 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 11000 + +Function that creates a visual representation of the area regime for the documentation. +""" +function p3_a_plot2( + colors::Vector{String}, + threshold_colors::Vector{String}, + len_D_range::Int64 = 11000, +) + D_range = range(0.00005, stop = 0.003, length = len_D_range) + + fig = Plt.Figure() + + lw = 3 + + ax = Plt.Axis( + fig[1:10, 1:11], + title = Plt.L"A(D) regime for $F_r = 0.5$", + xlabel = Plt.L"$D$ (mm)", + ylabel = Plt.L"$A$ ($m^2$)", + xscale = Plt.log10, + yscale = Plt.log10, + yminorticksvisible = true, + yminorticks = Plt.IntervalsBetween(2), + xticks = [0.1, 1], + xminorticksvisible = true, + xminorticks = Plt.IntervalsBetween(4), + aspect = 1.67, + ) + + fig200 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.5, P3.thresholds(200.0, 0.5)) for D in D_range], + color = colors[1], + linewidth = lw, + ) + + fig400 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.5, P3.thresholds(400.0, 0.5)) for D in D_range], + color = colors[2], + linewidth = lw, + ) + + fig800 = Plt.lines!( + ax, + D_range * 1e3, + [area(D, 0.5, P3.thresholds(800.0, 0.5)) for D in D_range], + color = colors[3], + linewidth = lw, + ) + + d_thb = Plt.lines!( + ax, + [D = D_th * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = "red", + linewidth = lw, + ) + + d_cr_200 = Plt.lines!( + ax, + [D = P3.thresholds(200.0, 0.5)[1] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = threshold_colors[1], + linewidth = lw, + ) + + d_cr_400 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.5)[1] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = threshold_colors[2], + linewidth = lw, + ) + + d_cr_800 = Plt.lines!( + ax, + [D = P3.thresholds(800.0, 0.5)[1] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "---", + color = threshold_colors[3], + linewidth = lw, + ) + + d_gr_200 = Plt.lines!( + ax, + [D = P3.thresholds(200.0, 0.5)[2] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "...", + color = threshold_colors[1], + linewidth = lw, + ) + + d_gr_400 = Plt.lines!( + ax, + [D = P3.thresholds(400.0, 0.8)[2] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "...", + color = threshold_colors[2], + linewidth = lw, + ) + + d_gr_800 = Plt.lines!( + ax, + [D = P3.thresholds(800.0, 0.5)[2] * 1e3 for D in D_range], + range(1e-9, stop = 1e-5, length = len_D_range), + linestyle = "...", + color = threshold_colors[3], + linewidth = lw, + ) + + leg1_b = Plt.Legend( + fig[11:12, 2], + [fig200, fig400, fig800], + [ + 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(fig[11:12, 3], [d_thb], [Plt.L"$D_{th}$"]) + + leg1_b_dcr = Plt.Legend( + fig[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( + fig[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("P3Scheme_Area_2.svg", fig) +end diff --git a/test/Project.toml b/test/Project.toml index 9adea3252..fdd36a9e0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" +CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 6ca3215eb..b88a416fa 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 # deleted GPU tests for P3 from this branch const liquid = CMT.LiquidType() const ice = CMT.IceType()