Merge pull request #613 from slayoo/freezing

introducing immersed surface area attribute

Author: slayoo

PySDM/attributes/ice/__init__.py

@@ -1,2 +1,2 @@
 from .freezing_temperature import FreezingTemperature
-from .nucleation_sites import NucleationSites
+from .immersed_surface_area import ImmersedSurfaceArea

PySDM/attributes/ice/immersed_surface_area.py

@@ -0,0 +1,5 @@
+from ..impl import ExtensiveAttribute
+
+class ImmersedSurfaceArea(ExtensiveAttribute):
+    def __init__(self, particles_builder):
+        super().__init__(particles_builder, name='immersed surface area')

PySDM/attributes/ice/nucleation_sites.py

@@ -3,7 +3,7 @@
 from PySDM.attributes.physics.hygroscopicity import Kappa, KappaTimesDryVolume
 from PySDM.attributes.physics import (Multiplicities, Volume, Radius, DryRadius,
                                       TerminalVelocity, Temperature, Heat, CriticalVolume)
-from PySDM.attributes.ice import FreezingTemperature, NucleationSites
+from PySDM.attributes.ice import FreezingTemperature, ImmersedSurfaceArea
 from PySDM.attributes.numerics import CellID, CellOrigin, PositionInCell
 from PySDM.attributes.chemistry import MoleAmount, Concentration, pH, HydrogenIonConcentration
 from PySDM.physics.aqueous_chemistry.support import AQUEOUS_COMPOUNDS
@@ -41,7 +41,7 @@
     'pH': lambda _: pH,
     'conc_H': lambda _: HydrogenIonConcentration,
     'freezing temperature': lambda _: FreezingTemperature,
-    'nucleation sites': lambda _: NucleationSites
+    'immersed surface area': lambda _: ImmersedSurfaceArea
 }
 
 

PySDM/attributes/impl/mapper.py

@@ -24,17 +24,17 @@ def freeze_singular_body(T_fz, v_wet, T, RH, cell):
         self.freeze_singular_body = freeze_singular_body
 
         @numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath})
-        def freeze_time_dependent_body(rand, nucleation_sites, volume, r, dt):
-            p = 1 - np.exp(-r * dt)
-            # TODO #599: assert if > 1?
-            for i in numba.prange(len(nucleation_sites)):
+        def freeze_time_dependent_body(rand, immersed_surface_area, volume, J_het, dt):
+            n_sd = len(volume)
+            for i in numba.prange(n_sd):
                 if _unfrozen(volume, i):
+                    p = 1 - np.exp(-J_het * immersed_surface_area[i] * dt)  # TODO #599: assert if > 1?
                     if rand[i] < p:
                         _freeze(volume, i)
         self.freeze_time_dependent_body = freeze_time_dependent_body
 
     def freeze_singular(self, T_fz, v_wet, T, RH, cell):
         self.freeze_singular_body(T_fz.data, v_wet.data, T.data, RH.data, cell.data)
 
-    def freeze_time_dependent(self, rand, nucleation_sites, volume, r, dt):
-        self.freeze_time_dependent_body(rand.data, nucleation_sites.data, volume.data, r, dt)
+    def freeze_time_dependent(self, rand, immersed_surface_area, volume, J_het, dt):
+        self.freeze_time_dependent_body(rand.data, immersed_surface_area.data, volume.data, J_het, dt)

PySDM/backends/numba/impl/freezing_methods.py

@@ -1,8 +1,8 @@
 class Freezing:
-    def __init__(self, *, singular=True, r=None):
+    def __init__(self, *, singular=True, J_het=None):
         self.singular = singular
         self.enable = True
-        self.r = r
+        self.J_het = J_het
         self.rand = None
         self.rng = None
 
@@ -11,9 +11,11 @@ def register(self, builder):
 
         builder.request_attribute("volume")
         if self.singular:
+            assert self.J_het is None
             builder.request_attribute("freezing temperature")
         else:
-            builder.request_attribute("nucleation sites")
+            assert self.J_het is not None
+            builder.request_attribute("immersed surface area")
             self.rand = self.core.Storage.empty(self.core.n_sd, dtype=float)
             self.rng = self.core.Random(self.core.n_sd, self.core.bck.formulae.seed)
 
