terminal velocity of ice particles

PySDM/attributes/physics/relative_fall_velocity.py

@@ -1,5 +1,5 @@
 """
-Attributes for tracking droplet velocity
+Attributes for tracking particle velocity
 """
 
 from PySDM.attributes.impl import (
@@ -29,13 +29,13 @@ def __init__(self, builder):
 class RelativeFallVelocity(DerivedAttribute):
     def __init__(self, builder):
         self.momentum = builder.get_attribute("relative fall momentum")
-        self.water_mass = builder.get_attribute("water mass")
+        self.signed_water_mass = builder.get_attribute("signed water mass")
 
         super().__init__(
             builder,
             name="relative fall velocity",
-            dependencies=(self.momentum, self.water_mass),
+            dependencies=(self.momentum, self.signed_water_mass),
         )
 
     def recalculate(self):
-        self.data.ratio(self.momentum.get(), self.water_mass.get())
+        self.data.ratio(self.momentum.get(), self.signed_water_mass.get())

PySDM/attributes/physics/terminal_velocity.py

@@ -13,12 +13,28 @@
 class TerminalVelocity(DerivedAttribute):
     def __init__(self, builder):
         self.radius = builder.get_attribute("radius")
-        dependencies = [self.radius]
+        self.signed_water_mass = builder.get_attribute("signed water mass")
+        self.cell_id = builder.get_attribute("cell id")
+        dependencies = [self.radius,self.signed_water_mass,self.cell_id]
         super().__init__(builder, name="terminal velocity", dependencies=dependencies)
 
-        self.approximation = builder.formulae.terminal_velocity_class(
+        self.approximation_liquid = builder.formulae.terminal_velocity_class(
+            builder.particulator
+        )
+        self.approximation_ice = builder.formulae.terminal_velocity_ice_class(
             builder.particulator
         )
 
     def recalculate(self):
-        self.approximation(self.data, self.radius.get())
+        self.approximation_liquid(self.data, self.radius.get())
+        # TODO: fix issue that order of functions calls changes result. approximation_liquid will override
+        #  approximation_ice since r < 0 is not a suitable test for ice particles with mixed-phase spheres shape active
+        if self.formulae.particle_shape_and_density.supports_mixed_phase():
+            self.approximation_ice(self.data,
+                                   self.signed_water_mass.get(),
+                                   self.cell_id.get(),
+                                   self.particulator.environment["T"],
+                                   self.particulator.environment["p"],
+                                   )
+
+

PySDM/backends/impl_numba/methods/terminal_velocity_methods.py

@@ -10,36 +10,41 @@
 
 
 class TerminalVelocityMethods(BackendMethods):
+    # TODO: give terminal velocity functions name of the parametrisation
     @cached_property
     def _interpolation_body(self):
         @numba.njit(**self.default_jit_flags)
         def body(output, radius, factor, b, c):
             for i in numba.prange(len(radius)):  # pylint: disable=not-an-iterable
-                if radius[i] < 0:
-                    output[i] = 0
-                else:
+                if radius[i] > 0:
                     r_id = int(factor * radius[i])
                     r_rest = ((factor * radius[i]) % 1) / factor
                     output[i] = b[r_id] + r_rest * c[r_id]
+                # TODO: check if output 0 for radius 0 is necessary
+                elif radius == 0:
+                    output[i] = 0
 
         return body
 
+    # TODO: give terminal velocity functions name of the parametrisation
     def interpolation(self, *, output, radius, factor, b, c):
         return self._interpolation_body(
             output.data, radius.data, factor, b.data, c.data
         )
-
+    # TODO: give terminal velocity functions name of the parametrisation
     @cached_property
     def _terminal_velocity_body(self):
         v_term = self.formulae.terminal_velocity.v_term
 
         @numba.njit(**self.default_jit_flags)
         def body(*, values, radius):
             for i in numba.prange(len(values)):  # pylint: disable=not-an-iterable
-                values[i] = v_term(radius[i])
+                if radius[i] >= 0.:
+                    values[i] = v_term(radius[i])
 
         return body
 
+    # TODO: give terminal velocity functions name of the parametrisation
     def terminal_velocity(self, *, values, radius):
         self._terminal_velocity_body(values=values, radius=radius)
 
@@ -62,3 +67,54 @@ def power_series(self, *, values, radius, num_terms, prefactors, powers):
             prefactors=prefactors,
             powers=powers,
         )
+
+
+    def terminal_velocity_columnar_ice_crystals(self, *, values, signed_water_mass, cell_id, temperature, pressure):
+        self._terminal_velocity_columnar_ice_crystals_body( values=values,
+                                                            signed_water_mass=signed_water_mass,
+                                                            cell_id=cell_id,
+                                                            temperature=temperature,
+                                                            pressure=pressure,
+                                                            )
+
+    @cached_property
+    def _terminal_velocity_columnar_ice_crystals_body(self):
+        v_base_term = self.formulae.terminal_velocity_ice.v_base_term
+        atmospheric_correction_factor = self.formulae.terminal_velocity_ice.atmospheric_correction_factor
+
+        @numba.njit(**self.default_jit_flags)
+        def body(*, values, signed_water_mass, cell_id, temperature, pressure):
+            for i in numba.prange(len(values)):  # pylint: disable=not-an-iterable
+                if signed_water_mass[i] < 0:
+                    cid = cell_id[i]
+                    correction = atmospheric_correction_factor(temperature[cid], pressure[cid])
+                    values[i] = v_base_term(-signed_water_mass[i]) * correction
+
+        return body
+
+
+
+    def terminal_velocity_ice_spheres(self, *, values, signed_water_mass, cell_id, temperature, pressure):
+        self._terminal_velocity_ice_spheres_body(values=values,
+                                                 signed_water_mass=signed_water_mass,
+                                                 cell_id=cell_id,
+                                                 temperature=temperature,
+                                                 )
+
+    @cached_property
+    def _terminal_velocity_ice_spheres_body(self):
+        v_base_term = self.formulae.terminal_velocity_ice.v_base_term
+        stokes_prefactor = self.formulae.terminal_velocity_ice.stokes_regime
+        formulae = self.formulae_flattened
+
+        def body(*, values, signed_water_mass, cell_id, temperature):
+            for i in numba.prange(len(values)):  # pylint: disable=not-an-iterable
+                if signed_water_mass[i] < 0:
+                    cid = cell_id[i]
+                    radius = formulae.particle_shape_and_density__mass_to_radius(signed_water_mass[i])
+                    dynamic_viscosity = formulae.air_dynamic_viscosity__eta_air(temperature[cid])
+                    prefactor = stokes_prefactor( radius, dynamic_viscosity)
+                    values[i] = v_base_term(prefactor, radius)
+                    print( dynamic_viscosity, signed_water_mass[i], radius, prefactor, values[i] )
+
+        return body

PySDM/dynamics/relaxed_velocity.py

@@ -63,7 +63,7 @@ def register(self, builder):
             "relative fall momentum"
         )
         self.terminal_vel_attr: Attribute = builder.get_attribute("terminal velocity")
