vapour deposition on ice (temporarily as a separate dynamic - to be combined with condensation); introducing physics.latent_heat_sublimation (and renamings in physics: latent_heat to latent_heat_vapourisation, condensation_coordinate to diffusion_coordinate) (#1389)

Co-authored-by: Tim Luettmer <[email protected]>

Author: slayoo

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/condensation_methods.py

@@ -290,7 +290,7 @@ def step_impl(  # pylint: disable=too-many-arguments,too-many-locals
                     rhod, T, water_vapour_mixing_ratio
                 )
                 pv = formulae.state_variable_triplet__pv(p, water_vapour_mixing_ratio)
-                lv = formulae.latent_heat__lv(T)
+                lv = formulae.latent_heat_vapourisation__lv(T)
                 pvs = formulae.saturation_vapour_pressure__pvs_water(T)
                 DTp = formulae.diffusion_thermics__D(T, p)
                 RH = pv / pvs
@@ -390,9 +390,9 @@ def minfun(  # pylint: disable=too-many-arguments,too-many-locals
             K,
             ventilation_factor,
         ):
-            if x_new > formulae.condensation_coordinate__x_max():
+            if x_new > formulae.diffusion_coordinate__x_max():
                 return x_old - x_new
-            mass_new = formulae.condensation_coordinate__mass(x_new)
+            mass_new = formulae.diffusion_coordinate__mass(x_new)
             volume_new = formulae.particle_shape_and_density__mass_to_volume(mass_new)
             r_new = formulae.trivia__radius(volume_new)
             RH_eq = formulae.hygroscopicity__RH_eq(
@@ -418,7 +418,7 @@ def minfun(  # pylint: disable=too-many-arguments,too-many-locals
             return (
                 x_old
                 - x_new
-                + timestep * formulae.condensation_coordinate__dx_dt(mass_new, dm_dt)
+                + timestep * formulae.diffusion_coordinate__dx_dt(mass_new, dm_dt)
             )
 
         @numba.njit(**jit_flags)
@@ -450,9 +450,9 @@ def calculate_ml_new(  # pylint: disable=too-many-branches,too-many-arguments,to
                 v_drop = formulae.particle_shape_and_density__mass_to_volume(
                     attributes.water_mass[drop]
                 )
-                x_old = formulae.condensation_coordinate__x(attributes.water_mass[drop])
+                x_old = formulae.diffusion_coordinate__x(attributes.water_mass[drop])
                 r_old = formulae.trivia__radius(v_drop)
-                x_insane = formulae.condensation_coordinate__x(
+                x_insane = formulae.diffusion_coordinate__x(
                     formulae.particle_shape_and_density__volume_to_mass(
                         attributes.vdry[drop] / 100
                     )
@@ -499,11 +499,11 @@ def calculate_ml_new(  # pylint: disable=too-many-branches,too-many-arguments,to
                         Kr,
                         ventilation_factor,
                     )
-                    mass_old = formulae.condensation_coordinate__mass(x_old)
+                    mass_old = formulae.diffusion_coordinate__mass(x_old)
                     dm_dt_old = formulae.particle_shape_and_density__dm_dt(
                         r=r_old, r_dr_dt=r_dr_dt_old
                     )
-                    dx_old = timestep * formulae.condensation_coordinate__dx_dt(
+                    dx_old = timestep * formulae.diffusion_coordinate__dx_dt(
                         mass_old, dm_dt_old
                     )
                 else:
@@ -572,7 +572,7 @@ def calculate_ml_new(  # pylint: disable=too-many-branches,too-many-arguments,to
                     else:
                         x_new = x_old
 
-                mass_new = formulae.condensation_coordinate__mass(x_new)
+                mass_new = formulae.diffusion_coordinate__mass(x_new)
                 mass_cr = formulae.particle_shape_and_density__volume_to_mass(
                     attributes.v_cr[drop]
                 )

PySDM/backends/impl_numba/methods/deposition_methods.py

@@ -0,0 +1,172 @@
+"""basic water vapor deposition on ice for CPU backend, for Howell factor see:
+[Howell 1949](https://doi.org/10.1175/1520-0469(1949)006%3C0134:TGOCDI%3E2.0.CO;2)
+"""
+
+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
+
+        @numba.njit(**{**self.default_jit_flags, **{"parallel": False}})
+        def body(  # pylint: disable=too-many-arguments
+            *,
+            multiplicity,
+            signed_water_mass,
+            current_temperature,
+            current_total_pressure,
+            current_relative_humidity,
+            current_water_activity,
+            current_vapour_mixing_ratio,
+            current_dry_air_density,
+            current_dry_potential_temperature,
+            cell_volume,
+            time_step,
+            cell_id,
+            reynolds_number,
+            schmidt_number,
+            # to be modified
+            predicted_vapour_mixing_ratio,
+            predicted_dry_potential_temperature,
+        ):
+            # pylint: disable=too-many-locals
+            n_sd = len(signed_water_mass)
+            for i in range(n_sd):
+                if not formulae.trivia__unfrozen(signed_water_mass[i]):
+                    ice_mass = -signed_water_mass[i]
+                    cid = cell_id[i]
+
+                    radius = formulae.particle_shape_and_density__mass_to_radius(
+                        signed_water_mass[i]
+                    )
+
+                    diameter = radius * 2.0
+
+                    temperature = current_temperature[cid]
+                    pressure = current_total_pressure[cid]
+                    rho = current_dry_air_density[cid]
+                    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
+                    )
+                    saturation_ratio_ice = (
+                        current_relative_humidity[cid] / current_water_activity[cid]
+                    )
+
+                    if saturation_ratio_ice == 1:
+                        continue
+
+                    howell_factor_x_diffcoef_x_rhovsice_x_icess = (
+                        formulae.drop_growth__r_dr_dt(
+                            RH_eq=1,
+                            T=temperature,
+                            RH=saturation_ratio_ice,
+                            lv=latent_heat_sub,
+                            pvs=pvs_ice,
+                            D=diffusion_coefficient,
+                            K=thermal_conductivity,
+                            ventilation_factor=ventilation_factor,
+                        )
+                        * formulae.constants.rho_w
+                    )
+
+                    dm_dt = (
+                        4
+                        * np.pi
+                        * capacity
+                        * howell_factor_x_diffcoef_x_rhovsice_x_icess
+                    )
+
+                    delta_rv_i = (
+                        -dm_dt * multiplicity[i] * time_step / (cell_volume * rho)
+                    )
+                    if -delta_rv_i > current_vapour_mixing_ratio[cid]:
+                        assert False
+                    predicted_vapour_mixing_ratio[cid] += delta_rv_i
+
+                    predicted_dry_potential_temperature[cid] += (
+                        formulae.state_variable_triplet__dthd_dt(
+                            rhod=current_dry_air_density[cid],
+                            thd=current_dry_potential_temperature[cid],
+                            T=temperature,
+                            d_water_vapour_mixing_ratio__dt=delta_rv_i / time_step,
+                            lv=latent_heat_sub,
+                        )
+                        * time_step
+                    )
+
+                    x_old = formulae.diffusion_coordinate__x(ice_mass)
+                    dx_dt_old = formulae.diffusion_coordinate__dx_dt(ice_mass, dm_dt)
+                    x_new = formulae.trivia__explicit_euler(x_old, time_step, dx_dt_old)
+                    signed_water_mass[i] = -formulae.diffusion_coordinate__mass(x_new)
+
+        return body
+
+    def deposition(  # pylint: disable=too-many-locals
+        self,
+        *,
+        multiplicity,
+        signed_water_mass,
+        current_temperature,
+        current_total_pressure,
+        current_relative_humidity,
+        current_water_activity,
+        current_vapour_mixing_ratio,
+        current_dry_air_density,
+        current_dry_potential_temperature,
+        cell_volume,
+        time_step,
+        cell_id,
+        reynolds_number,
+        schmidt_number,
+        predicted_vapour_mixing_ratio,
+        predicted_dry_potential_temperature,
+    ):
+        self._deposition(
+            multiplicity=multiplicity.data,
+            signed_water_mass=signed_water_mass.data,
+            current_temperature=current_temperature.data,
+            current_total_pressure=current_total_pressure.data,
+            current_relative_humidity=current_relative_humidity.data,
+            current_water_activity=current_water_activity.data,
+            current_vapour_mixing_ratio=current_vapour_mixing_ratio.data,
+            current_dry_air_density=current_dry_air_density.data,
+            current_dry_potential_temperature=current_dry_potential_temperature.data,
+            cell_volume=cell_volume,
+            time_step=time_step,
+            cell_id=cell_id.data,
+            reynolds_number=reynolds_number.data,
+            schmidt_number=schmidt_number.data,
+            predicted_vapour_mixing_ratio=predicted_vapour_mixing_ratio.data,
+            predicted_dry_potential_temperature=predicted_dry_potential_temperature.data,
+        )

PySDM/backends/impl_numba/test_helpers/scipy_ode_condensation_solver.py

@@ -92,7 +92,7 @@ def _make_solve(formulae):  # pylint: disable=too-many-statements,too-many-local
 
     @numba.njit(**{**JIT_FLAGS, **{"parallel": False}})
     def _liquid_water_mixing_ratio(n, x, m_d_mean):
-        return np.sum(n * jit_formulae.condensation_coordinate__mass(x)) / m_d_mean
+        return np.sum(n * jit_formulae.diffusion_coordinate__mass(x)) / m_d_mean
 
     @numba.njit(**{**JIT_FLAGS, **{"parallel": False}})
     def _impl(  # pylint: disable=too-many-arguments,too-many-locals
@@ -127,7 +127,7 @@ def _impl(  # pylint: disable=too-many-arguments,too-many-locals
         )
         sum_n_dm_dt = 0
         for i, x_i in enumerate(x):
-            m = jit_formulae.condensation_coordinate__mass(x_i)
+            m = jit_formulae.diffusion_coordinate__mass(x_i)
             v = jit_formulae.particle_shape_and_density__mass_to_volume(m)
             r = jit_formulae.trivia__radius(v)
             Dr = jit_formulae.diffusion_kinetics__D(DTp, r, lambdaD)
@@ -152,7 +152,7 @@ def _impl(  # pylint: disable=too-many-arguments,too-many-locals
                 ventilation_factor,
             )
             dm_dt = jit_formulae.particle_shape_and_density__dm_dt(r, r_dr_dt)
-            dy_dt[idx_x + i] = jit_formulae.condensation_coordinate__dx_dt(m, dm_dt)
+            dy_dt[idx_x + i] = jit_formulae.diffusion_coordinate__dx_dt(m, dm_dt)
             sum_n_dm_dt += n[i] * dm_dt
         dy_dt[idx_thd] = dot_thd + jit_formulae.state_variable_triplet__dthd_dt(
             rhod, thd, T, dot_water_vapour_mixing_ratio - sum_n_dm_dt / m_d_mean, lv
@@ -211,7 +211,7 @@ def _odesys(  # pylint: disable=too-many-arguments,too-many-locals
             air_dynamic_viscosity,
             rhod,
             pvs,
-            jit_formulae.latent_heat__lv(T),
+            jit_formulae.latent_heat_vapourisation__lv(T),
         )
         return dy_dt
 
@@ -234,7 +234,7 @@ def solve(  # pylint: disable=too-many-arguments,too-many-locals,unused-argument
         n_sd_in_cell = len(cell_idx)
         y0 = np.empty(n_sd_in_cell + idx_x)
         y0[idx_thd] = thd
-        y0[idx_x:] = jit_formulae.condensation_coordinate__x(
+        y0[idx_x:] = jit_formulae.diffusion_coordinate__x(
             attributes.water_mass[cell_idx]
         )
         total_water_mixing_ratio = (
@@ -282,7 +282,7 @@ def solve(  # pylint: disable=too-many-arguments,too-many-locals,unused-argument
         m_new = 0
         for i in range(n_sd_in_cell):
             attributes.water_mass[cell_idx[i]] = (
-                jit_formulae.condensation_coordinate__mass(y1[idx_x + i])
+                jit_formulae.diffusion_coordinate__mass(y1[idx_x + i])
             )
             m_new += (
                 attributes.multiplicity[cell_idx[i]]

PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py

@@ -100,10 +100,10 @@ def __update_drop_masses(self):
             struct Minfun {{
                 static __device__ real_type value(real_type x_new, void* args_p) {{
                     auto args = static_cast<real_type*>(args_p);
-                    if (x_new > {phys.condensation_coordinate.x_max.c_inline()}) {{
+                    if (x_new > {phys.diffusion_coordinate.x_max.c_inline()}) {{
                         return {args("x_old")} - x_new;
                     }}
-                    auto m_new = {phys.condensation_coordinate.mass.c_inline(x="x_new")};
+                    auto m_new = {phys.diffusion_coordinate.mass.c_inline(x="x_new")};
                     auto v_new = {phys.particle_shape_and_density.mass_to_volume.c_inline(mass="m_new")};
                     auto r_new = {phys.trivia.radius.c_inline(volume="v_new")};
                     auto sgm = {phys.surface_tension.sigma.c_inline(
@@ -133,7 +133,7 @@ def __update_drop_masses(self):
                         r="r_new", r_dr_dt="r_dr_dt"
                     )};
                     return {args("x_old")} - x_new + {args("dt")} * {
-                        phys.condensation_coordinate.dx_dt.c_inline(m="m_new", dm_dt="dm_dt")
+                        phys.diffusion_coordinate.dx_dt.c_inline(m="m_new", dm_dt="dm_dt")
                     };
                 }}
             }};
@@ -152,10 +152,10 @@ def __update_drop_masses(self):
             auto v_old = {phys.particle_shape_and_density.mass_to_volume.c_inline(
                 mass="water_mass[i]"
             )};
-            auto x_old = {phys.condensation_coordinate.x.c_inline(mass="water_mass[i]")};
+            auto x_old = {phys.diffusion_coordinate.x.c_inline(mass="water_mass[i]")};
             auto r_old = {phys.trivia.radius.c_inline(volume="v_old")};
             auto m_insane = {phys.particle_shape_and_density.volume_to_mass.c_inline(volume="vdry[i] / 100")};
-            auto x_insane = {phys.condensation_coordinate.x.c_inline(mass="m_insane")};
+            auto x_insane = {phys.diffusion_coordinate.x.c_inline(mass="m_insane")};
             auto rd3 = vdry[i] / {const.PI_4_3};
             auto sgm = {phys.surface_tension.sigma.c_inline(
                 T="_T", v_wet="v", v_dry="vdry[i]", f_org="_f_org[i]"
@@ -194,7 +194,7 @@ def __update_drop_masses(self):
                     ventilation_factor="ventilation_factor",
                 )};
                 dm_dt_old = {phys.particle_shape_and_density.dm_dt.c_inline(r="r_old", r_dr_dt="r_dr_dt_old")};
-                dx_old = dt * {phys.condensation_coordinate.dx_dt.c_inline(
+                dx_old = dt * {phys.diffusion_coordinate.dx_dt.c_inline(
                     m="water_mass[i]", dm_dt="dm_dt_old"
                 )};
             }}
@@ -249,7 +249,7 @@ def __update_drop_masses(self):
                     x_new = x_old;
                 }}
             }}
-            water_mass[i] = {phys.condensation_coordinate.mass.c_inline(x="x_new")};
+            water_mass[i] = {phys.diffusion_coordinate.mass.c_inline(x="x_new")};
         """.replace(
                 "real_type", self._get_c_type()
             ),
@@ -323,7 +323,7 @@ def __pre(self):
             )};
             pv[i] = {phys.state_variable_triplet.pv.c_inline(
                 p='p[i]', water_vapour_mixing_ratio='predicted_water_vapour_mixing_ratio[i]')};
-            lv[i] = {phys.latent_heat.lv.c_inline(
+            lv[i] = {phys.latent_heat_vapourisation.lv.c_inline(
                 T='T[i]')};
             pvs[i] = {phys.saturation_vapour_pressure.pvs_water.c_inline(
                 T='T[i]')};

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
@@ -78,3 +79,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()