@@ -36,9 +38,9 @@ def __call__(self):
             self.rand.urand(self.rng)
             self.core.bck.freeze_time_dependent(
                 rand=self.rand,
-                nucleation_sites=self.core.particles['nucleation sites'],
+                immersed_surface_area=self.core.particles['immersed surface area'],
                 volume=self.core.particles['volume'],
-                r=self.r,
+                J_het=self.J_het,
                 dt=self.core.dt
             )
 

PySDM/dynamics/freezing.py

@@ -49,6 +49,9 @@ class Lognormal(Spectrum):
     def __init__(self, norm_factor: float, m_mode: float, s_geom: float):
         super().__init__(lognorm, (math.log(s_geom), 0, m_mode), norm_factor)
 
+    @property
+    def s_geom(self):
+        return math.exp(self.distribution_params[0])
 
 class Sum:
 

PySDM/physics/spectra.py

@@ -19,7 +19,7 @@ class TestFreezingMethods:
     #     pass
 
     @staticmethod
-    def test_freeze_time_dependent(plot=False):
+    def test_freeze_time_dependent(plot=True):
         # Arrange
         cases = (
             {'dt': 5e5, 'N':  1},
@@ -29,17 +29,18 @@ def test_freeze_time_dependent(plot=False):
             {'dt': 5e5, 'N': 32},
             {'dt': 1e6, 'N': 32},
         )
-        r = 1e-9
+        rate = 1e-9
+        immersed_surface_area = 1
+        J_het = rate / immersed_surface_area
         number_of_real_droplets = 1024
-        nucleation_sites_per_particle = 1024
         total_time = 2e9  # effectively interpretted here as seconds, i.e. cycle = 1 * si.s
 
         # dummy (but must-be-set) values
         vol = 44  # just to enable sign flipping (ice water uses negative volumes), actual value does not matter
         dv = 666  # products use concentration, just dividing there and multiplying back here, actual value does not matter
 
-        hgh = lambda t: np.exp(-0.8 * r * (t - total_time / 10))
-        low = lambda t: np.exp(-1.2 * r * (t + total_time / 10))
+        hgh = lambda t: np.exp(-0.8 * rate * (t - total_time / 10))
+        low = lambda t: np.exp(-1.2 * rate * (t + total_time / 10))
 
         # Act
         output = {}
@@ -54,10 +55,10 @@ def test_freeze_time_dependent(plot=False):
 
             builder = Builder(n_sd=n_sd, backend=CPU)
             builder.set_environment(Box(dt=case['dt'], dv=dv))
-            builder.add_dynamic(Freezing(singular=False, r=r))
+            builder.add_dynamic(Freezing(singular=False, J_het=J_het))
             attributes = {
                 'n': np.full(n_sd, int(case['N'])),
-                'nucleation sites': np.full(n_sd, nucleation_sites_per_particle),
+                'immersed surface area': np.full(n_sd, immersed_surface_area),
                 'volume': np.full(n_sd, vol)
             }
             products = (IceWaterContent(specific=False),)
@@ -76,7 +77,7 @@ def test_freeze_time_dependent(plot=False):
         # Plot
         if plot:
             fit_x = np.linspace(0, total_time, num=100)
-            fit_y = np.exp(-r * fit_x)
+            fit_y = np.exp(-rate * fit_x)
 
             for key in output.keys():
                 pylab.step(

PySDM_tests/unit_tests/backends/test_freezing_methods.py

@@ -75,7 +75,7 @@ beyond the latest release).
 For development purposes, we suggest cloning the repository and installing it using ``pip -e``.
 Test-time dependencies are listed in the ``test-time-requirements.txt`` file.
 
-PySDM examples listed below are hosted in a separate repository and constitute 
+PySDM examples are hosted in a separate repository and constitute 
 the [``PySDM_examples``](https://github.com/atmos-cloud-sim-uj/PySDM-examples) package.
 The examples have additional dependencies listed in [``PySDM_examples`` package ``setup.py``](https://github.com/atmos-cloud-sim-uj/PySDM-examples/blob/main/setup.py) file.
 Running the examples requires the ``PySDM_examples`` package to be installed.
@@ -87,7 +87,7 @@ pip install -e .
 jupyter-notebook
 ```
 Alternatively, one can also install the examples package from pypi.org by 
-using `pip install PySDM-examples``.
+using ``pip install PySDM-examples``.
 
 ## PySDM examples (Jupyter notebooks reproducing results from literature):