minor changes in docs, fix perf test

CliMA · Aug 14, 2023 · 09a7242 · 09a7242
1 parent 11b790f
commit 09a7242
Show file tree

Hide file tree

Showing 3 changed files with 78 additions and 46 deletions.
diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md
@@ -11,28 +11,38 @@ Traditionally, cloud ice microphysics schemes use various predefined
  a single ice category and evolving its properties. This simplification
  aids in attempts to constrain the scheme's free parameters.
 
+The prognostic variables are:
+ - ``N_{ice}`` - number concentration 1/m3
+ - ``q_{ice}`` - ice mass density kg/m3
+ - ``q_{rim}`` - rime mass density kg/m3
+ - ``B_{rim}`` - rime volume - TODO - units?
+
+!!! note
+    TODO - At some point we should switch to specific humidities...
+
 ## 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 
+ 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),
+ - ``\mu`` is the shape parameter (dimensionless),
  - ``\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
+We assume ``\mu \ = 0.00191 \; \lambda \ ^{0.8} - 2``.
+Following [MorrisonGrabowski2008](@cite) we limit ``\mu \ \in (0,6)``.
+A non-zero ``\mu`` can occur only for very small mean particle sizes``\frac{1}{\lambda} < ~0.17 mm``.
+``N_0`` and ``\lambda`` can be found using different moments of the PSD,
+namely the total number concentration ``N`` and mass mixing ratio ``q``, where
 
 ```math
 N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D
@@ -42,44 +52,47 @@ N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D
 q = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D
 ```
 
+!!! note
+    TODO - write the formulae for ``N_0`` and ``\lambda``
+
 
 ## 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
+ is a piecewise function with variable thresholds 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:
+| 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}}``             |
+|graupel (completely rimed, spherical) | ``q_{rim} > 0`` and ``D_{gr} < D < D_{cr}``  | ``\frac{\pi}{6} \rho_g \ D^3``               |
+|partially rimed ice                   | ``q_{rim} > 0`` and ``D > D_{cr}``           | ``\frac{\alpha_{va}}{1-F_r} D^{\beta_{va}}`` |
 
