Homogeneous freezing

PySDM/backends/impl_common/freezing_attributes.py

@@ -19,3 +19,11 @@ class TimeDependentAttributes(
     """groups attributes required in time-dependent regime"""
 
     __slots__ = ()
+
+
+class TimeDependentHomogeneousAttributes(
+    namedtuple("TimeDependentHomogeneousAttributes", ("volume", "water_mass"))
+):
+    """groups attributes required in time-dependent regime for homogeneous freezing"""
+
+    __slots__ = ()

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -1,5 +1,5 @@
 """
-CPU implementation of backend methods for freezing (singular and time-dependent immersion freezing)
+CPU implementation of backend methods for homogeneous freezing and heterogeneous freezing (singular and time-dependent immersion freezing)
 """
 
 import numba
@@ -10,13 +10,15 @@
 from ...impl_common.freezing_attributes import (
     SingularAttributes,
     TimeDependentAttributes,
+    TimeDependentHomogeneousAttributes,
 )
 
 
 class FreezingMethods(BackendMethods):
     def __init__(self):
         BackendMethods.__init__(self)
         unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated
+        unfrozen_and_ice_saturated = self.formulae.trivia.unfrozen_and_ice_saturated
         frozen_and_above_freezing_point = (
             self.formulae.trivia.frozen_and_above_freezing_point
         )
@@ -92,6 +94,48 @@ def freeze_time_dependent_body(  # pylint: disable=unused-argument,too-many-argu
 
         self.freeze_time_dependent_body = freeze_time_dependent_body
 
+
+        j_hom = self.formulae.homogeneous_ice_nucleation_rate.j_hom
+
+        @numba.njit(**self.default_jit_flags)
+        def freeze_time_dependent_homogeneous_body(  # pylint: disable=unused-argument,too-many-arguments
+            rand,
+            attributes,
+            timestep,
+            cell,
+            a_w_ice,
+            temperature,
+            relative_humidity_ice,
+            record_freezing_temperature,
+            freezing_temperature,
+            thaw,
+        ):
+
+            n_sd = len(attributes.water_mass)
+            for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
+                cell_id = cell[i]
+                #relative_humidity_ice = relative_humidity[cell_id] / a_w_ice[cell_id]
+                if thaw and frozen_and_above_freezing_point(
+                    attributes.water_mass[i], temperature[cell_id]
+                ):
+                    _thaw(attributes.water_mass, i)
+                elif unfrozen_and_ice_saturated(
+                    attributes.water_mass[i], relative_humidity_ice[cell_id]
+                ):
+                    d_a_w_ice = (relative_humidity_ice[cell_id] - 1) * a_w_ice[cell_id]
+                    rate = j_hom(temperature, d_a_w_ice)
+                    # TODO #594: this assumes constant T throughout timestep, can we do better?
+                    prob = 1 - np.exp(  # TODO #599: common code for Poissonian prob
+                        -rate * attributes.volume[i] * timestep
+                    )
+                    if rand[i] < prob:
+                        _freeze(attributes.water_mass, i)
+                        # if record_freezing_temperature:
+                        #     freezing_temperature[i] = temperature[cell_id]
+
+
+        self.freeze_time_dependent_homogeneous_body = freeze_time_dependent_homogeneous_body
+
     def freeze_singular(
         self, *, attributes, temperature, relative_humidity, cell, thaw: bool
     ):
@@ -106,6 +150,8 @@ def freeze_singular(
             thaw=thaw,
         )
 
+
+
     def freeze_time_dependent(
         self,
         *,
@@ -137,3 +183,38 @@ def freeze_time_dependent(
             ),
             thaw=thaw,
         )
+
+
+
+    def freeze_time_dependent_homogeneous(
+        self,
+        *,
+        rand,
+        attributes,
+        timestep,
+        cell,
+        a_w_ice,
+        temperature,
+        relative_humidity_ice,
+        record_freezing_temperature,
+        freezing_temperature,
+        thaw: bool
+    ):
+        self.freeze_time_dependent_homogeneous_body(
+            rand.data,
+            TimeDependentHomogeneousAttributes(
+                volume=attributes.volume.data,
+                water_mass=attributes.water_mass.data,
+            ),
+            timestep,
+            cell.data,
+            a_w_ice.data,
+            temperature.data,
+            relative_humidity_ice.data,
+            record_freezing_temperature=record_freezing_temperature,
+            freezing_temperature=(
+                freezing_temperature.data if record_freezing_temperature else None
+            ),
+            thaw=thaw,
+        )
+

PySDM/dynamics/freezing.py

@@ -1,20 +1,23 @@
 """
-immersion freezing using either singular or time-dependent formulation
+homogeneous freezing and heterogeneous freezing (singular and time-dependent immersion freezing)
 """
 
 from PySDM.backends.impl_common.freezing_attributes import (
     SingularAttributes,
     TimeDependentAttributes,
+    TimeDependentHomogeneousAttributes,
 )
 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):
+    def __init__(self, *, singular=True, homogeneous_freezing=False, immersion_freezing=True, record_freezing_temperature=False, thaw=False):
         assert not (record_freezing_temperature and singular)
         self.singular = singular
+        self.homogeneous_freezing = homogeneous_freezing
+        self.immersion_freezing = immersion_freezing
         self.record_freezing_temperature = record_freezing_temperature
         self.thaw = thaw
         self.enable = True
@@ -35,14 +38,17 @@ def register(self, builder):
             assert not isinstance(
                 builder.formulae.heterogeneous_ice_nucleation_rate, Null
             )
-            builder.request_attribute("immersed surface area")
+            if self.immersion_freezing:
+                builder.request_attribute("immersed surface area")
             self.rand = self.particulator.Storage.empty(
                 self.particulator.n_sd, dtype=float
             )
             self.rng = self.particulator.Random(
                 self.particulator.n_sd, self.particulator.backend.formulae.seed
             )
 
+
+
     def __call__(self):
         if "Coalescence" in self.particulator.dynamics:
             # TODO #594
@@ -53,34 +59,60 @@ def __call__(self):
         if not self.enable:
             return
 
-        if self.singular:
-            self.particulator.backend.freeze_singular(
-                attributes=SingularAttributes(
-                    freezing_temperature=self.particulator.attributes[
-                        "freezing temperature"
-                    ],
-                    water_mass=self.particulator.attributes["water mass"],
-                ),
-                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)
-            self.particulator.backend.freeze_time_dependent(
-                rand=self.rand,
-                attributes=TimeDependentAttributes(
+        if self.immersion_freezing:
+
+            if self.singular:
+                self.particulator.backend.freeze_singular(
+                    attributes=SingularAttributes(
+                        freezing_temperature=self.particulator.attributes[
+                            "freezing temperature"
+                        ],
+                        water_mass=self.particulator.attributes["water mass"],
+                    ),
+                    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)
+                self.particulator.backend.freeze_time_dependent(
+                    rand=self.rand,
+                    attributes=TimeDependentAttributes(
                     immersed_surface_area=self.particulator.attributes[
                         "immersed surface area"
                     ],
                     water_mass=self.particulator.attributes["water mass"],
+                    ),
+                    timestep=self.particulator.dt,
+                    cell=self.particulator.attributes["cell id"],
+                    a_w_ice=self.particulator.environment["a_w_ice"],
+                    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,
+                )
+
+
+        if self.homogeneous_freezing:
+
+            self.rand.urand(self.rng)
+            self.particulator.backend.freeze_time_dependent_homogeneous(
+                rand=self.rand,
+                attributes=TimeDependentHomogeneousAttributes(
+                volume=self.particulator.attributes["volume"],
+                water_mass=self.particulator.attributes["water mass"],
                 ),
                 timestep=self.particulator.dt,
                 cell=self.particulator.attributes["cell id"],
                 a_w_ice=self.particulator.environment["a_w_ice"],
                 temperature=self.particulator.environment["T"],
-                relative_humidity=self.particulator.environment["RH"],
+                relative_humidity_ice=self.particulator.environment["RH_ice"],
                 record_freezing_temperature=self.record_freezing_temperature,
                 freezing_temperature=(
                     self.particulator.attributes["freezing temperature"]
@@ -90,6 +122,8 @@ def __call__(self):
                 thaw=self.thaw,
             )
 
+
+
         self.particulator.attributes.mark_updated("water mass")
         if self.record_freezing_temperature:
             self.particulator.attributes.mark_updated("freezing temperature")

PySDM/formulae.py

@@ -44,6 +44,7 @@ def __init__(  # pylint: disable=too-many-locals
         hydrostatics: str = "Default",
         freezing_temperature_spectrum: str = "Null",
         heterogeneous_ice_nucleation_rate: str = "Null",
+        homogeneous_ice_nucleation_rate: str = "Null",
         fragmentation_function: str = "AlwaysN",
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_meteoric_water_line_excess: str = "Null",
@@ -76,6 +77,7 @@ def __init__(  # pylint: disable=too-many-locals
         self.hydrostatics = hydrostatics
         self.freezing_temperature_spectrum = freezing_temperature_spectrum
         self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate
+        self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate
         self.fragmentation_function = fragmentation_function
         self.isotope_equilibrium_fractionation_factors = (
             isotope_equilibrium_fractionation_factors

PySDM/physics/__init__.py

@@ -25,6 +25,7 @@
     fragmentation_function,
     freezing_temperature_spectrum,
     heterogeneous_ice_nucleation_rate,
+    homogeneous_ice_nucleation_rate,
     hydrostatics,
     hygroscopicity,
     impl,

PySDM/physics/constants_defaults.py

@@ -185,7 +185,26 @@
 ABIFM_C = np.inf
 ABIFM_UNIT = 1 / si.cm**2 / si.s
 
+KOOP_2000_C1 = -906.7
+KOOP_2000_C2 = 8502
+KOOP_2000_C3 = 26924
+KOOP_2000_C4 = 29180
+KOOP_UNIT = 1 / si.cm**3 / si.s
+
+KOOP_CORR = -1.522
+
+KOOP_MURRAY_C0 = -3020.684
+KOOP_MURRAY_C1 = -425.921
+KOOP_MURRAY_C2 = -25.9779
+KOOP_MURRAY_C3 = -0.868451
+KOOP_MURRAY_C4 = -1.66203e-2
+KOOP_MURRAY_C5 = -1.71736e-4
+KOOP_MURRAY_C6 = -7.46953e-7
+
+T_hom_freeze = 235. * si.kelvin
+
 J_HET = np.nan
+J_HOM = np.nan
 
 STRAUB_E_D1 = 0.04 * si.cm
 STRAUB_MU2 = 0.095 * si.cm

PySDM/physics/homogeneous_ice_nucleation_rate/__init__.py

@@ -0,0 +1,7 @@
+"""
+homogeneous-freezing rate (aka J_hom) formulations
+"""
+
+from .constant import Constant
+from .null import Null
+from .koop import KOOP

PySDM/physics/homogeneous_ice_nucleation_rate/constant.py

@@ -0,0 +1,14 @@
+"""
+constant rate formulation (for tests)
+"""
+
+import numpy as np
+
+
+class Constant:  # pylint: disable=too-few-public-methods
+    def __init__(self, const):
+        assert np.isfinite(const.J_HOM)
+
+    @staticmethod
+    def j_hom(const, T, a_w_ice):  # pylint: disable=unused-argument
+        return const.J_HOM

PySDM/physics/homogeneous_ice_nucleation_rate/koop.py

@@ -0,0 +1,14 @@
+"""
+Koop homogeneous nucleation rate parameterization for solution droplets 
+valid for 0.26 < da_w_ice < 0.34
+ ([Koop et al. 2000](https://doi.org/10.1038/35020537))
+"""
+
+
+class KOOP:  # pylint: disable=too-few-public-methods
+    def __init__(self, const):
+        pass
+
+    @staticmethod
+    def j_hom(const, T, da_w_ice):
+        return 10**(const.KOOP_2000_C1 + const.KOOP_2000_C2 * da_w_ice - const.KOOP_2000_C3 * da_w_ice**2. + const.KOOP_2000_C4 * da_w_ice**3.) * const.KOOP_UNIT

PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py

@@ -0,0 +1,15 @@
+"""
+Koop homogeneous nucleation rate parameterization for solution droplets [Koop et al. 2000] corrected 
+such that it coincides with homogeneous nucleation rate parameterization for pure water droplets 
+[Koop and Murray 2016] at water saturation between 235K < T < 240K
+ ([Spichtinger et al. 2023](https://doi.org/10.5194/acp-23-2035-2023))
+"""
+
+
+class KOOP_CORR:  # pylint: disable=too-few-public-methods
+    def __init__(self, const):
+        pass
+
+    @staticmethod
+    def j_hom(const, T, da_w_ice):
+        return 10**(const.KOOP_2000_C1 + const.KOOP_2000_C2 * da_w_ice - const.KOOP_2000_C3 * da_w_ice**2. + const.KOOP_2000_C4 * da_w_ice**3. + const.KOOP_CORR) * const.KOOP_UNIT

PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray.py

@@ -0,0 +1,15 @@
+"""
+Koop and Murray homogeneous nucleation rate parameterization for pure water droplets 
+at water saturation
+ ([Koop and Murray 2016](https://doi.org/10.1063/1.4962355))
+"""
+
+
+class KOOP_MURRAY:  # pylint: disable=too-few-public-methods
+    def __init__(self, const):
+        pass
+
+    @staticmethod
+    def j_hom(const, T, da_w_ice):
+        T_diff = T - const.T_tri
+        return 10**(const.KOOP_MURRAY_C0 + const.KOOP_MURRAY_C1*T_diff + const.KOOP_MURRAY_C2*T_diff**2. + const.KOOP_MURRAY_C3*T_diff**3. + const.KOOP_MURRAY_C4*T_diff**4. + const.KOOP_MURRAY_C5*T_diff**5. + const.KOOP_MURRAY_C6*T_diff**6.) * const.KOOP_UNIT