recording temperature of last freezing for time-dependent immersion freezing (#1566)

Author: slayoo

PySDM/attributes/ice/freezing_temperature.py

@@ -1,12 +1,42 @@
 """
-particle freezing temperature (assigned at initialisation, modified through collisions only,
- used in singular regime)
+particle freezing temperature
 """
 
 from PySDM.attributes.impl import MaximumAttribute, register_attribute
+from ..impl import DerivedAttribute
 
 
 @register_attribute()
 class FreezingTemperature(MaximumAttribute):
+    """singular variant: assigned at initialisation, modified through collisions only"""
+
     def __init__(self, builder):
         super().__init__(builder, name="freezing temperature")
+
+
+@register_attribute()
+class TemperatureOfLastFreezing(DerivedAttribute):
+    """time-dependent variant: assigned upon freezing"""
+
+    def __init__(self, builder):
+        assert "Freezing" in builder.particulator.dynamics
+        assert not builder.particulator.dynamics["Freezing"].singular
+        self.signed_water_mass = builder.get_attribute("signed water mass")
+        self.cell_id = builder.get_attribute("cell id")
+        super().__init__(
+            builder,
+            name="temperature of last freezing",
+            dependencies=(self.signed_water_mass, self.cell_id),
+        )
+        builder.particulator.observers.append(self)
+
+    def notify(self):
+        self.update()
+
+    def recalculate(self):
+        self.particulator.backend.record_freezing_temperatures(
+            data=self.data,
+            cell_id=self.cell_id.data,
+            temperature=self.particulator.environment["T"],
+            signed_water_mass=self.signed_water_mass.data,
+        )

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -2,7 +2,10 @@
 CPU implementation of backend methods for freezing (singular and time-dependent immersion freezing)
 """
 
+from functools import cached_property
+
 import numba
+import numpy as np
 
 from PySDM.backends.impl_common.backend_methods import BackendMethods
 
@@ -56,16 +59,14 @@ def freeze_singular_body(
         prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
 
         @numba.njit(**self.default_jit_flags)
-        def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-arguments
+        def freeze_time_dependent_body(  # pylint: disable=too-many-arguments
             rand,
             attributes,
             timestep,
             cell,
             a_w_ice,
             temperature,
             relative_humidity,
-            record_freezing_temperature,
-            freezing_temperature,
             thaw,
         ):
             n_sd = len(attributes.signed_water_mass)
@@ -88,8 +89,6 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
                     )
                     if rand[i] < prob:
                         _freeze(attributes.signed_water_mass, i)
-                        # if record_freezing_temperature:
-                        #     freezing_temperature[i] = temperature[cell_id]
 
         self.freeze_time_dependent_body = freeze_time_dependent_body
 
@@ -117,8 +116,6 @@ def freeze_time_dependent(
         a_w_ice,
         temperature,
         relative_humidity,
-        record_freezing_temperature,
-        freezing_temperature,
         thaw: bool
     ):
         self.freeze_time_dependent_body(
@@ -132,9 +129,31 @@ def freeze_time_dependent(
             a_w_ice.data,
             temperature.data,
             relative_humidity.data,
-            record_freezing_temperature=record_freezing_temperature,
-            freezing_temperature=(
-                freezing_temperature.data if record_freezing_temperature else None
-            ),
             thaw=thaw,
         )
+
+    @cached_property
+    def _record_freezing_temperatures_body(self):
+        ff = self.formulae_flattened
+
+        @numba.njit(**{**self.default_jit_flags, "fastmath": False})
+        def body(data, cell_id, temperature, signed_water_mass):
+            for drop_id in numba.prange(len(data)):  # pylint: disable=not-an-iterable
+                if ff.trivia__unfrozen(signed_water_mass[drop_id]):
+                    if data[drop_id] > 0:
+                        data[drop_id] = np.nan
+                else:
+                    if np.isnan(data[drop_id]):
+                        data[drop_id] = temperature[cell_id[drop_id]]
+
+        return body
+
+    def record_freezing_temperatures(
+        self, *, data, cell_id, temperature, signed_water_mass
+    ):
+        self._record_freezing_temperatures_body(
+            data=data.data,
+            cell_id=cell_id.data,
+            temperature=temperature.data,
+            signed_water_mass=signed_water_mass.data,
+        )

PySDM/backends/impl_thrust_rtc/methods/freezing_methods.py

@@ -118,8 +118,6 @@ def freeze_time_dependent(  # pylint: disable=unused-argument
         a_w_ice,
         temperature,
         relative_humidity,
-        record_freezing_temperature,
-        freezing_temperature,
         thaw,
     ):
         n_sd = len(attributes.immersed_surface_area)

PySDM/dynamics/freezing.py

@@ -2,16 +2,13 @@
 immersion freezing using either singular or time-dependent formulation
 """
 
-from PySDM.physics.heterogeneous_ice_nucleation_rate import Null
 from PySDM.dynamics.impl import register_dynamic
 
 
 @register_dynamic()
 class Freezing:
-    def __init__(self, *, singular=True, record_freezing_temperature=False, thaw=False):
-        assert not (record_freezing_temperature and singular)
+    def __init__(self, *, singular=True, thaw=False):
         self.singular = singular
-        self.record_freezing_temperature = record_freezing_temperature
         self.thaw = thaw
         self.enable = True
         self.rand = None
@@ -26,12 +23,13 @@ def register(self, builder):
         )
 
         builder.request_attribute("signed water mass")
-        if self.singular or self.record_freezing_temperature:
+        if self.singular:
             builder.request_attribute("freezing temperature")
 
         if not self.singular:
