introducing ice thawing to Freezing dynamic (still disabled by default) (#993)

Author: slayoo

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -18,22 +18,37 @@ def __init__(self):
         BackendMethods.__init__(self)
         const = self.formulae.constants
         unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated
+        frozen_and_above_freezing_point = (
+            self.formulae.trivia.frozen_and_above_freezing_point
+        )
 
         @numba.njit(
             **{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath, "parallel": False}
         )
         def _freeze(volume, i):
             volume[i] = -1 * volume[i] * const.rho_w / const.rho_i
             # TODO #599: change thd (latent heat)!
-            # TODO #599: handle the negative volume in tests, attributes, products, dynamics, ...
+
+        @numba.njit(
+            **{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath, "parallel": False}
+        )
+        def _thaw(volume, i):
+            volume[i] = -1 * volume[i] * const.rho_i / const.rho_w
+            # TODO #599: change thd (latent heat)!
 
         @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath})
-        def freeze_singular_body(attributes, temperature, relative_humidity, cell):
+        def freeze_singular_body(
+            attributes, temperature, relative_humidity, cell, thaw
+        ):
             n_sd = len(attributes.freezing_temperature)
             for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 if attributes.freezing_temperature[i] == 0:
                     continue
-                if (
+                if thaw and frozen_and_above_freezing_point(
+                    attributes.wet_volume[i], temperature[cell[i]]
+                ):
+                    _thaw(attributes.wet_volume, i)
+                elif (
                     unfrozen_and_saturated(
                         attributes.wet_volume[i], relative_humidity[cell[i]]
                     )
@@ -56,13 +71,18 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
             relative_humidity,
             record_freezing_temperature,
             freezing_temperature,
+            thaw,
         ):
             n_sd = len(attributes.wet_volume)
             for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 if attributes.immersed_surface_area[i] == 0:
                     continue
                 cell_id = cell[i]
-                if unfrozen_and_saturated(
+                if thaw and frozen_and_above_freezing_point(
+                    attributes.wet_volume[i], temperature[cell_id]
+                ):
+                    _thaw(attributes.wet_volume, i)
+                elif unfrozen_and_saturated(
                     attributes.wet_volume[i], relative_humidity[cell_id]
                 ):
                     rate = j_het(a_w_ice[cell_id])
@@ -77,7 +97,9 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
 
         self.freeze_time_dependent_body = freeze_time_dependent_body
 
-    def freeze_singular(self, *, attributes, temperature, relative_humidity, cell):
+    def freeze_singular(
+        self, *, attributes, temperature, relative_humidity, cell, thaw: bool
+    ):
         self.freeze_singular_body(
             SingularAttributes(
                 freezing_temperature=attributes.freezing_temperature.data,
@@ -86,6 +108,7 @@ def freeze_singular(self, *, attributes, temperature, relative_humidity, cell):
             temperature.data,
             relative_humidity.data,
             cell.data,
+            thaw=thaw,
         )
 
     def freeze_time_dependent(
@@ -99,7 +122,8 @@ def freeze_time_dependent(
         temperature,
         relative_humidity,
         record_freezing_temperature,
-        freezing_temperature
+        freezing_temperature,
+        thaw: bool
     ):
         self.freeze_time_dependent_body(
             rand.data,
@@ -110,10 +134,11 @@ def freeze_time_dependent(
             timestep,
             cell.data,
             a_w_ice.data,
-            temperature.data if record_freezing_temperature else None,
+            temperature.data,
             relative_humidity.data,
             record_freezing_temperature=record_freezing_temperature,
-            freezing_temperature=freezing_temperature.data
-            if record_freezing_temperature
-            else None,
+            freezing_temperature=(
+                freezing_temperature.data if record_freezing_temperature else None
+            ),
+            thaw=thaw,
         )

PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py

@@ -23,13 +23,20 @@ def __init__(self):
                 "cell",
                 "a_w_ice",
                 "relative_humidity",
+                "thaw",
+                "temperature",
             ),
             name_iter="i",
             body=f"""
                 if (immersed_surface_area[i] == 0) {{
                     return;
                 }}
-                if ({self.formulae.trivia.unfrozen_and_saturated.c_inline(
+                if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline(
+                    volume="wet_volume[i]",
+                    temperature="temperature[cell[i]]"
+                )}) {{
+                    wet_volume[i] = -1 * wet_volume[i] * {const.rho_i} / {const.rho_w};
+                }} else if ({self.formulae.trivia.unfrozen_and_saturated.c_inline(
                         volume="wet_volume[i]",
                         relative_humidity="relative_humidity[cell[i]]"
                     )}) {{
@@ -55,13 +62,19 @@ def __init__(self):
                 "temperature",
                 "relative_humidity",
                 "cell",
+                "thaw",
             ),
             name_iter="i",
             body=f"""
                 if (freezing_temperature[i] == 0) {{
                     return;
                 }}
-                if (
+                if (thaw && {self.formulae.trivia.frozen_and_above_freezing_point.c_inline(
+                    volume="wet_volume[i]",
+                    temperature="temperature[cell[i]]"
+                )}) {{
+                    wet_volume[i] = -1 * wet_volume[i] * {const.rho_i} / {const.rho_w};
+                }} else if (
                     {self.formulae.trivia.unfrozen_and_saturated.c_inline(
                         volume="wet_volume[i]",
                         relative_humidity="relative_humidity[cell[i]]"
@@ -75,7 +88,9 @@ def __init__(self):
         )
 
     @nice_thrust(**NICE_THRUST_FLAGS)
-    def freeze_singular(self, *, attributes, temperature, relative_humidity, cell):
+    def freeze_singular(
+        self, *, attributes, temperature, relative_humidity, cell, thaw
+    ):
         n_sd = len(attributes.freezing_temperature)
         self.freeze_singular_body.launch_n(
             n=n_sd,
@@ -85,6 +100,7 @@ def freeze_singular(self, *, attributes, temperature, relative_humidity, cell):
                 temperature.data,
                 relative_humidity.data,
                 cell.data,
+                trtc.DVBool(thaw),
             ),
         )
 
@@ -101,6 +117,7 @@ def freeze_time_dependent(  # pylint: disable=unused-argument
         relative_humidity,
         record_freezing_temperature,
         freezing_temperature,
+        thaw,
     ):
         n_sd = len(attributes.immersed_surface_area)
         self.freeze_time_dependent_body.launch_n(
@@ -113,5 +130,7 @@ def freeze_time_dependent(  # pylint: disable=unused-argument
                 cell.data,
                 a_w_ice.data,
                 relative_humidity.data,
+                trtc.DVBool(thaw),
+                temperature.data,
             ),
         )

PySDM/dynamics/freezing.py

@@ -9,10 +9,11 @@
 
 
 class Freezing:
-    def __init__(self, *, singular=True, record_freezing_temperature=False):
+    def __init__(self, *, singular=True, record_freezing_temperature=False, thaw=False):
         assert not (record_freezing_temperature and singular)
         self.singular = singular
         self.record_freezing_temperature = record_freezing_temperature
+        self.thaw = thaw
         self.enable = True
         self.rand = None
         self.rng = None
@@ -58,6 +59,7 @@ def __call__(self):
                 temperature=self.particulator.environment["T"],
                 relative_humidity=self.particulator.environment["RH"],
                 cell=self.particulator.attributes["cell id"],
+                thaw=self.thaw,
             )
         else:
             self.rand.urand(self.rng)
@@ -72,18 +74,15 @@ def __call__(self):
                 timestep=self.particulator.dt,
                 cell=self.particulator.attributes["cell id"],
                 a_w_ice=self.particulator.environment["a_w_ice"],
-                temperature=(
-                    self.particulator.environment["T"]
-                    if self.record_freezing_temperature
-                    else None
-                ),
+                temperature=self.particulator.environment["T"],
                 relative_humidity=self.particulator.environment["RH"],
                 record_freezing_temperature=self.record_freezing_temperature,
                 freezing_temperature=(
                     self.particulator.attributes["freezing temperature"]
                     if self.record_freezing_temperature
                     else None
                 ),
+                thaw=self.thaw,
             )
 
         self.particulator.attributes.mark_updated("volume")

PySDM/physics/heterogeneous_ice_nucleation_rate/null.py

@@ -10,4 +10,4 @@ def __init__(self, _):
 
     @staticmethod
     def j_het(const, a_w_ice):
-        return const.NaN
+        return 0

PySDM/physics/trivia.py

@@ -75,3 +75,7 @@ def th_std(const, p, T):
     @staticmethod
     def unfrozen_and_saturated(_, volume, relative_humidity):
         return volume > 0 and relative_humidity > 1
+
+    @staticmethod
+    def frozen_and_above_freezing_point(const, volume, temperature):
+        return volume < 0 and temperature > const.T0

test-time-requirements.txt

@@ -6,7 +6,7 @@ ghapi
 pytest
 
 # note: if cloning both PySDM and PySDM examples, consider "pip install -e"
-PySDM-examples @ git+https://github.com/open-atmos/PySDM-examples@151a869#egg=PySDM-examples
+PySDM-examples @ git+https://github.com/open-atmos/PySDM-examples@885b82a#egg=PySDM-examples
 PyMPDATA @ git+https://github.com/open-atmos/PyMPDATA@63f1c11#egg=PyMPDATA
 
 PyPartMC==0.0.24

tests/unit_tests/backends/test_freezing_methods.py

@@ -1,5 +1,6 @@
 # pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
 import numpy as np
+import pytest
 from matplotlib import pyplot
 
 from PySDM import Builder, Formulae
@@ -23,6 +24,47 @@ def test_record_freezing_temperature_on_time_dependent_freeze(self):
     def test_no_subsaturated_freezing(self):
         pass
 
+    @staticmethod
+    @pytest.mark.parametrize("singular", (True, False))
+    @pytest.mark.parametrize("thaw", (True, False))
+    @pytest.mark.parametrize("epsilon", (0, 1e-5))
+    # pylint: disable=redefined-outer-name
+    def test_thaw(backend_class, singular, thaw, epsilon):
+        # arrange
+        formulae = Formulae()
+        builder = Builder(n_sd=1, backend=backend_class(formulae=formulae))
+        env = Box(dt=1 * si.s, dv=1 * si.m**3)
+        builder.set_environment(env)
+        builder.add_dynamic(Freezing(singular=singular, thaw=thaw))
+        particulator = builder.build(
+            products=(IceWaterContent(),),
+            attributes={
+                "n": np.ones(builder.particulator.n_sd),
+                "volume": -1 * np.ones(builder.particulator.n_sd) * si.um**3,
+                **(
+                    {"freezing temperature": np.full(builder.particulator.n_sd, -1)}
+                    if singular
+                    else {
+                        "immersed surface area": np.full(builder.particulator.n_sd, -1)
+                    }
+                ),
+            },
+        )
+        env["T"] = formulae.constants.T0 + epsilon
+        env["RH"] = np.nan
+        if not singular:
+            env["a_w_ice"] = np.nan
+        assert particulator.products["ice water content"].get() > 0
+
+        # act
+        particulator.run(steps=1)
+
+        # assert
+        if thaw and epsilon > 0:
+            assert particulator.products["ice water content"].get() == 0
+        else:
+            assert particulator.products["ice water content"].get() > 0
+
     @staticmethod
     # pylint: disable=redefined-outer-name
     def test_freeze_singular(backend_class):
@@ -121,6 +163,7 @@ def low(t):
             particulator = builder.build(attributes=attributes, products=products)
             env["RH"] = 1.0001
             env["a_w_ice"] = np.nan
+            env["T"] = np.nan
 
             cell_id = 0
             for i in range(int(total_time / case["dt"]) + 1):