-  - ``\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):
+where:
+ - ``D_{th}``, ``D_{gr}``, ``D_{cr}`` are particle size thresholds in ``m``,
+ - ``\rho_i`` is cloud ice density in ``kg m^{-3}``,
+ - ``\beta_{va} = 1.9`` is a dimensionless parameter from [BrownFrancis1995](@cite) (based on measurements of vapor diffusion and aggregation in midlatitude cirrus),
+ - ``\alpha_{va} = 7.38 \; 10^{-11} \; 10^{6 \beta_{va} - 3}`` in ``kg m^{-β_{va}}`` is a parameter from [BrownFrancis1995](@cite) in base SI units,
+ - ``\rho_g`` is the bulk density of graupel in ``kg m^{-3}``.
+
+The first threshold is solely determined by the free parameters:
+  ``D_{th} = (\frac{\pi \rho_i}{6\alpha_{va}})^{\frac{1}{\beta_{va} - 3}}``.
+The remaining thresholds: ``D_{gr}``, ``D_{cr}``, as well as the
+  bulk density of graupel ``\rho_{g}``,
+  and the bulk density of the unrimed part ``\rho_d``
+  form a nonlinear system:
+ - ``D_{gr} = (\frac{6\alpha_{va}}{\pi \rho_g})^{\frac{1}{3 - \beta_{va}}}``
+ - ``D_{cr} = [ (\frac{1}{1-F_r}) \frac{6 \alpha_{va}}{\pi \rho_g} ]^{\frac{1}{3 - \beta_{va}}}``
+ - ``\rho_g = \rho_r F_r + (1 - F_r) \rho_d``
+ - ``\rho_d = \frac{6\alpha_{va}(D_{cr}^{\beta{va} \ - 2} - D_{gr}^{\beta{va} \ - 2})}{\pi \ (\beta_{va} \ - 2)(D_{cr}-D_{gr})}``
+where
+ - ``F_r = \frac{q_{rim}}{q_{ice}}`` is the rime mass fraction,
+ - ``\rho_{r} = \frac{q_{rim}}{B_{rim}}`` is the predicted rime density.
+The system is solved using [`NonlinearSolve.jl`](https://docs.sciml.ai/NonlinearSolve/stable/).
+
+Below we show the m(D) regime, replicating Figures 1 (a) and (b) from [MorrisonMilbrandt2015](@cite).
 
 ```@example
 include("P3SchemePlots.jl")
@@ -89,4 +102,4 @@ p3_m_plot2(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4"
 
 ![](MorrisonandMilbrandtFig1a.svg)
 
-![](MorrisonandMilbrandtFig1b.svg)
+![](MorrisonandMilbrandtFig1b.svg)
diff --git a/docs/src/P3SchemePlots.jl b/docs/src/P3SchemePlots.jl
@@ -13,6 +13,7 @@ thermo_params = CMP.thermodynamics_params(param_set)
 
 # bulk density of ice
 const ρ_i::FT = CMP.ρ_cloud_ice(param_set)
+# TODO - use from CLIMAParameters
 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))
@@ -106,6 +107,8 @@ function p3_m_plot1(
 )
     D_range = range(1e-5, stop = 1e-2, length = len_D_range)
 
+    lw = 4
+
     fig1_a = Plt.Figure()
 
     ax1_a = Plt.Axis(
@@ -128,20 +131,23 @@ function p3_m_plot1(
         D_range * 1e3,
         [mass(D, 0.0) for D in D_range],
         color = colors[1],
+        linewidth = lw,
     )
 
     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],
+        linewidth = lw,
     )
 
     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],
+        linewidth = lw,
     )
 
     d_cr_5 = Plt.lines!(
@@ -150,6 +156,7 @@ function p3_m_plot1(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = threshold_colors[2],
+        linewidth = lw,
     )
 
     d_cr_8 = Plt.lines!(
@@ -158,6 +165,7 @@ function p3_m_plot1(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = threshold_colors[3],
+        linewidth = lw,
     )
 
     d_gr_5 = Plt.lines!(
@@ -166,7 +174,7 @@ function p3_m_plot1(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "...",
         color = threshold_colors[2],
-        linewidth = 1.5,
+        linewidth = lw,
     )
 
     d_gr_8 = Plt.lines!(
@@ -175,7 +183,7 @@ function p3_m_plot1(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "...",
         color = threshold_colors[3],
-        linewidth = 1.5,
+        linewidth = lw,
     )
 
     d_tha = Plt.lines!(
@@ -184,6 +192,7 @@ function p3_m_plot1(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = "red",
+        linewidth = lw,
     )
 
     leg1_a = Plt.Legend(
@@ -227,6 +236,8 @@ function p3_m_plot2(
 
     fig1_b = Plt.Figure()
 
+    lw = 4
+
     ax1_b = Plt.Axis(
         fig1_b[1:10, 1:11],
         title = Plt.L"m(D) regime for $F_r = 0.95$",
@@ -247,20 +258,23 @@ function p3_m_plot2(
         D_range * 1e3,
         [mass(D, 0.95, P3.thresholds(200.0, 0.95)) for D in D_range],
         color = colors[1],
+        linewidth = lw,
     )
 
     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],
+        linewidth = lw,
     )
 
     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],
+        linewidth = lw,
     )
 
     d_thb = Plt.lines!(
@@ -269,6 +283,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = "red",
+        linewidth = lw,
     )
 
     d_cr_200 = Plt.lines!(
@@ -277,6 +292,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = threshold_colors[1],
+        linewidth = lw,
     )
 
     d_cr_400 = Plt.lines!(
@@ -285,6 +301,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = threshold_colors[2],
+        linewidth = lw,
     )
 
     d_cr_800 = Plt.lines!(
@@ -293,6 +310,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "---",
         color = threshold_colors[3],
+        linewidth = lw,
     )
 
     d_gr_200 = Plt.lines!(
@@ -301,7 +319,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "...",
         color = threshold_colors[1],
-        linewidth = 1.5,
+        linewidth = lw,
     )
 
     d_gr_400 = Plt.lines!(
@@ -310,6 +328,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "...",
         color = threshold_colors[2],
+        linewidth = lw,
     )
 
     d_gr_800 = Plt.lines!(
@@ -318,7 +337,7 @@ function p3_m_plot2(
         range(1e-14, stop = 1e-4, length = len_D_range),
         linestyle = "...",
         color = threshold_colors[3],
-        linewidth = 1.5,
+        linewidth = lw,
     )
 
     leg1_b = Plt.Legend(

diff --git a/test/performance_tests.jl b/test/performance_tests.jl
@@ -22,7 +22,7 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl"))
 
 @info "Performance Tests"
 
-function bench_press(foo, args, min_run_time)
+function bench_press(foo, args, min_run_time, min_memory=0.0, min_allocs=0.0)
 
     println("Testing ", "$foo")
     # Calling foo once before benchmarking
@@ -34,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 <= 4e6
-    TT.@test trail.allocs <= 4e4
+    TT.@test trail.memory <= min_memory
+    TT.@test trail.allocs <= min_allocs
 end
 
 function benchmark_test(FT)
@@ -90,7 +90,7 @@ function benchmark_test(FT)
     Delta_a_w = FT(0.27)
 
     # P3 scheme
-    bench_press(P3.thresholds, (ρ_r, F_r), 12e6)
+    bench_press(P3.thresholds, (ρ_r, F_r), 12e6, 4e6, 4e4)
 
     # aerosol activation
     bench_press(