P3 aspect ratio as option, docs on regimes and density, The Duck

CliMA · Aug 12, 2024 · af77c95 · af77c95
1 parent 5b1fd4b
commit af77c95
Show file tree

Hide file tree

Showing 10 changed files with 576 additions and 128 deletions.
diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md
@@ -68,6 +68,12 @@ where
     TODO - Check units, see in [issue #151](https://github.com/CliMA/CloudMicrophysics.jl/issues/151)
 
 Below we show the m(D) and a(D) regimes replicating Figures 1 (a) and (b) from [MorrisonMilbrandt2015](@cite).
+We also show the density as a function of D.
+Note that because graupel is completely filled with rime,
+  the density (``\rho_{g}``) is independent of ``D`` between ``D_{gr}`` and ``D_{cr}``.
+Following [MorrisonMilbrandt2015](@cite), for nonspherical particles
+  ``\rho_{ice}``is assumed to be equal to the mass of the particle
+  divided by the volume of a sphere with the same particle size
 ```@example
 include("plots/P3SchemePlots.jl")
 ```
@@ -115,7 +121,7 @@ N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D = \int_{0}^{\infty} \! N_{0} D^
 
 ``L_{ice}`` depends on the variable mass-size relation ``m(D)`` defined above.
 We solve for ``L_{ice}`` in a piece-wise fashion defined by the same thresholds as ``m(D)``.
-As a result ``L_{ice}`` can be expressed as a sum of inclomplete gamma functions,
+As a result ``L_{ice}`` can be expressed as a sum of incomplete gamma functions,
   and the shape parameters are found using iterative solver.
 
 |      condition(s)                            |    ``L_{ice} = \int \! m(D) N'(D) \mathrm{d}D``                                          |         gamma representation          |
@@ -175,6 +181,14 @@ V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D}
 where ``\phi = (16 \rho_{ice}^2 A(D)^3) / (9 \pi m(D)^2)`` is the aspect ratio,
 and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters.
 
+Note that ``\phi = 1`` corresponds to spherical particles
+  (small spherical ice (``D < D_{th}``) and graupel (``D_{gr} < D < D_{cr}``)).
+The plot provided below helps to visualize the transitions between spherical and nonspherical regimes.
+```@example
+include("plots/P3AspectRatioPlot.jl")
+```
+![](P3Scheme_aspect_ratio.svg)
+
 The mass-weighted fall speed (``V_m``) and the number-weighted fall speed (``V_n``) are calculated as
 ```math
 V_m = \frac{\int_{0}^{\infty} \! V(D) m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D}
@@ -188,9 +202,21 @@ We also plot the mass-weighted mean particle size ``D_m`` which is given by:
 D_m = \frac{\int_{0}^{\infty} \! D m(D) N'(D) \mathrm{d}D}{\int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D}
 ```
 
-Below w show these relationships for small, medium, and large ``D_m``
+Below we provide plots of these relationships for small, medium, and large ``D_m``:
+  the first row highlights the particle size regime,
+  the second displays ``D_m`` of the particles,
+  the third shows the aspect ratio ``\phi (D_m)``,
+  and the final row exhibits ``V_m``.
 They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite).
 ```@example
 include("plots/P3TerminalVelocityPlots.jl")
 ```
 ![](MorrisonandMilbrandtFig2.svg)
