updates to the Grabowski & Pawlowska 2023 example

examples/PySDM_examples/Grabowski_and_Pawlowska_2023/settings.py

@@ -2,9 +2,12 @@
 from pystrict import strict
 
 from PySDM import Formulae
+from PySDM.dynamics import condensation
 from PySDM.initialisation.spectra import Lognormal, Sum
 from PySDM.physics import si
 
+condensation_tolerance = condensation.DEFAULTS.rtol_thd / 100
+
 
 @strict
 class Settings:
@@ -20,6 +23,8 @@ def __init__(
         initial_relative_humidity: float = 0.97,
         displacement: float = 1000 * si.m,
         mass_accommodation_coefficient: float = 0.3,
+        rtol_thd: float = condensation_tolerance,
+        rtol_x: float = condensation_tolerance
     ):
         self.formulae = Formulae(constants={"MAC": mass_accommodation_coefficient})
         self.n_sd = n_sd
@@ -64,6 +69,8 @@ def __init__(
         self.t_max = displacement / vertical_velocity
         self.timestep = dt
         self.output_interval = self.timestep
+        self.rtol_thd = rtol_thd
+        self.rtol_x = rtol_x
 
     @property
     def initial_air_density(self):

examples/PySDM_examples/Grabowski_and_Pawlowska_2023/simulation.py

@@ -30,7 +30,9 @@ def __init__(
         builder = Builder(n_sd=settings.n_sd, backend=CPU(formulae=settings.formulae))
         builder.set_environment(env)
         builder.add_dynamic(AmbientThermodynamics())
-        builder.add_dynamic(Condensation())
+        builder.add_dynamic(
+            Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x)
+        )
         for attribute in (
             "critical supersaturation",
             "equilibrium supersaturation",

tests/smoke_tests/parcel/grabowski_and_pawlowska_2023/test_condensation_tolerance.py

@@ -0,0 +1,62 @@
+"""
+test for condensation tolerance parameters: rtol_thd and rtol_x,
+checking if supersaturation has more than one local maximum
+"""
+
+import numpy as np
+import pytest
+from PySDM_examples.Grabowski_and_Pawlowska_2023 import Settings, Simulation
+from PySDM_examples.Grabowski_and_Pawlowska_2023.settings import condensation_tolerance
+from scipy import signal
+
+from PySDM.dynamics import condensation
+from PySDM.physics import si
+from PySDM.products import AmbientRelativeHumidity
+
+PRODUCTS = [
+    AmbientRelativeHumidity(name="S_max", var="RH"),
+]
+
+N_SD = 25
+DZ = 10 * si.m
+VELOCITIES_CM_PER_S = (25, 100, 400)
+AEROSOLS = ("polluted",)
+
+
+@pytest.mark.parametrize(
+    "rtol_cond",
+    (
+        condensation_tolerance,
+        pytest.param(
+            condensation.DEFAULTS.rtol_thd, marks=pytest.mark.xfail(strict=True)
+        ),
+    ),
+)
+@pytest.mark.parametrize("aerosol", AEROSOLS)
+@pytest.mark.parametrize("w_cm_per_s", VELOCITIES_CM_PER_S)
+def test_condensation_tolerance(
+    rtol_cond: float,
+    w_cm_per_s: int,
+    aerosol: str,
+):
+    # arrange
+    output = Simulation(
+        Settings(
+            n_sd=25,
+            dt=10 * si.s,
+            vertical_velocity=w_cm_per_s * si.cm / si.s,
+            aerosol=aerosol,
+            rtol_thd=rtol_cond,
+            rtol_x=rtol_cond,
+        ),
+        products=PRODUCTS,
+    ).run()
+
+    assert (
+        signal.find_peaks(
+            np.asarray(output["products"]["S_max"]),
+            height=(1, None),
+            prominence=(1e-5, None),
+        )[0].shape[0]
+        == 1
+    )

tests/smoke_tests/parcel/grabowski_and_pawlowska_2023/test_ripening_rate.py

@@ -0,0 +1,73 @@
+"""
+test for values of the ripening rate
+"""
+
+import numpy as np
+import pytest
+from PySDM_examples.Grabowski_and_Pawlowska_2023 import Settings, Simulation
+
+from PySDM.physics import si
+from PySDM.products import RipeningRate
+
+PRODUCTS = [
+    RipeningRate(
+        name="ripening",
+    )
+]
+
+N_SD = 200
+
+VELOCITIES_CM_PER_S = (25, 100)
+
+AEROSOLS = ("pristine", "polluted")
+
+DT = 1 * si.s
+
+RTOL = 0.1
+
+
+@pytest.fixture(scope="session", name="outputs")
+def outputs_fixture():
+    outputs = {}
+    for aerosol in AEROSOLS:
+        outputs[aerosol] = {}
+        for w_cm_per_s in VELOCITIES_CM_PER_S:
+            vertical_velocity = w_cm_per_s * si.cm / si.s
+            outputs[aerosol][w_cm_per_s] = Simulation(
+                Settings(
+                    n_sd=N_SD,
+                    dt=DT,
+                    vertical_velocity=vertical_velocity,
+                    aerosol=aerosol,
+                ),
+                products=PRODUCTS,
+            ).run()
+    return outputs
+
+
+@pytest.mark.parametrize("aerosol", AEROSOLS)
+@pytest.mark.parametrize("w_cm_per_s", VELOCITIES_CM_PER_S)
+def test_ripening_rate(
+    outputs: dict,
+    w_cm_per_s: int,
+    aerosol: str,
+):
+    # arrange
+    expected = {
+        "pristine": {
+            25: 0,
+            100: 0,
+        },
+        "polluted": {
+            25: 2.8 * 1e8,
+            100: 4 * 1e8,
+        },
+    }[aerosol][w_cm_per_s]
+    output = outputs[aerosol][w_cm_per_s]["products"]["ripening"]
+
+    # assert
+    assert np.isclose(
+        np.max(output),
+        expected,
+        rtol=RTOL,
+    ).all()