-        self.water_mass_attr: Attribute = builder.get_attribute("water mass")
+        self.water_mass_attr: Attribute = builder.get_attribute("signed water mass")
         self.sqrt_radius_attr: Attribute = builder.get_attribute(
             "square root of radius"
         )

PySDM/dynamics/terminal_velocity/__init__.py

@@ -5,3 +5,5 @@
 from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import GunnKinzer1949
 from PySDM.dynamics.terminal_velocity.power_series import PowerSeries
 from PySDM.dynamics.terminal_velocity.rogers_and_yau import RogersYau
+from PySDM.dynamics.terminal_velocity.columnar_ice_crystal import ColumnarIceCrystal
+from PySDM.dynamics.terminal_velocity.spheres_ice import IceSphere

PySDM/dynamics/terminal_velocity/columnar_ice_crystal.py

@@ -0,0 +1,19 @@
+"""
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)
+Eq. (18) .Assumed shape is columnar based on empirical parameterizations of
+[Heymsfield & Iaquinta (2000)](https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2)
+[Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009)
+"""
+
+class ColumnarIceCrystal:  # pylint: disable=too-few-public-methods
+    def __init__(self, particulator):
+        self.particulator = particulator
+
+    def __call__(self, output, signed_water_mass, cell_id, temperature, pressure):
+        self.particulator.backend.terminal_velocity_columnar_ice_crystals(
+            values=output.data,
+            signed_water_mass=signed_water_mass.data,
+            cell_id=cell_id.data,
+            temperature=temperature.data,
+            pressure=pressure.data,
+        )

PySDM/dynamics/terminal_velocity/spheres_ice.py

