new example: Hill et al, 2023

docs/bibliography.json

@@ -946,5 +946,13 @@
     ],
     "title": "A physically constrained classical description of the homogeneous nucleation of ice in water",
     "label": "Koop and Murray 2016 (J. Chem. Phys. 145)"
+  },
+  "https://doi.org/10.1175/JAS-D-21-0275.1": {
+    "usages": [
+      "examples/PySDM_examples/Hill_et_al_2023/__init__.py",
+      "examples/PySDM_examples/Hill_et_al_2023/figs_1_2_4.ipynb"
+    ],
+    "title": "Toward a Numerical Benchmark for Warm Rain Processes",
+    "label": "Hill et al. 2023 (JAS 80)"
   }
 }

examples/PySDM_examples/Hill_et_al_2023/__init__.py

@@ -0,0 +1,8 @@
+# pylint: disable=invalid-name
+"""
+single-column prescribed-flow constant-temperature example from
+[Hill et al., 2023](https://doi.org/10.1175/JAS-D-21-0275.1)
+
+figs_1_2_4.ipynb
+.. include:: ./figs_1_2_4.ipynb.badges.md
+"""

examples/PySDM_examples/Hill_et_al_2023/cross_sections.py

@@ -0,0 +1,70 @@
+from PySDM.physics import si
+import numpy as np
+
+
+def compute_dvol_andsigma_time_series(output, settings, cloud_base):
+    num_conc = (output["nc"] + output["nr"])[cloud_base]
+    rho_w = settings.formulae.constants.rho_w
+    D_vol = (6 * output["LWC"][cloud_base] / (np.pi * num_conc * rho_w)) ** (1 / 3)
+
+    r_bin_edges = settings.r_bins_edges
+    r_bin_centers = (r_bin_edges[:-1] + r_bin_edges[1:]) / 2
+    r_bin_volumes = settings.formulae.trivia.volume(r_bin_centers)
+
+    nc = output["wet spectrum"][cloud_base] + output["dry spectrum"][cloud_base]
+    mass_bins = (nc * r_bin_volumes[:, np.newaxis]) * rho_w
+    Dk = r_bin_centers * 2
+    Dvolmass = np.sum(mass_bins * Dk[:, np.newaxis], axis=0) / np.sum(mass_bins, axis=0)
+    Nd = np.sum(nc)
+    sigma = np.sqrt(
+        (1 / (Nd - 1)) * np.sum(nc * (Dk[:, np.newaxis] - Dvolmass) ** 2, axis=0)
+    )
+    return D_vol, sigma
+
+
+def compute_dvol_andsigma_vertical(output, settings, zslice, time):
+    t_idx = np.where(output["t"] == time)[0][0]
+    num_conc = (output["nc"] + output["nr"])[zslice, t_idx]
+    rho_w = settings.formulae.constants.rho_w
+
+    D_vol = (6 * output["LWC"][zslice, t_idx] / (np.pi * num_conc * rho_w)) ** (1 / 3)
+
+    r_bin_edges = settings.r_bins_edges
+    r_bin_centers = (r_bin_edges[:-1] + r_bin_edges[1:]) / 2
+    r_bin_volumes = settings.formulae.trivia.volume(r_bin_centers)
+    nc = (
+        output["wet spectrum"][zslice, :, t_idx]
+        + output["dry spectrum"][zslice, :, t_idx]
+    )
+    mass_bins = (nc * r_bin_volumes) * rho_w
+    Dk = r_bin_centers * 2
+    Dvolmass = np.sum(mass_bins * Dk, axis=1) / np.sum(mass_bins, axis=1)
+    Nd = np.sum(nc, axis=1)
+    sigma = np.sqrt(
+        (1 / (Nd - 1)) * np.sum(nc * (Dk - Dvolmass[:, np.newaxis]) ** 2, axis=1)
+    )
+    return D_vol, sigma
+
+
+def compute_LWP_and_nc_time_series(output, settings, z_slice):
+    # LWP = np.mean(output["LWC"][z_slice], axis=0) * settings.z_max
+    LWP = np.sum(output["LWC"][z_slice], axis=0) * settings.dz
+    num_conc = (output["nc"] + output["nr"])[z_slice]
+    mask = output["LWC"][z_slice] > 1e-5
+    masked = num_conc * mask
+    masked[masked == 0] = np.nan
+    mean_num_conc = np.nanmean(masked, axis=0)
+    mean_num_conc[mean_num_conc == np.nan] = 0
+    return LWP, mean_num_conc
+
+
+def compute_LWP_and_nc_vertical(output, z_slice, time):
+    t_idx = np.where(output["t"] == time)[0][0]
+    LWC = output["LWC"][z_slice, t_idx]
+    num_conc = (output["nc"] + output["nr"] + output["na"])[z_slice, t_idx]
+    mask = output["LWC"][z_slice, t_idx] > 1e-5
+    masked = num_conc * mask
+    masked[masked == 0] = np.nan
+    mean_num_conc = masked
+    mean_num_conc[mean_num_conc == np.nan] = 0
+    return LWC, num_conc

examples/PySDM_examples/Shipway_and_Hill_2012/settings.py

@@ -50,6 +50,9 @@ def __init__(
         save_spec_and_attr_times=(),
         collision_kernel=None,
         old_buggy_density_formula=False,
+        cloud_water_radius_range=(1 * si.um, 50 * si.um),
+        rain_water_radius_range=(50 * si.um, np.inf),
+        ignore_moisture_profile_in_density_calc=False,
     ):
         self.formulae = formulae or Formulae()
         self.n_sd_per_gridbox = n_sd_per_gridbox
@@ -112,9 +115,11 @@ def drhod_dz(z_above_reservoir, rhod):
             water_vapour_mixing_ratio = self.water_vapour_mixing_ratio(
                 z_above_reservoir
             )
-            d_water_vapour_mixing_ratio__dz = Derivative(
-                self.water_vapour_mixing_ratio
-            )(z_above_reservoir)
+            d_water_vapour_mixing_ratio__dz = (
+                0
+                if ignore_moisture_profile_in_density_calc
+                else Derivative(self.water_vapour_mixing_ratio)(z_above_reservoir)
+            )
             T = self.formulae.state_variable_triplet.T(
                 rhod[0], self.thd(z_above_reservoir)
             )
@@ -164,8 +169,8 @@ def drhod_dz(z_above_reservoir, rhod):
             self.number_of_bins + 1,
             endpoint=True,
         )
-        self.cloud_water_radius_range = [1 * si.um, 50 * si.um]
-        self.rain_water_radius_range = [50 * si.um, np.inf * si.um]
+        self.cloud_water_radius_range = cloud_water_radius_range
+        self.rain_water_radius_range = rain_water_radius_range
         self.save_spec_and_attr_times = save_spec_and_attr_times
 
     @property

examples/PySDM_examples/Shipway_and_Hill_2012/simulation.py

@@ -141,6 +141,7 @@ def zZ_to_z_above_reservoir(zZ):
             PySDM_products.AmbientWaterVapourMixingRatio(
                 name="water_vapour_mixing_ratio"
             ),
+            PySDM_products.LiquidWaterContent(name="LWC"),
         ]
         if settings.enable_condensation:
             self.products.extend(