terminal velocity of ice particles (#1562)

Co-authored-by: Sylwester Arabas <[email protected]>

Author: tluettm

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,27 @@
 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 #1605 order of functions calls changes result. approximation_liquid will override
+        #  approximation_ice when 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,38 +10,40 @@
 
 
 class TerminalVelocityMethods(BackendMethods):
+
     @cached_property
-    def _interpolation_body(self):
+    def _gunn_and_kinzer_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]
+                elif radius[i] == 0:
+                    output[i] = 0
 
         return body
 
-    def interpolation(self, *, output, radius, factor, b, c):
-        return self._interpolation_body(
+    def gunn_and_kinzer_interpolation(self, *, output, radius, factor, b, c):
+        return self._gunn_and_kinzer_interpolation_body(
             output.data, radius.data, factor, b.data, c.data
         )
 
     @cached_property
-    def _terminal_velocity_body(self):
+    def _rogers_and_yau_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.0:
+                    values[i] = v_term(radius[i])
 
         return body
 
-    def terminal_velocity(self, *, values, radius):
-        self._terminal_velocity_body(values=values, radius=radius)
+    def rogers_and_yau_terminal_velocity(self, *, values, radius):
+        self._rogers_and_yau_terminal_velocity_body(values=values, radius=radius)
 
     @cached_property
     def _power_series_body(self):
@@ -62,3 +64,64 @@ 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
+    ):  # pylint: disable=unused-argument
+        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)
+
+        return body

PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py

@@ -69,7 +69,7 @@ def __linear_collection_efficiency_body(self):
         )
 
     @cached_property
-    def __interpolation_body(self):
+    def __gunn_and_kinzer_interpolation_body(self):
         # TODO #599 r<0
         return trtc.For(
             ("output", "radius", "factor", "a", "b"),
@@ -126,9 +126,9 @@ def linear_collection_efficiency(
         )
 
     @nice_thrust(**NICE_THRUST_FLAGS)
-    def interpolation(self, *, output, radius, factor, b, c):
+    def gunn_and_kinzer_interpolation(self, *, output, radius, factor, b, c):
         factor_device = trtc.DVInt64(factor)
-        self.__interpolation_body.launch_n(
+        self.__gunn_and_kinzer_interpolation_body.launch_n(
             len(radius), (output.data, radius.data, factor_device, b.data, c.data)
         )
 
@@ -157,7 +157,7 @@ def power_series(self, *, values, radius, num_terms, prefactors, powers):
         )
 
     @cached_property
-    def __terminal_velocity_body(self):
+    def __rogers_and_yau_terminal_velocity_body(self):
         return trtc.For(
             param_names=("values", "radius"),
             name_iter="i",
@@ -171,5 +171,7 @@ def __terminal_velocity_body(self):
         )
 
     @nice_thrust(**NICE_THRUST_FLAGS)
-    def terminal_velocity(self, *, values, radius):
-        self.__terminal_velocity_body.launch_n(n=values.size(), args=[values, radius])
+    def rogers_and_yau_terminal_velocity(self, *, values, radius):
+        self.__rogers_and_yau_terminal_velocity_body.launch_n(
+            n=values.size(), args=[values, radius]
+        )

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,20 @@
+"""
+[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%3C0916:CCTV%3E2.0.CO;2)
+[Barthazy & Schefold (2006)](https://doi.org/10.1016/j.atmosres.2005.12.009)
+"""
+
+
+class ColumnarIceCrystal:  # pylint: disable=too-few-public-methods,too-many-arguments
+    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/gunn_and_kinzer.py

@@ -131,7 +131,7 @@ def __call__(self, output, radius):
                 f"Radii can be interpolated up to {self.maximum_radius} m"
                 + f" (max value of {r_max} m within input data)"
             )
-        self.particulator.backend.interpolation(
+        self.particulator.backend.gunn_and_kinzer_interpolation(
             output=output, radius=radius, factor=self.factor, b=self.a, c=self.b
         )
 

PySDM/dynamics/terminal_velocity/rogers_and_yau.py

@@ -8,7 +8,7 @@ def __init__(self, particulator):
         self.particulator = particulator
 
     def __call__(self, output, radius):
-        self.particulator.backend.terminal_velocity(
+        self.particulator.backend.rogers_and_yau_terminal_velocity(
             values=output.data,
             radius=radius.data,
         )

PySDM/dynamics/terminal_velocity/spheres_ice.py

@@ -0,0 +1,17 @@
+"""
+terminal velocity of solid ice spheres
+"""
+
+
+class IceSphere:  # pylint: disable=too-few-public-methods,too-many-arguments
+    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,
+        )