work in progress on mass-based diffusion coordinate for deposition

PySDM/backends/impl_numba/methods/deposition_methods.py

@@ -11,13 +11,19 @@ def _deposition(self):
         assert self.formulae.particle_shape_and_density.supports_mixed_phase()
 
         formulae = self.formulae_flattened
-        liquid = formulae.trivia__unzfroze
+        liquid = formulae.trivia__unfrozen
 
         diffusion_coefficient_function = self.formulae.diffusion_thermics.D
 
         @numba.jit(**self.default_jit_flags)
         def body(
-            water_mass, ambient_temperature, ambient_total_pressure, time_step, cell_id
+            water_mass,
+            ambient_temperature,
+            ambient_total_pressure,
+            time_step,
+            cell_id,
+            reynolds_number,
+            schmidt_number,
         ):
 
             n_sd = len(water_mass)
@@ -32,6 +38,16 @@ def body(
 
                     temperature = ambient_temperature[cid]
                     pressure = ambient_total_pressure[cid]
+                    capacity = formulae.particle_shape_and_density__ice_mass_to_radius(
+                        water_mass[i]
+                    )
+
+                    ventilation_factor = formulae.ventilation__ventilation_coefficient(
+                        sqrt_re_times_cbrt_sc=formulae.trivia__sqrt_re_times_cbrt_sc(
+                            Re=reynolds_number[i],
+                            Sc=schmidt_number[cid],
+                        )
+                    )
 
                     diffusion_coefficient = diffusion_coefficient_function(
                         temperature, pressure
@@ -41,7 +57,13 @@ def body(
                         volume, temperature, pressure, diffusion_coefficient, time_step
                     )
 
-                    water_mass[i] *= 1.1
+                    dm_dt = 1.1
+
+                    x_old = formulae.diffusion_coordinate__x(water_mass[i])
+                    dx_dt_old = formulae.diffusion_coordinate__dx_dt(x_old, dm_dt)
+                    x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old)
+
+                    water_mass[i] = formulae.diffusion_coordinate__mass(x_new)
 
         return body
 
@@ -52,12 +74,16 @@ def deposition(
         ambient_temperature,
         ambient_total_pressure,
         time_step,
-        cell_id
+        cell_id,
+        reynolds_number,
+        schmidt_number,
     ):
         self._deposition(
             water_mass=water_mass.data,
             ambient_temperature=ambient_temperature.data,
             ambient_total_pressure=ambient_total_pressure.data,
             time_step=time_step,
             cell_id=cell_id.data,
+            reynolds_number=reynolds_number.data,
+            schmidt_number=schmidt_number.data,
         )

PySDM/dynamics/vapour_deposition_on_ice.py

@@ -12,6 +12,7 @@ def __init__(self):
     def register(self, *, builder):
         """called by the builder"""
         self.particulator = builder.particulator
+        builder.request_attribute("Reynolds number")
 
     def __call__(self):
         """called by the particulator during simulation"""

PySDM/particulator.py

@@ -487,5 +487,7 @@ def deposition(self):
             ambient_total_pressure=self.environment["P"],
             time_step=self.dt,
             cell_id=self.attributes["cell id"],
+            reynolds_number=self.attributes["Reynolds number"],
+            schmidt_number=self.environment["Schmidt number"],
         )
         self.attributes.mark_updated("water mass")

PySDM/physics/diffusion_coordinate/__init__.py

@@ -4,3 +4,5 @@
 
 from .volume import Volume
 from .volume_logarithm import VolumeLogarithm
+from .mass import Mass
+from .mass_logarithm import MassLogarithm

PySDM/physics/diffusion_coordinate/mass.py

@@ -0,0 +1,22 @@
+"""
+particle mass as diffusion coordinate (i.e. no transformation)
+"""
+
+import numpy as np
+
+
+class Mass:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def dx_dt(const, x, dm_dt):
+        return dm_dt
+
+    @staticmethod
+    def mass(_, x):
+        return x
+
+    @staticmethod
+    def x(_, mass):
+        return mass

PySDM/physics/diffusion_coordinate/mass_logarithm.py

@@ -0,0 +1,23 @@
+"""
+logarithm of particle mass as diffusion coordinate
+(ensures non-negative values)
+"""
+
+import numpy as np
+
+
+class MassLogarithm:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def dx_dt(const, x, dm_dt):
+        return dm_dt / np.exp(x)
+
+    @staticmethod
+    def mass(x):
+        return np.exp(x)
+
+    @staticmethod
+    def x(mass):
+        return np.log(mass)

PySDM/physics/particle_shape_and_density/mixed_phase_spheres.py

@@ -30,3 +30,7 @@ def volume_to_mass(const, volume):
     @staticmethod
     def radius_to_mass(const, radius):
         raise NotImplementedError()
+
+    @staticmethod
+    def ice_mass_to_radius(const, ice_mass):
+        return np.power(-ice_mass / const.PI_4_3 / const.rho_i, const.ONE_THIRD)

tests/unit_tests/dynamics/test_vapour_deposition_on_ice.py

@@ -4,7 +4,7 @@
 
 import pytest
 
-from PySDM.physics import si
+from PySDM.physics import si, diffusion_coordinate
 from PySDM.backends import CPU
 from PySDM import Builder
 from PySDM import Formulae
@@ -16,15 +16,20 @@
 @pytest.mark.parametrize("dt", (1 * si.s, 0.1 * si.s))
 @pytest.mark.parametrize("water_mass", (-si.ng, -si.ug, -si.mg, si.mg))
 @pytest.mark.parametrize("fastmath", (True, False))
-def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3):
+@pytest.mark.parametrize("diffusion_coordinate", ("Mass", "MassLogarithm"))
+def test_iwc_lower_after_timestep(
+    dt, water_mass, fastmath, diffusion_coordinate, dv=1 * si.m**3
+):
     # arrange
     n_sd = 1
     builder = Builder(
         n_sd=n_sd,
         environment=Box(dt=dt, dv=dv),
         backend=CPU(
             formulae=Formulae(
-                fastmath=fastmath, particle_shape_and_density="MixedPhaseSpheres"
+                fastmath=fastmath,
+                particle_shape_and_density="MixedPhaseSpheres",
+                diffusion_coordinate=diffusion_coordinate,
             )
         ),
     )
@@ -39,6 +44,8 @@ def test_iwc_lower_after_timestep(dt, water_mass, fastmath, dv=1 * si.m**3):
     )
     particulator.environment["T"] = 250 * si.K
     particulator.environment["P"] = 500 * si.hPa
+    particulator.environment["Schmidt number"] = 1
+
     # act
     iwc_old = particulator.products["ice water content"].get().copy()
     particulator.run(steps=1)