Homogeneous freezing (as a new option in Freezing dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: Spichtinger_et_al_2023 (#1488)

Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl>

Author: 2 people

PySDM/backends/impl_common/freezing_attributes.py

@@ -37,6 +37,7 @@ class TimeDependentHomogeneousAttributes(
     """groups attributes required in time-dependent regime for homogeneous freezing"""
 
     __slots__ = ()
+<<<<<<< HEAD
 
 
 class ThresholdHomogeneousAndThawAttributes(
@@ -48,3 +49,5 @@ class ThresholdHomogeneousAndThawAttributes(
     """groups attributes required in time-dependent regime for homogeneous freezing"""
 
     __slots__ = ()
+=======
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))

PySDM/backends/impl_numba/methods/freezing_methods.py

@@ -14,7 +14,10 @@
     SingularAttributes,
     TimeDependentAttributes,
     TimeDependentHomogeneousAttributes,
+<<<<<<< HEAD
     ThresholdHomogeneousAndThawAttributes,
+=======
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
 )
 
 
@@ -38,6 +41,7 @@ def body(signed_water_mass, i):
         return body
 
     @cached_property
+<<<<<<< HEAD
     def _thaw_instantaneous_body(self):
         _thaw = self._thaw
         frozen_and_above_freezing_point = (
@@ -67,6 +71,18 @@ def body(
             relative_humidity,
             cell,
         ):
+=======
+    def _freeze_singular_body(self):
+        _thaw = self._thaw
+        _freeze = self._freeze
+        frozen_and_above_freezing_point = (
+            self.formulae.trivia.frozen_and_above_freezing_point
+        )
+        unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated
+
+        @numba.njit(**self.default_jit_flags)
+        def body(attributes, temperature, relative_humidity, cell, thaw):
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             n_sd = len(attributes.freezing_temperature)
             for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
                 if attributes.freezing_temperature[i] == 0:
@@ -82,8 +98,17 @@ def body(
         return body
 
     @cached_property
+<<<<<<< HEAD
     def _immersion_freezing_time_dependent_body(self):
         _freeze = self._freeze
+=======
+    def _freeze_time_dependent_body(self):
+        _thaw = self._thaw
+        _freeze = self._freeze
+        frozen_and_above_freezing_point = (
+            self.formulae.trivia.frozen_and_above_freezing_point
+        )
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
         unfrozen_and_saturated = self.formulae.trivia.unfrozen_and_saturated
         j_het = self.formulae.heterogeneous_ice_nucleation_rate.j_het
         prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
@@ -115,6 +140,66 @@ def body(  # pylint: disable=too-many-arguments
                         _freeze(attributes.signed_water_mass, i)
 
         return body
+<<<<<<< HEAD
+=======
+
+    @cached_property
+    def _freeze_time_dependent_homogeneous_body(self):
+        _thaw = self._thaw
+        _freeze = self._freeze
+        frozen_and_above_freezing_point = (
+            self.formulae.trivia.frozen_and_above_freezing_point
+        )
+        unfrozen_and_ice_saturated = self.formulae.trivia.unfrozen_and_ice_saturated
+        j_hom = self.formulae.homogeneous_ice_nucleation_rate.j_hom
+        prob_zero_events = self.formulae.trivia.poissonian_avoidance_function
+        d_a_w_ice_within_range = (
+            self.formulae.homogeneous_ice_nucleation_rate.d_a_w_ice_within_range
+        )
+        d_a_w_ice_maximum = (
+            self.formulae.homogeneous_ice_nucleation_rate.d_a_w_ice_maximum
+        )
+
+        @numba.njit(**self.default_jit_flags)
+        def body(  # pylint: disable=unused-argument,too-many-arguments
+            rand,
+            attributes,
+            timestep,
+            cell,
+            a_w_ice,
+            temperature,
+            relative_humidity_ice,
+            thaw,
+        ):
+
+            n_sd = len(attributes.signed_water_mass)
+            for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
+                cell_id = cell[i]
+                if thaw and frozen_and_above_freezing_point(
+                    attributes.signed_water_mass[i], temperature[cell_id]
+                ):
+                    _thaw(attributes.signed_water_mass, i)
+                elif unfrozen_and_ice_saturated(
+                    attributes.signed_water_mass[i], relative_humidity_ice[cell_id]
+                ):
+                    d_a_w_ice = (relative_humidity_ice[cell_id] - 1.0) * a_w_ice[
+                        cell_id
+                    ]
+
+                    if d_a_w_ice_within_range(d_a_w_ice):
+                        d_a_w_ice = d_a_w_ice_maximum(d_a_w_ice)
+                        rate_assuming_constant_temperature_within_dt = (
+                            j_hom(temperature[cell_id], d_a_w_ice)
+                            * attributes.volume[i]
+                        )
+                        prob = 1 - prob_zero_events(
+                            r=rate_assuming_constant_temperature_within_dt, dt=timestep
+                        )
+                        if rand[i] < prob:
+                            _freeze(attributes.signed_water_mass, i)
+
+        return body
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
 
     @cached_property
     def _homogeneous_freezing_time_dependent_body(self):
@@ -190,6 +275,7 @@ def thaw_instantaneous(
         cell,
         temperature,
     ):
+<<<<<<< HEAD
         self._thaw_instantaneous_body(
             ThresholdHomogeneousAndThawAttributes(
                 signed_water_mass=attributes.signed_water_mass.data,
@@ -202,6 +288,9 @@ def immersion_freezing_singular(
         self, *, attributes, temperature, relative_humidity, cell
     ):
         self._immersion_freezing_singular_body(
+=======
+        self._freeze_singular_body(
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             SingularAttributes(
                 freezing_temperature=attributes.freezing_temperature.data,
                 signed_water_mass=attributes.signed_water_mass.data,
@@ -221,7 +310,11 @@ def immersion_freezing_time_dependent(
         a_w_ice,
         relative_humidity,
     ):
+<<<<<<< HEAD
         self._immersion_freezing_time_dependent_body(
+=======
+        self._freeze_time_dependent_body(
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             rand.data,
             TimeDependentAttributes(
                 immersed_surface_area=attributes.immersed_surface_area.data,
@@ -274,6 +367,32 @@ def homogeneous_freezing_time_dependent(
             relative_humidity_ice.data,
         )
 
+    def freeze_time_dependent_homogeneous(
+        self,
+        *,
+        rand,
+        attributes,
+        timestep,
+        cell,
+        a_w_ice,
+        temperature,
+        relative_humidity_ice,
+        thaw: bool,
+    ):
+        self._freeze_time_dependent_homogeneous_body(
+            rand.data,
+            TimeDependentHomogeneousAttributes(
+                volume=attributes.volume.data,
+                signed_water_mass=attributes.signed_water_mass.data,
+            ),
+            timestep,
+            cell.data,
+            a_w_ice.data,
+            temperature.data,
+            relative_humidity_ice.data,
+            thaw=thaw,
+        )
+
     @cached_property
     def _record_freezing_temperatures_body(self):
         ff = self.formulae_flattened

PySDM/dynamics/freezing.py

@@ -13,6 +13,7 @@ class Freezing:  # pylint: disable=too-many-instance-attributes
     def __init__(
         self,
         *,
+<<<<<<< HEAD
         homogeneous_freezing: Optional[str] = None,
         immersion_freezing: Optional[str] = None,
         thaw: Optional[str] = None,
@@ -30,6 +31,14 @@ def __init__(
         )
         assert thaw is None or thaw == "instantaneous"
 
+=======
+        singular=True,
+        homogeneous_freezing=False,
+        immersion_freezing=True,
+        thaw=False,
+    ):
+        self.singular = singular
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
         self.homogeneous_freezing = homogeneous_freezing
         self.immersion_freezing = immersion_freezing
         self.thaw = thaw
@@ -49,24 +58,36 @@ def register(self, builder):
         if self.immersion_freezing == "singular":
             builder.request_attribute("freezing temperature")
 
+<<<<<<< HEAD
         if self.immersion_freezing == "time-dependent":
+=======
+        if not self.singular and self.immersion_freezing:
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             assert (
                 self.particulator.formulae.heterogeneous_ice_nucleation_rate.__name__
                 != "Null"
             )
             builder.request_attribute("immersed surface area")
 
+<<<<<<< HEAD
         if self.homogeneous_freezing == "time-dependent":
+=======
+        if self.homogeneous_freezing:
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             assert (
                 self.particulator.formulae.homogeneous_ice_nucleation_rate.__name__
                 != "Null"
             )
             builder.request_attribute("volume")
 
+<<<<<<< HEAD
         if (
             self.homogeneous_freezing == "time-dependent"
             or self.immersion_freezing == "time-dependent"
         ):
+=======
+        if self.homogeneous_freezing or not self.singular:
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             self.rand = self.particulator.Storage.empty(
                 self.particulator.n_sd, dtype=float
             )
@@ -84,11 +105,25 @@ def __call__(self):
         if not self.enable:
             return
 
+<<<<<<< HEAD
         if self.immersion_freezing == "singular":
             self.particulator.immersion_freezing_singular()
         elif self.immersion_freezing == "time-dependent":
+=======
+        if self.immersion_freezing:
+            if self.singular:
+                self.particulator.immersion_freezing_singular(thaw=self.thaw)
+            else:
+                self.rand.urand(self.rng)
+                self.particulator.immersion_freezing_time_dependent(
+                    rand=self.rand,
+                    thaw=self.thaw,
+                )
+
+        if self.homogeneous_freezing:
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             self.rand.urand(self.rng)
-            self.particulator.immersion_freezing_time_dependent(
+            self.particulator.homogeneous_freezing_time_dependent(
                 rand=self.rand,
             )
 

PySDM/particulator.py

@@ -9,7 +9,10 @@
     SingularAttributes,
     TimeDependentAttributes,
     TimeDependentHomogeneousAttributes,
+<<<<<<< HEAD
     ThresholdHomogeneousAndThawAttributes,
+=======
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
 )
 from PySDM.backends.impl_common.index import make_Index
 from PySDM.backends.impl_common.indexed_storage import make_IndexedStorage
@@ -533,8 +536,13 @@ def immersion_freezing_singular(self):
         )
         self.attributes.mark_updated("signed water mass")
 
+<<<<<<< HEAD
     def homogeneous_freezing_time_dependent(self, *, rand: Storage):
         self.backend.homogeneous_freezing_time_dependent(
+=======
+    def homogeneous_freezing_time_dependent(self, *, thaw: bool, rand: Storage):
+        self.backend.freeze_time_dependent_homogeneous(
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
             rand=rand,
             attributes=TimeDependentHomogeneousAttributes(
                 volume=self.attributes["volume"],
@@ -545,6 +553,7 @@ def homogeneous_freezing_time_dependent(self, *, rand: Storage):
             a_w_ice=self.environment["a_w_ice"],
             temperature=self.environment["T"],
             relative_humidity_ice=self.environment["RH_ice"],
+<<<<<<< HEAD
         )
 
     def homogeneous_freezing_threshold(self):
@@ -564,4 +573,7 @@ def thaw_instantaneous(self):
             ),
             cell=self.attributes["cell id"],
             temperature=self.environment["T"],
+=======
+            thaw=thaw,
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
         )

PySDM/physics/constants_defaults.py

@@ -104,10 +104,13 @@
 """ thermal accommodation coefficient for vapour deposition as recommended in
 [Pruppacher & Klett](https://doi.org/10.1007/978-0-306-48100-0) """
 
+<<<<<<< HEAD
 C_cunn = 0.7
 """ Cunningham correction factor as used in
 [Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) """
 
+=======
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
 ARM_C1 = 6.1094 * si.hectopascal
 """ [August](https://doi.org/10.1002/andp.18280890511) Roche Magnus formula coefficients
 (values from [Alduchov & Eskridge 1996](https://doi.org/10.1175%2F1520-0450%281996%29035%3C0601%3AIMFAOS%3E2.0.CO%3B2))

PySDM/physics/homogeneous_ice_nucleation_rate/constant.py

@@ -5,7 +5,11 @@
 import numpy as np
 
 
+<<<<<<< HEAD
 class Constant:
+=======
+class Constant:  # pylint: disable=too-few-public-methods
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
     def __init__(self, const):
         assert np.isfinite(const.J_HOM)
 
@@ -18,5 +22,9 @@ def d_a_w_ice_maximum(const, da_w_ice):  # pylint: disable=unused-argument
         return da_w_ice
 
     @staticmethod
+<<<<<<< HEAD
     def j_hom(const, T, da_w_ice):  # pylint: disable=unused-argument
+=======
+    def j_hom(const, T, a_w_ice):  # pylint: disable=unused-argument
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
         return const.J_HOM

PySDM/physics/homogeneous_ice_nucleation_rate/koop.py

@@ -7,7 +7,11 @@
 import numpy as np
 
 
+<<<<<<< HEAD
 class Koop2000:
+=======
+class Koop2000:  # pylint: disable=too-few-public-methods
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
     def __init__(self, const):
         pass
 

PySDM/physics/homogeneous_ice_nucleation_rate/koop_corr.py

@@ -8,7 +8,11 @@
 import numpy as np
 
 
+<<<<<<< HEAD
 class Koop_Correction:
+=======
+class Koop_Correction:  # pylint: disable=too-few-public-methods
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
     def __init__(self, const):
         pass
 

PySDM/physics/homogeneous_ice_nucleation_rate/koop_murray.py

@@ -7,7 +7,11 @@
 import numpy as np
 
 
+<<<<<<< HEAD
 class KoopMurray2016:
+=======
+class KoopMurray2016:  # pylint: disable=too-few-public-methods
+>>>>>>> 561d4927 (Homogeneous freezing (as a new option in `Freezing` dynamic, disabled by default + new physics formulae for hom. nucl. rate); new example: `Spichtinger_et_al_2023` (#1488))
     def __init__(self, const):
         pass