+
+## Acknowledgments
+
+Click on the P3 mascot duck to be taken to the repository
+  in which the authors of [MorrisonMilbrandt2015](@cite) and others
+  have implemented the P3 scheme in Fortran!
+
+[![P3 mascot](assets/p3_mascot.png)](https://github.com/P3-microphysics/P3-microphysics)
diff --git a/docs/src/assets/p3_mascot.png b/docs/src/assets/p3_mascot.png
diff --git a/docs/src/plots/P3AspectRatioPlot.jl b/docs/src/plots/P3AspectRatioPlot.jl
@@ -0,0 +1,105 @@
+import CairoMakie as CMK
+import CloudMicrophysics as CM
+import ClimaParams as CP
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.P3Scheme as P3
+
+FT = Float64
+
+const PSP3 = CMP.ParametersP3
+p3 = CMP.ParametersP3(FT)
+
+function define_axis(
+    fig,
+    row_range,
+    col_range,
+    title,
+    ylabel,
+    yticks,
+    aspect;
+    logscale = true,
+)
+    return CMK.Axis(
+        fig[row_range, col_range],
+        title = title,
+        xlabel = CMK.L"$D$ (mm)",
+        ylabel = ylabel,
+        xscale = CMK.log10, #ifelse(logscale, CMK.log10, CMK.identity),
+        yscale = ifelse(logscale, CMK.log10, CMK.identity),
+        xminorticksvisible = true,
+        xminorticks = CMK.IntervalsBetween(5),
+        yminorticksvisible = true,
+        yminorticks = CMK.IntervalsBetween(3),
+        xticks = [0.01, 0.1, 1, 10],
+        yticks = yticks,
+        aspect = aspect,
+        limits = ((0.02, 10.0), nothing),
+    )
+end
+
+#! format: off
+function p3_aspect()
+
+    D_range = range(3e-5, stop = 1e-2, length = Int(1e4))
+    logocolors = CMK.Colors.JULIA_LOGO_COLORS
+    cl = [logocolors.blue, logocolors.green, logocolors.red, logocolors.purple]
+    lw = 3
+
+    fig = CMK.Figure()
+
+    # define plot axis
+    #[row, column]
+    ax1 = define_axis(fig, 1:7,  1:9,   CMK.L"ϕᵢ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.9, logscale = false)
+    ax3 = define_axis(fig, 1:7,  10:18, CMK.L"ϕᵢ(D) regime for $F_r = 0.4$",          CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.8, logscale = false)
+
+    # Get thresholds
+    sol4_0 = P3.thresholds(p3, 400.0, 0.0)
+    sol4_5 = P3.thresholds(p3, 400.0, 0.5)
+    sol4_8 = P3.thresholds(p3, 400.0, 0.8)
+    # ϕᵢ(D)
+    fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.0, sol4_0) for D in D_range], color = cl[1], linewidth = lw)
+    fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.ϕᵢ(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, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
+    for ax in [ax1]
+        global d_tha  = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
+        global d_cr_5 = CMK.vlines!(ax, sol4_5[1]          * 1e3, linestyle = :dot,  color = cl[2], linewidth = lw)
+        global d_cr_8 = CMK.vlines!(ax, sol4_8[1]          * 1e3, linestyle = :dot,  color = cl[3], linewidth = lw)
+        global d_gr_5 = CMK.vlines!(ax, sol4_5[2]          * 1e3, linestyle = :dash, color = cl[2], linewidth = lw)
+        global d_gr_8 = CMK.vlines!(ax, sol4_8[2]          * 1e3, linestyle = :dash, color = cl[3], linewidth = lw)
+    end
+
+    # get thresholds
+    sol_2 = P3.thresholds(p3, 200.0, 0.4)
+    sol_4 = P3.thresholds(p3, 400.0, 0.4)
+    sol_8 = P3.thresholds(p3, 800.0, 0.4)
+    # ϕᵢ(D)
+    fig2_200 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    fig2_400 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    fig2_800 = CMK.lines!(ax3, D_range * 1e3, [P3.ϕᵢ(p3, D, 0.4, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    # plot verical lines
+    for ax in [ax3]
+        global d_thb    = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
+        global d_cr_200 = CMK.vlines!(ax, sol_2[1] * 1e3,           linestyle = :dot,  color = cl[1], linewidth = lw)
+        global d_cr_400 = CMK.vlines!(ax, sol_4[1] * 1e3,           linestyle = :dot,  color = cl[2], linewidth = lw)
+        global d_cr_800 = CMK.vlines!(ax, sol_8[1] * 1e3,           linestyle = :dot,  color = cl[3], linewidth = lw)
+        global d_gr_200 = CMK.vlines!(ax, sol_2[2] * 1e3,           linestyle = :dash, color = cl[1], linewidth = lw)
+        global d_gr_400 = CMK.vlines!(ax, sol_4[2] * 1e3,           linestyle = :dash, color = cl[2], linewidth = lw)
+        global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3,           linestyle = :dash, color = cl[3], linewidth = lw)
+    end
+    # add legend
+    CMK.Legend(fig[8:9, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false)
+    CMK.Legend(fig[8:9, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false)
+    CMK.Legend(fig[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)
+    CMK.Legend(fig[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)
+
+    CMK.Legend(fig[8:9, 13], [fig2_200, fig2_400, fig2_800],    [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)
+    CMK.Legend(fig[8:9, 14], [d_thb],                           [CMK.L"$D_{th}$"], framevisible = false)
+    CMK.Legend(fig[8:9, 17], [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)
+    CMK.Legend(fig[8:9, 16], [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)
+
+    CMK.resize_to_layout!(fig)
+    CMK.save("P3Scheme_aspect_ratio.svg", fig)
+end
+#! format: on
+
+p3_aspect()
diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl
@@ -9,14 +9,23 @@ FT = Float64
 const PSP3 = CMP.ParametersP3
 p3 = CMP.ParametersP3(FT)
 
-function define_axis(fig, row_range, col_range, title, ylabel, yticks, aspect)
+function define_axis(
+    fig,
+    row_range,
+    col_range,
+    title,
+    ylabel,
+    yticks,
+    aspect;
+    logscale = true,
+)
     return CMK.Axis(
         fig[row_range, col_range],
         title = title,
         xlabel = CMK.L"$D$ (mm)",
         ylabel = ylabel,
-        xscale = CMK.log10,
-        yscale = CMK.log10,
+        xscale = CMK.log10, #ifelse(logscale, CMK.log10, CMK.identity),
+        yscale = ifelse(logscale, CMK.log10, CMK.identity),
         xminorticksvisible = true,
         xminorticks = CMK.IntervalsBetween(5),
         yminorticksvisible = true,
@@ -44,8 +53,11 @@ function p3_relations_plot()
     ax3 = define_axis(fig, 1:7,  10:18, CMK.L"m(D) regime for $F_r = 0.95$",          CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8)
     ax2 = define_axis(fig, 8:15, 1:9,   CMK.L"A(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.7)
     ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_r = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6)
+    ax5 = define_axis(fig, 16:22, 1:9,   CMK.L"ρ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 2, logscale = false)
+    ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_r = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false)
 
     # Get thresholds
+    sol4_0 = P3.thresholds(p3, 400.0, 0.0)
     sol4_5 = P3.thresholds(p3, 400.0, 0.5)
     sol4_8 = P3.thresholds(p3, 400.0, 0.8)
     # m(D)
@@ -56,8 +68,12 @@ function p3_relations_plot()
     fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0        ) for D in D_range], color = cl[1], linewidth = lw)
     fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
     fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
+    # ρ(D)
+    density_0 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.0, sol4_0) for D in D_range], color = cl[1], linewidth = lw)
+    density_5 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw)
+    density_8 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw)
     # plot verical lines
-    for ax in [ax1, ax2]
+    for ax in [ax1, ax2, ax5]
         global d_tha  = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
         global d_cr_5 = CMK.vlines!(ax, sol4_5[1]          * 1e3, linestyle = :dot,  color = cl[2], linewidth = lw)
         global d_cr_8 = CMK.vlines!(ax, sol4_8[1]          * 1e3, linestyle = :dot,  color = cl[3], linewidth = lw)
@@ -77,8 +93,12 @@ function p3_relations_plot()
     fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_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, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
     fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
+    # ρ(D)
+    density_200 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw)
+    density_400 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw)
+    density_800 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw)
     # plot verical lines
