Merge pull request #635 from slayoo/freezing

ABIFM implementation; set/getitem in Box env; T and a_w_ice dependence in J_het

PySDM/backends/numba/impl/freezing_methods.py

@@ -26,17 +26,17 @@ def freeze_singular_body(T_fz, v_wet, T, RH, cell):
         J_het = self.formulae.heterogeneous_ice_nucleation_rate.J_het
 
         @numba.njit(**{**conf.JIT_FLAGS, 'fastmath': self.formulae.fastmath})
-        def freeze_time_dependent_body(rand, immersed_surface_area, volume, dt):
+        def freeze_time_dependent_body(rand, immersed_surface_area, volume, dt, cell, a_w_ice):
             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?
+                    p = 1 - np.exp(-J_het(a_w_ice[cell[i]]) * 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, immersed_surface_area, volume, dt):
-        self.freeze_time_dependent_body(rand.data, immersed_surface_area.data, volume.data, dt)
+    def freeze_time_dependent(self, rand, immersed_surface_area, volume, dt, cell, a_w_ice):
+        self.freeze_time_dependent_body(rand.data, immersed_surface_area.data, volume.data, dt, cell.data, a_w_ice.data)

PySDM/dynamics/freezing.py

@@ -40,7 +40,9 @@ def __call__(self):
                 rand=self.rand,
                 immersed_surface_area=self.particulator.attributes['immersed surface area'],
                 volume=self.particulator.attributes['volume'],
-                dt=self.particulator.dt
+                dt=self.particulator.dt,
+                cell=self.particulator.attributes['cell id'],
+                a_w_ice=self.particulator.environment['a_w_ice']
             )
 
         self.particulator.attributes.attributes['volume'].mark_updated()

PySDM/environments/box.py

@@ -3,6 +3,7 @@
 """
 
 from PySDM.state.mesh import Mesh
+import numpy as np
 
 
 class Box:
@@ -11,6 +12,13 @@ def __init__(self, dt, dv):
         self.dt = dt
         self.mesh = Mesh.mesh_0d(dv)
         self.particulator = None
+        self._ambient_air = {}
+
+    def __getitem__(self, item):
+        return self._ambient_air[item]
+
+    def __setitem__(self, key, value):
+        self._ambient_air[key] = self.particulator.backend.Storage.from_ndarray(np.array([value]))
 
     def register(self, builder):
         self.particulator = builder.particulator

PySDM/physics/formulae.py

@@ -115,7 +115,7 @@ def __init__(self, *,
         self.hydrostatics = _magick(hydrostatics, physics.hydrostatics, fastmath)
         self.freezing_temperature_spectrum = _magick(freezing_temperature_spectrum, physics.freezing_temperature_spectrum, fastmath)
         self.heterogeneous_ice_nucleation_rate = _magick(heterogeneous_ice_nucleation_rate, physics.heterogeneous_ice_nucleation_rate, fastmath)
-        self.check()
+        self.__check()
 
     def __str__(self):
         description = []
@@ -128,7 +128,7 @@ def __str__(self):
                 description.append(f"{attr}: {value}")
         return ', '.join(description)
 
-    def check(self):
+    def __check(self):
         for attr in dir(self):
             if not attr.startswith('__'):
                 if hasattr(getattr(self, attr), '_check'):

PySDM/physics/heterogeneous_ice_nucleation_rate/__init__.py

@@ -1,3 +1,3 @@
 from .constant import Constant
 from .null import Null
-from .ABIFM import ABIFM
+from .abifm import ABIFM

PySDM/physics/heterogeneous_ice_nucleation_rate/abifm.py

@@ -0,0 +1,14 @@
+import numpy as np
+
+m = np.inf
+c = np.inf
+
+class ABIFM:
+    @staticmethod
+    def _check():
+        assert np.isfinite(m)
+        assert np.isfinite(c)
+
+    @staticmethod
+    def J_het(a_w_ice):
+        return 10**(m * (1 - a_w_ice) + c)

PySDM/physics/heterogeneous_ice_nucleation_rate/constant.py

@@ -8,5 +8,5 @@ def _check():
         assert np.isfinite(J_het)
 
     @staticmethod
-    def J_het():
+    def J_het(a_w_ice):
         return J_het

PySDM/physics/heterogeneous_ice_nucleation_rate/null.py

@@ -1,4 +1,4 @@
 class Null:
     @staticmethod
-    def J_het():
+    def J_het(a_w_ice):
         pass

PySDM/physics/saturation_vapour_pressure/flatau_walko_cotton.py

@@ -29,3 +29,11 @@ def ice_Celsius(T):
                 const.FWC_I7 + T *
                 const.FWC_I8
         ))))))))
+
+    @staticmethod
+    def a_w_ice(T):
+        return (
+            FlatauWalkoCotton.ice_Celsius(T - const.T0)
+            /
+            FlatauWalkoCotton.pvs_Celsius(T - const.T0)
+        )

test-time-requirements.txt

@@ -4,4 +4,4 @@ ghapi
 pytest
 
 # note: if cloning both PySDM and PySDM examples, consider "pip install -e"
-PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@0be0f1e#egg=PySDM-examples
+PySDM-examples @ git+git://github.com/slayoo/PySDM-examples@b37041a#egg=PySDM-examples

tests/unit_tests/backends/test_freezing_methods.py

@@ -57,7 +57,8 @@ def test_freeze_time_dependent(plot=False):
 
             formulae = Formulae(heterogeneous_ice_nucleation_rate='Constant')
             builder = Builder(n_sd=n_sd, backend=CPU, formulae=formulae)
-            builder.set_environment(Box(dt=case['dt'], dv=dv))
+            env = Box(dt=case['dt'], dv=dv)
+            builder.set_environment(env)
             builder.add_dynamic(Freezing(singular=False))
             attributes = {
                 'n': np.full(n_sd, int(case['N'])),
@@ -67,6 +68,8 @@ def test_freeze_time_dependent(plot=False):
             products = (IceWaterContent(specific=False),)
             particulator = builder.build(attributes=attributes, products=products)
 
+            env['a_w_ice'] = np.nan
+
             cell_id = 0
             for i in range(int(total_time / case['dt']) + 1):
                 particulator.run(0 if i == 0 else 1)