@@ -0,0 +1,15 @@
+
+
+
+class IceSphere:  # pylint: disable=too-few-public-methods
+    def __init__(self, particulator):
+        self.particulator = particulator
+
+    def __call__(self, output, signed_water_mass, cell_id, temperature, pressure):
+        self.particulator.backend.terminal_velocity_ice_spheres(
+            values=output.data,
+            signed_water_mass=signed_water_mass.data,
+            cell_id=cell_id.data,
+            temperature=temperature.data,
+            pressure=pressure.data,
+        )

PySDM/formulae.py

@@ -19,7 +19,7 @@
 
 from PySDM import physics
 from PySDM.backends.impl_numba import conf
-from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau
+from PySDM.dynamics.terminal_velocity import GunnKinzer1949, PowerSeries, RogersYau, ColumnarIceCrystal, IceSphere
 from PySDM.dynamics.terminal_velocity.gunn_and_kinzer import TpDependent
 
 
@@ -58,6 +58,7 @@ def __init__(  # pylint: disable=too-many-locals
         optical_depth: str = "Null",
         particle_shape_and_density: str = "LiquidSpheres",
         terminal_velocity: str = "GunnKinzer1949",
+        terminal_velocity_ice: str = "ColumnarIceCrystal",
         air_dynamic_viscosity: str = "ZografosEtAl1987",
         bulk_phase_partitioning: str = "Null",
         handle_all_breakups: bool = False,
@@ -96,6 +97,7 @@ def __init__(  # pylint: disable=too-many-locals
         self.particle_shape_and_density = particle_shape_and_density
         self.air_dynamic_viscosity = air_dynamic_viscosity
         self.terminal_velocity = terminal_velocity
+        self.terminal_velocity_ice = terminal_velocity_ice
         self.bulk_phase_partitioning = bulk_phase_partitioning
 
         self._components = tuple(
@@ -158,6 +160,10 @@ def __init__(  # pylint: disable=too-many-locals
             "TpDependent": TpDependent,
             "PowerSeries": PowerSeries,
         }[terminal_velocity]
+        self.terminal_velocity_ice_class = {
+            "ColumnarIceCrystal": ColumnarIceCrystal,
+            "IceSphere": IceSphere,
+        }[terminal_velocity_ice]
 
     def __str__(self):
         description = []

PySDM/physics/__init__.py

@@ -49,6 +49,7 @@
     ventilation,
     air_dynamic_viscosity,
     terminal_velocity,
+    terminal_velocity_ice,
     bulk_phase_partitioning,
 )
 from .constants import convert_to, in_unit, si

PySDM/physics/constants_defaults.py

@@ -580,6 +580,39 @@
 ROGERS_YAU_TERM_VEL_MEDIUM_R_LIMIT = 600 * si.um
 """ 〃 """
 
+SPICHTINGER_GIERENS_TERM_VEL_LIMIT_0 = 2.146e-13 * si.kg
+""" empirical terminal velocity formulation  
+for columnar ice crystals from Table 2. in
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) """
+SPICHTINGER_GIERENS_TERM_VEL_LIMIT_1 = 2.166e-9 * si.kg
+""" 〃 """
+SPICHTINGER_GIERENS_TERM_VEL_LIMIT_2 = 4.264e-8 * si.kg
+""" 〃 """
+SPICHTINGER_TERM_DELTA_COEFF0 = 0.42
+""" 〃 """
+SPICHTINGER_TERM_DELTA_COEFF1 = 0.57
+""" 〃 """
+SPICHTINGER_TERM_DELTA_COEFF2 = 0.31
+""" 〃 """
+SPICHTINGER_TERM_DELTA_COEFF3 = 0.096
+""" 〃 """
+SPICHTINGER_TERM_GAMMA_COEFF0 = 735.4 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF0
+""" 〃 """
+SPICHTINGER_TERM_GAMMA_COEFF1 = 63292.4 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF1
+""" 〃 """
+SPICHTINGER_TERM_GAMMA_COEFF2 = 329.8 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF2
+""" 〃 """
+SPICHTINGER_TERM_GAMMA_COEFF3 = 8.8 * si.m / si.s / si.kg**SPICHTINGER_TERM_DELTA_COEFF3
+""" 〃 """
+SPICHTINGER_CORRECTION_P0 = 300 * si.hectopascal
+""" 〃 """
+SPICHTINGER_CORRECTION_P_EXPO = -0.178
+""" 〃 """
+SPICHTINGER_CORRECTION_T0 = 233 * si.kelvin
+""" 〃 """
+SPICHTINGER_CORRECTION_T_EXPO = -0.394
+""" 〃 """
+
 W76W_G0 = -2.9912729e3 * si.K**2
 """ [Wexler 1976](https://doi.org/10.6028/jres.080A.071) saturation vapour pressure """
 W76W_G1 = -6.0170128e3 * si.K

PySDM/physics/terminal_velocity/__init__.py

@@ -1,4 +1,4 @@
-"""terminal velocity formulae"""
+"""terminal velocity formulae for liquid droplets"""
 
 from .rogers_yau import RogersYau
 from .gunn_kinzer_1949 import GunnKinzer1949

PySDM/physics/terminal_velocity_ice/__init__.py

@@ -0,0 +1,4 @@
+"""terminal velocity formulae for ice particles"""
+
+from .columnar_ice_crystal import ColumnarIceCrystal
+from .spheres_ice import IceSphere

PySDM/physics/terminal_velocity_ice/columnar_ice_crystal.py

@@ -0,0 +1,34 @@
+"""
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)
+Eq. (18) and (19) Assumed shape is columnar based on empirical parameterizations of
+[Heymsfield & Iaquinta (2000)](https://doi.org/10.1175/1520-0469(2000)057<0916:CCTV>2.0.CO;2)
+[Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009)
+"""
+
+import numpy as np
+
+
+class ColumnarIceCrystal:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def v_base_term(const, mass):
+        return np.where(
+            mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_0,
+            const.SPICHTINGER_TERM_GAMMA_COEFF0 * mass**const.SPICHTINGER_TERM_DELTA_COEFF0,
+            np.where(
+                mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_1,
+                const.SPICHTINGER_TERM_GAMMA_COEFF1 * mass**const.SPICHTINGER_TERM_DELTA_COEFF1,
+                np.where(  mass < const.SPICHTINGER_GIERENS_TERM_VEL_LIMIT_2,
+                           const.SPICHTINGER_TERM_GAMMA_COEFF2 * mass**const.SPICHTINGER_TERM_DELTA_COEFF2,
+                           const.SPICHTINGER_TERM_GAMMA_COEFF3 * mass**const.SPICHTINGER_TERM_DELTA_COEFF3
+                )
+            ),
+        )
+
+    @staticmethod
+    def atmospheric_correction_factor(const, temperature, pressure):
+        return( (pressure / const.SPICHTINGER_CORRECTION_P0)**const.SPICHTINGER_CORRECTION_P_EXPO *
+                (const.SPICHTINGER_CORRECTION_T0 / temperature)**const.SPICHTINGER_CORRECTION_T_EXPO
+                )

PySDM/physics/terminal_velocity_ice/spheres_ice.py

@@ -0,0 +1,20 @@
+"""
+terminal velocity of smooth ice spheres
+"""
+
+class IceSphere:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def v_base_term(radius, prefactor):
+        return( prefactor * 2 * radius)
+
+    @staticmethod
+    def stokes_regime(const, radius, dynamic_viscosity):
+        return( const.g_std * const.rho_i * radius / 9 / dynamic_viscosity )
+
+    # TODO: find parametrisation for general functional relationship between reynolds number and drag coefficient
+    @staticmethod
+    def general_flow_regime(const):
+        pass

tests/unit_tests/attributes/test_fall_velocity.py

@@ -77,7 +77,7 @@ def test_fall_velocity_calculation(default_attributes, backend_instance):
         assert np.allclose(
             particulator.attributes["relative fall velocity"].to_ndarray(),
             particulator.attributes["relative fall momentum"].to_ndarray()
-            / (particulator.attributes["water mass"].to_ndarray()),
+            / (particulator.attributes["signed water mass"].to_ndarray()),
         )
 
     @staticmethod
@@ -151,7 +151,7 @@ def test_attribute_selection(backend_instance):
         _ = builder.build(
             attributes={
                 "multiplicity": np.ones(1),
-                "water mass": np.zeros(1),
+                "signed water mass": np.zeros(1),
                 "relative fall momentum": np.zeros(1),
             },
             products=(),
@@ -170,7 +170,7 @@ def test_attribute_selection(backend_instance):
             _ = builder.build(
                 attributes={
                     "multiplicity": np.ones(1),
-                    "water mass": np.zeros(1),
+                    "signed water mass": np.zeros(1),
                 },
                 products=(),
             )