-            assert not isinstance(
-                self.particulator.formulae.heterogeneous_ice_nucleation_rate, Null
+            assert (
+                self.particulator.formulae.heterogeneous_ice_nucleation_rate.__name__
+                != "Null"
             )
             builder.request_attribute("immersed surface area")
             self.rand = self.particulator.Storage.empty(
@@ -57,10 +55,5 @@ def __call__(self):
             self.rand.urand(self.rng)
             self.particulator.immersion_freezing_time_dependent(
                 rand=self.rand,
-                record_freezing_temperature=self.record_freezing_temperature,
                 thaw=self.thaw,
             )
-
-        self.particulator.attributes.mark_updated("signed water mass")
-        if self.record_freezing_temperature:
-            self.particulator.attributes.mark_updated("freezing temperature")

PySDM/particulator.py

@@ -509,9 +509,7 @@ def deposition(self):
         # TODO #1524 - should we update here?
         # self.update_TpRH(only_if_not_last='VapourDepositionOnIce')
 
-    def immersion_freezing_time_dependent(
-        self, *, thaw: bool, record_freezing_temperature: bool, rand: Storage
-    ):
+    def immersion_freezing_time_dependent(self, *, thaw: bool, rand: Storage):
         self.backend.freeze_time_dependent(
             rand=rand,
             attributes=TimeDependentAttributes(
@@ -523,14 +521,9 @@ def immersion_freezing_time_dependent(
             a_w_ice=self.environment["a_w_ice"],
             temperature=self.environment["T"],
             relative_humidity=self.environment["RH"],
-            record_freezing_temperature=record_freezing_temperature,
-            freezing_temperature=(
-                self.attributes["freezing temperature"]
-                if record_freezing_temperature
-                else None
-            ),
             thaw=thaw,
         )
+        self.attributes.mark_updated("signed water mass")
 
     def immersion_freezing_singular(self, *, thaw: bool):
         self.backend.freeze_singular(
@@ -543,3 +536,4 @@ def immersion_freezing_singular(self, *, thaw: bool):
             cell=self.attributes["cell id"],
             thaw=thaw,
         )
+        self.attributes.mark_updated("signed water mass")

tests/unit_tests/dynamics/test_immersion_freezing.py

@@ -8,12 +8,78 @@
 from PySDM.environments import Box
 from PySDM.physics import si
 from PySDM.products import IceWaterContent
+from PySDM.backends import GPU
+
+VERY_BIG_J_HET = 1e20
+EPSILON_RH = 1e-3
 
 
 class TestImmersionFreezing:
-    # TODO #599
-    def test_record_freezing_temperature_on_time_dependent_freeze(self):
-        pass
+    @staticmethod
+    @pytest.mark.parametrize(
+        "record_freezing_temperature",
+        (pytest.param(True, id="recording"), pytest.param(False, id="not recording")),
+    )
+    def test_record_freezing_temperature_on_time_dependent_freeze(
+        backend_class, record_freezing_temperature
+    ):
+        if backend_class is GPU and record_freezing_temperature:
+            pytest.skip("TODO #1495")
+
+        # arrange
+        formulae = Formulae(
+            particle_shape_and_density="MixedPhaseSpheres",
+            heterogeneous_ice_nucleation_rate="Constant",
+            constants={"J_HET": VERY_BIG_J_HET},
+        )
+        builder = Builder(
+            n_sd=1,
+            backend=backend_class(formulae=formulae),
+            environment=Box(dt=1 * si.s, dv=1 * si.m**3),
+        )
+        builder.add_dynamic(Freezing(singular=False, thaw=True))
+        if record_freezing_temperature:
+            builder.request_attribute("temperature of last freezing")
+        particulator = builder.build(
+            attributes={
+                "multiplicity": np.asarray([1]),
+                "signed water mass": np.asarray([1 * si.ug]),
+                "immersed surface area": np.asarray([1 * si.um**2]),
+            }
+        )
+
+        temp_1 = 200 * si.K
+        temp_2 = 250 * si.K
+        particulator.environment["a_w_ice"] = np.nan
+        particulator.environment["T"] = temp_1
+
+        # act & assert
+        attr_name = "temperature of last freezing"
+        if not record_freezing_temperature:
+            assert attr_name not in particulator.attributes
+        else:
+            # never frozen yet
+            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+
+            # still not frozen since RH not greater than 100%
+            particulator.environment["RH"] = 1.0
+            particulator.run(steps=1)
+            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+
+            # should freeze and record T1
+            particulator.environment["RH"] += EPSILON_RH
+            particulator.run(steps=1)
+            assert all(particulator.attributes[attr_name].to_ndarray() == temp_1)
+
+            # should thaw
+            particulator.environment["T"] = 300 * si.K
+            particulator.run(steps=1)
+            np.isnan(particulator.attributes[attr_name].to_ndarray()).all()
+
+            # should re-freeze and record T2
+            particulator.environment["T"] = temp_2
+            particulator.run(steps=1)
+            assert all(particulator.attributes[attr_name].to_ndarray() == temp_2)
 
     # TODO #599
     def test_no_subsaturated_freezing(self):
@@ -25,7 +91,17 @@ def test_no_subsaturated_freezing(self):
     @pytest.mark.parametrize("epsilon", (0, 1e-5))
     def test_thaw(backend_class, singular, thaw, epsilon):
         # arrange
-        formulae = Formulae(particle_shape_and_density="MixedPhaseSpheres")
+        formulae = Formulae(
+            particle_shape_and_density="MixedPhaseSpheres",
+            **(
+                {}
+                if singular
+                else {
+                    "heterogeneous_ice_nucleation_rate": "Constant",
+                    "constants": {"J_HET": 0},
+                }
+            ),
+        )
         env = Box(dt=1 * si.s, dv=1 * si.m**3)
         builder = Builder(
             n_sd=1, backend=backend_class(formulae=formulae), environment=env