Merge pull request #853 from sajjadazimi/condensation_thd_update_option

an option to fix thd in condensation dynamics is added

Author: slayoo

PySDM/dynamics/condensation.py

@@ -27,6 +27,7 @@ def __init__(
         dt_cond_range: tuple = DEFAULTS.cond_range,
         schedule: str = DEFAULTS.schedule,
         max_iters: int = 16,
+        update_thd: bool = True,
     ):
 
         self.particulator = None
@@ -47,6 +48,8 @@ def __init__(
 
         self.cell_order = None
 
+        self.update_thd = update_thd
+
     def register(self, builder):
         self.particulator = builder.particulator
 
@@ -100,6 +103,10 @@ def __call__(self):
             )
             if not self.success.all():
                 raise RuntimeError("Condensation failed")
+            if not self.update_thd:
+                self.particulator.environment.get_predicted("thd").ravel(
+                    self.particulator.environment.get_thd()
+                )
             # note: this makes order of dynamics matter
             #       (e.g., condensation after chemistry or before)
             self.particulator.update_TpRH()

test-time-requirements.txt

@@ -4,5 +4,5 @@ ghapi
 pytest
 
 # note: if cloning both PySDM and PySDM examples, consider "pip install -e"
-PySDM-examples @ git+https://github.com/atmos-cloud-sim-uj/PySDM-examples@050ab12#egg=PySDM-examples
+PySDM-examples @ git+https://github.com/sajjadazimi/PySDM-examples@e8797c6#egg=PySDM-examples
 PyMPDATA @ git+https://github.com/atmos-cloud-sim-uj/PyMPDATA@e7b73a7#egg=PyMPDATA

tests/smoke_tests/shipway_and_hill_2012/test_few_steps.py

@@ -28,6 +28,7 @@ def test_few_steps_no_precip(particle_reservoir_depth, plot=False):
     )
     settings.particle_reservoir_depth = particle_reservoir_depth
     settings.t_max = 50 * settings.dt
+    settings.condensation_update_thd = True
     simulation = Simulation(settings)
 
     # Act
@@ -80,3 +81,38 @@ def mean_profile_over_last_steps(var, smooth=True):
 
     assert max(mean_profile_over_last_steps("ripening rate")) > 0
     assert max(mean_profile_over_last_steps("deactivating rate")) > 0
+
+
+def test_fixed_thd():
+    # Arrange
+    n_sd_per_gridbox = 128
+    settings = Settings(
+        n_sd_per_gridbox=n_sd_per_gridbox,
+        dt=30 * si.s,
+        dz=60 * si.m,
+        precip=False,
+        rho_times_w_1=2 * si.m / si.s * si.kg / si.m**3,
+    )
+    settings.t_max = 50 * settings.dt
+    settings.condensation_update_thd = False
+    simulation = Simulation(settings)
+
+    # Act
+    output = simulation.run().products
+
+    # Assert
+    mean_profile_over_last_steps = lambda var: np.mean(
+        output[var][output["z"] >= 0, -10:], axis=1
+    )
+    sd_prof = mean_profile_over_last_steps("super droplet count per gridbox")
+    assert 0.5 * n_sd_per_gridbox < min(sd_prof) < 1.5 * n_sd_per_gridbox
+    assert 0.5 * n_sd_per_gridbox < max(sd_prof) < 1.5 * n_sd_per_gridbox
+
+    assert 0.01 < max(mean_profile_over_last_steps("peak supersaturation")) < 0.1
+    assert min(mean_profile_over_last_steps("ql")) < 1e-10
+    assert 0.6 < max(mean_profile_over_last_steps("ql")) < 0.9
+    assert max(mean_profile_over_last_steps("activating rate")) == 0
+    assert max(mean_profile_over_last_steps("ripening rate")) == 0
+    assert max(mean_profile_over_last_steps("deactivating rate")) == 0
+
+    assert sum(np.amin(output["thd"], axis=1)) == sum(np.amax(output["thd"], axis=1))