new dynamic: vapour deposition on ice

PySDM/backends/impl_numba/methods/__init__.py

@@ -13,3 +13,4 @@
 from .physics_methods import PhysicsMethods
 from .terminal_velocity_methods import TerminalVelocityMethods
 from .seeding_methods import SeedingMethods
+from .deposition_methods import DepositionMethods

PySDM/backends/impl_numba/methods/deposition_methods.py

@@ -0,0 +1,147 @@
+""" basic water vapor deposition on ice for cpu backend """
+
+from functools import cached_property
+import numba
+import numpy as np
+from PySDM.backends.impl_common.backend_methods import BackendMethods
+
+
+class DepositionMethods(BackendMethods):  # pylint:disable=too-few-public-methods
+    @cached_property
+    def _deposition(self):
+        assert self.formulae.particle_shape_and_density.supports_mixed_phase()
+
+        formulae = self.formulae_flattened
+        liquid = formulae.trivia__unfrozen
+
+
+        @numba.jit(**self.default_jit_flags)
+        def body(
+            multiplicity,
+            water_mass,
+            ambient_temperature,
+            ambient_total_pressure,
+            ambient_humidity,
+            ambient_water_activity,
+            ambient_vapour_mixing_ratio,
+            ambient_dry_air_density,
+            cell_volume,
+            time_step,
+            cell_id,
+            reynolds_number,
+            schmidt_number,
+        ):
+            n_sd = len(water_mass)
+            ## for i in numba.prange(n_sd):  # pylint: disable=not-an-iterable
+            for i in range(n_sd):
+                if not liquid(water_mass[i]):
+                    ice_mass = -water_mass[i]
+                    cid = cell_id[i]
+
+                    radius = formulae.particle_shape_and_density__ice_mass_to_radius(
+                        water_mass[i]
+                    )
+                    diameter = radius * 2.
+
+                    temperature = ambient_temperature[cid]
+                    pressure = ambient_total_pressure[cid]
+                    rho = ambient_dry_air_density[cid]
+                    Rv = formulae.constants.Rv
+                    pvs_ice =  formulae.saturation_vapour_pressure__pvs_ice(temperature)
+                    latent_heat_sub = formulae.latent_heat_sublimation__ls(temperature)
+
+                    capacity = formulae.diffusion_ice_capacity__capacity( diameter )
+
+                    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],
+                        )
+                    )
+
+                    Dv_const = formulae.diffusion_thermics__D(temperature, pressure)
+                    lambdaD = formulae.diffusion_ice_kinetics__lambdaD(temperature, pressure)
+                    diffusion_coefficient = formulae.diffusion_ice_kinetics__D(Dv_const, radius, lambdaD, temperature)
+
+                    Ka_const = formulae.diffusion_thermics__K(temperature, pressure)
+                    lambdaK = formulae.diffusion_ice_kinetics__lambdaK(temperature, pressure)
+                    thermal_conductivity = formulae.diffusion_ice_kinetics__K(Ka_const,  radius, lambdaK, temperature, rho)  
+
+
+                    howell_factor = 1. / ((latent_heat_sub / Rv / temperature - 1.) * latent_heat_sub * diffusion_coefficient / temperature / thermal_conductivity +  Rv * temperature  /  pvs_ice)
+                    # print( f" {howell_factor=}, {latent_heat_sub= }, {Rv=}, {Dv_const=}, {diffusion_coefficient=}, {Ka_const=}, {thermal_conductivity=}, {pvs_ice=}")
+
+                    saturation_ratio_ice = (
+                        ambient_humidity[cid] / ambient_water_activity[cid]
+                    )
+
+                    rho_vs_ice = (
+                        pvs_ice
+                        / Rv
+                        / temperature
+                    )
+
+                    dm_dt = (
+                        4
+                        * np.pi
+                        * capacity
+                        * diffusion_coefficient
+                        * (saturation_ratio_ice - 1)
+                    ) * rho_vs_ice
+
+                    if dm_dt == 0:
+                        continue
+
+                    delta_rv_i = (
+                        -dm_dt
+                        * multiplicity[i]
+                        * time_step
+                        / (cell_volume * rho)
+                    )
+                    if -delta_rv_i > ambient_vapour_mixing_ratio[cid]:
+                        assert False
+                    ambient_vapour_mixing_ratio[cid] += delta_rv_i
+
+                    delta_T = -delta_rv_i * latent_heat_sub / formulae.constants.c_pd
+                    ambient_temperature[cid] += delta_T 
+
+                    x_old = formulae.diffusion_coordinate__x(ice_mass)
+                    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
+
+    def deposition(
+        self,
+        *,
+        multiplicity,
+        water_mass,
+        ambient_temperature,
+        ambient_total_pressure,
+        ambient_humidity,
+        ambient_water_activity,
+        ambient_vapour_mixing_ratio,
+        ambient_dry_air_density,
+        cell_volume,
+        time_step,
+        cell_id,
+        reynolds_number,
+        schmidt_number,
+    ):
+        self._deposition(
+            multiplicity=multiplicity.data,
+            water_mass=water_mass.data,
+            ambient_temperature=ambient_temperature.data,
+            ambient_total_pressure=ambient_total_pressure.data,
+            ambient_humidity=ambient_humidity.data,
+            ambient_water_activity=ambient_water_activity.data,
+            ambient_vapour_mixing_ratio=ambient_vapour_mixing_ratio.data,
+            ambient_dry_air_density=ambient_dry_air_density.data,
+            cell_volume=cell_volume,
+            time_step=time_step,
+            cell_id=cell_id.data,
+            reynolds_number=reynolds_number.data,
+            schmidt_number=schmidt_number.data,
+        )