-    for ax in [ax3, ax4]
+    for ax in [ax3, ax4, ax6]
         global d_thb    = CMK.vlines!(ax, P3.D_th_helper(p3) * 1e3, linestyle = :dash, color = cl[4], linewidth = lw)
         global d_cr_200 = CMK.vlines!(ax, sol_2[1] * 1e3,           linestyle = :dot,  color = cl[1], linewidth = lw)
         global d_cr_400 = CMK.vlines!(ax, sol_4[1] * 1e3,           linestyle = :dot,  color = cl[2], linewidth = lw)
@@ -88,15 +108,15 @@ function p3_relations_plot()
         global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3,           linestyle = :dash, color = cl[3], linewidth = lw)
     end
     # add legend
-    CMK.Legend(fig[16:17, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false)
-    CMK.Legend(fig[16:17, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false)
-    CMK.Legend(fig[16:17, 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)
-    CMK.Legend(fig[16:17, 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)
+    CMK.Legend(fig[24:25, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false)
+    CMK.Legend(fig[24:25, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false)
+    CMK.Legend(fig[24:25, 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)
+    CMK.Legend(fig[24:25, 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)
 
-    CMK.Legend(fig[16:17, 13], [fig3_200, fig3_400, fig3_800],    [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)
-    CMK.Legend(fig[16:17, 14], [d_thb],                           [CMK.L"$D_{th}$"], framevisible = false)
-    CMK.Legend(fig[16:17, 17], [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)
-    CMK.Legend(fig[16:17, 16], [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)
+    CMK.Legend(fig[24:25, 13], [fig3_200, fig3_400, fig3_800],    [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)
+    CMK.Legend(fig[24:25, 14], [d_thb],                           [CMK.L"$D_{th}$"], framevisible = false)
+    CMK.Legend(fig[24:25, 17], [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)
+    CMK.Legend(fig[24:25, 16], [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)
 
     CMK.resize_to_layout!(fig)
     CMK.save("P3Scheme_relations.svg", fig)