PySDM/backends/numba.py

@@ -29,6 +29,7 @@ class Numba(  # pylint: disable=too-many-ancestors,duplicate-code
     methods.TerminalVelocityMethods,
     methods.IsotopeMethods,
     methods.SeedingMethods,
+    methods.DepositionMethods,
 ):
     Storage = ImportedStorage
     Random = ImportedRandom
@@ -77,3 +78,4 @@ def __init__(self, formulae=None, double_precision=True, override_jit_flags=None
         methods.TerminalVelocityMethods.__init__(self)
         methods.IsotopeMethods.__init__(self)
         methods.SeedingMethods.__init__(self)
+        methods.DepositionMethods.__init__(self)

PySDM/dynamics/__init__.py

@@ -16,3 +16,4 @@
 from PySDM.dynamics.freezing import Freezing
 from PySDM.dynamics.relaxed_velocity import RelaxedVelocity
 from PySDM.dynamics.seeding import Seeding
+from PySDM.dynamics.vapour_deposition_on_ice import VapourDepositionOnIce

PySDM/dynamics/vapour_deposition_on_ice.py

@@ -0,0 +1,20 @@
+""" basic water vapor deposition on ice """
+
+from PySDM.dynamics.impl import register_dynamic
+
+
+@register_dynamic()
+class VapourDepositionOnIce:
+    def __init__(self):
+        """called by the user while building a particulator"""
+        self.particulator = None
+
+    def register(self, *, builder):
+        """called by the builder"""
+        self.particulator = builder.particulator
+        assert builder.formulae.particle_shape_and_density.supports_mixed_phase()
+        builder.request_attribute("Reynolds number")
+
+    def __call__(self):
+        """called by the particulator during simulation"""
+        self.particulator.deposition()

PySDM/formulae.py

@@ -30,13 +30,17 @@ def __init__(  # pylint: disable=too-many-locals
         constants: Optional[dict] = None,
         seed: int = None,
         fastmath: bool = True,
+        diffusion_coordinate: str = "MassLogarithm",
         condensation_coordinate: str = "VolumeLogarithm",
         saturation_vapour_pressure: str = "FlatauWalkoCotton",
         latent_heat: str = "Kirchhoff",
+        latent_heat_sublimation: str = "MurphyKoop",
         hygroscopicity: str = "KappaKoehlerLeadingTerms",
         drop_growth: str = "Mason1971",
         surface_tension: str = "Constant",
         diffusion_kinetics: str = "FuchsSutugin",
+        diffusion_ice_kinetics: str = "Standard",
+        diffusion_ice_capacity: str = "Spherical",
         diffusion_thermics: str = "Neglect",
         ventilation: str = "Neglect",
         state_variable_triplet: str = "LibcloudphPlusPlus",
@@ -62,13 +66,17 @@ def __init__(  # pylint: disable=too-many-locals
         # in PyCharm and alike, all these fields are later overwritten within this ctor
         self.optical_albedo = optical_albedo
         self.optical_depth = optical_depth
+        self.diffusion_coordinate = diffusion_coordinate
         self.condensation_coordinate = condensation_coordinate
         self.saturation_vapour_pressure = saturation_vapour_pressure
         self.hygroscopicity = hygroscopicity
         self.drop_growth = drop_growth
         self.surface_tension = surface_tension
         self.diffusion_kinetics = diffusion_kinetics
+        self.diffusion_ice_kinetics = diffusion_ice_kinetics
+        self.diffusion_ice_capacity = diffusion_ice_capacity
         self.latent_heat = latent_heat
+        self.latent_heat_sublimation = latent_heat_sublimation
         self.diffusion_thermics = diffusion_thermics
         self.ventilation = ventilation
         self.state_variable_triplet = state_variable_triplet

PySDM/particulator.py

@@ -479,3 +479,21 @@ def seeding(
         self.attributes.mark_updated("multiplicity")
         for key in self.attributes.get_extensive_attribute_keys():
             self.attributes.mark_updated(key)
+
+    def deposition(self):
+        self.backend.deposition(
+            multiplicity=self.attributes["multiplicity"],
+            water_mass=self.attributes["water mass"],
+            ambient_temperature=self.environment["T"],
+            ambient_total_pressure=self.environment["P"],
+            ambient_humidity=self.environment["RH"],
+            ambient_water_activity=self.environment["a_w_ice"],
+            ambient_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"],
+            ambient_dry_air_density=self.environment["rhod"],
+            cell_volume=self.environment.mesh.dv,
+            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/__init__.py

@@ -17,9 +17,12 @@
 """
 
 from . import (
+    diffusion_coordinate,
     condensation_coordinate,
     constants_defaults,
     diffusion_kinetics,
+    diffusion_ice_kinetics,
+    diffusion_ice_capacity,
     diffusion_thermics,
     drop_growth,
     fragmentation_function,
@@ -34,6 +37,7 @@
     isotope_diffusivity_ratios,
     isotope_relaxation_timescale,
     latent_heat,
+    latent_heat_sublimation,
     optical_albedo,
     optical_depth,
     particle_advection,

PySDM/physics/condensation_coordinate/__init__.py

@@ -1,6 +1,2 @@
-"""
-definitions of particle-size coordinates for the condensation solver
-"""
-
 from .volume import Volume
 from .volume_logarithm import VolumeLogarithm

PySDM/physics/condensation_coordinate/volume.py

@@ -1,5 +1,5 @@
 """
-particle volume as condensation coordinate (i.e. no transformation)
+particle volume as diffusion coordinate (i.e. no transformation)
 """
 
 import numpy as np

PySDM/physics/condensation_coordinate/volume_logarithm.py

@@ -1,5 +1,6 @@
 """
-logarithm of particle volume as coordinate (ensures non-negative values)
+logarithm of particle volume as diffusion coordinate
+(ensures non-negative values)
 """
 
 import numpy as np

PySDM/physics/constants_defaults.py

@@ -62,10 +62,14 @@
 
 K0 = 2.4e-2 * si.joules / si.metres / si.seconds / si.kelvins
 
-# mass and heat accommodation coefficients
+# mass and heat accommodation coefficients for condensation
 MAC = 1.0
 HAC = 1.0
 
+# mass and heat accommodation coefficients for vapour deposition on ice
+MAC_ice = 0.5
+HAC_ice = 0.7
+
 p1000 = 1000 * si.hectopascals
 c_pd = 1005 * si.joule / si.kilogram / si.kelvin
 c_pv = 1850 * si.joule / si.kilogram / si.kelvin
@@ -127,6 +131,7 @@
 l_l19_a = 0.167 * si.dimensionless
 l_l19_b = 3.65e-4 / si.kelvin
 
+
 # Thermal diffusivity constants from Lowe et al. (2019)
 k_l19_a = 4.2e-3 * si.joules / si.metres / si.seconds / si.kelvins
 k_l19_b = 1.0456 * si.dimensionless
@@ -135,6 +140,10 @@
 # Delta v for diffusivity in Pruppacher & Klett eq. 13-14
 dv_pk05 = 0.0 * si.metres
 
+# mean free path of lambda of water molecules
+# in pruppacher und klett 
+lmbd_w_0 = 6.6e-8 * si.metre
+
 # Seinfeld & Pandis eq. 15.65 (Hall & Pruppacher 1976)
 d_l19_a = 0.211e-4 * si.metre**2 / si.second
 d_l19_b = 1.94
@@ -160,6 +169,13 @@
 MK05_LIQ_C12 = 1 * si.K
 MK05_LIQ_C13 = 0.014025 / si.K
 
+# Latent heat of sublimation from Murphy and Koop (2005)
+MK05_SUB_C1 = 46782.5 * si.joule /  si.mole
+MK05_SUB_C2 = 35.8925 * si.joule /  si.mole / si.kelvin
+MK05_SUB_C3 = 0.07414 * si.joule /  si.mole / si.kelvin**2
+MK05_SUB_C4 = 541.5 * si.joule /  si.mole
+MK05_SUB_C5 = 123.75 * si.kelvin
+
 # standard pressure and temperature (ICAO)
 T_STP = (sci.zero_Celsius + 15) * si.kelvin
 p_STP = 101325 * si.pascal

PySDM/physics/diffusion_coordinate/__init__.py

@@ -0,0 +1,6 @@
+"""
+definitions of particle-size coordinates for the vapour diffusion solvers
+"""
+
+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)