new examples & smoke tests: Gedzelman & Arnold 1994 fig 2, Miyake et al. 1968 fig 19; isotope diffusivity ratios in Formulae (Graham's law; Stewart 1975; Hellmann & Harvey); isotope relaxation timescales formulae (Miyake); new ventilation test based on Kinzer & Gunn 1951 table 1; null refactor in Formulae; access to RogersYau terminal vel. from Formulae

.github/workflows/tests+artifacts+pypi.yml

@@ -132,7 +132,7 @@ jobs:
             python-version: "3.8"
       fail-fast: false
     runs-on: ${{ matrix.platform }}
-    timeout-minutes: 40
+    timeout-minutes: 45
     steps:
       - uses: actions/[email protected]
         with:

PySDM/backends/impl_numba/methods/terminal_velocity_methods.py

@@ -30,22 +30,17 @@ def interpolation(self, *, output, radius, factor, b, c):
 
     @cached_property
     def terminal_velocity_body(self):
+        v_term = self.formulae.terminal_velocity.v_term
+
         @numba.njit(**{**conf.JIT_FLAGS, "fastmath": self.formulae.fastmath})
-        def terminal_velocity_body(*, values, radius, k1, k2, k3, r1, r2):
+        def terminal_velocity_body(*, values, radius):
             for i in numba.prange(len(values)):  # pylint: disable=not-an-iterable
-                if radius[i] < r1:
-                    values[i] = k1 * radius[i] ** 2
-                elif radius[i] < r2:
-                    values[i] = k2 * radius[i]
-                else:
-                    values[i] = k3 * radius[i] ** (1 / 2)
+                values[i] = v_term(radius[i])
 
         return terminal_velocity_body
 
-    def terminal_velocity(self, *, values, radius, k1, k2, k3, r1, r2):
-        self.terminal_velocity_body(
-            values=values, radius=radius, k1=k1, k2=k2, k3=k3, r1=r1, r2=r2
-        )
+    def terminal_velocity(self, *, values, radius):
+        self.terminal_velocity_body(values=values, radius=radius)
 
     @cached_property
     def power_series_body(self):

PySDM/backends/impl_thrust_rtc/methods/physics_methods.py

@@ -68,28 +68,6 @@ def __critical_volume_body(self):
             ),
         )
 
-    @cached_property
-    def __terminal_velocity_body(self):
-        return trtc.For(
-            ("values", "radius", "k1", "k2", "k3", "r1", "r2"),
-            "i",
-            """
-            if (radius[i] < r1) {
-                values[i] = k1 * radius[i] * radius[i];
-            }
-            else {
-                if (radius[i] < r2) {
-                    values[i] = k2 * radius[i];
-                }
-                else {
-                    values[i] = k3 * pow(radius[i], (real_type)(.5));
-                }
-            }
-            """.replace(
-                "real_type", self._get_c_type()
-            ),
-        )
-
     @cached_property
     def __volume_of_mass_body(self):
         return trtc.For(
@@ -145,17 +123,6 @@ def temperature_pressure_RH(
             ),
         )
 
-    @nice_thrust(**NICE_THRUST_FLAGS)
-    def terminal_velocity(self, *, values, radius, k1, k2, k3, r1, r2):
-        k1 = self._get_floating_point(k1)
-        k2 = self._get_floating_point(k2)
-        k3 = self._get_floating_point(k3)
-        r1 = self._get_floating_point(r1)
-        r2 = self._get_floating_point(r2)
-        self.__terminal_velocity_body.launch_n(
-            values.size(), [values, radius, k1, k2, k3, r1, r2]
-        )
-
     @nice_thrust(**NICE_THRUST_FLAGS)
     def explicit_euler(self, y, dt, dy_dt):
         dt = self._get_floating_point(dt)

PySDM/backends/impl_thrust_rtc/methods/terminal_velocity_methods.py

@@ -155,3 +155,21 @@ def power_series(self, *, values, radius, num_terms, prefactors, powers):
         self.__power_series_body.launch_n(
             values.size(), (values, radius, num_terms, prefactors, powers)
         )
+
+    @cached_property
+    def __terminal_velocity_body(self):
+        return trtc.For(
+            param_names=("values", "radius"),
+            name_iter="i",
+            body=f"""
+            values[i] = {self.formulae.terminal_velocity.v_term.c_inline(
+                radius="radius[i]"
+            )};
+            """.replace(
+                "real_type", self._get_c_type()
+            ),
+        )
+
+    @nice_thrust(**NICE_THRUST_FLAGS)
+    def terminal_velocity(self, *, values, radius):
+        self.__terminal_velocity_body.launch_n(n=values.size(), args=[values, radius])

PySDM/dynamics/terminal_velocity/rogers_and_yau.py

@@ -2,35 +2,13 @@
 Rogers & Yau, equations: (8.5), (8.6), (8.8)
 """
 
-from PySDM.physics import constants as const
-
 
 class RogersYau:  # pylint: disable=too-few-public-methods
-    def __init__(
-        self,
-        particulator,
-        *,
-        small_k=None,
-        medium_k=None,
-        large_k=None,
-        small_r_limit=None,
-        medium_r_limit=None,
-    ):
-        si = const.si
+    def __init__(self, particulator):
         self.particulator = particulator
-        self.small_k = small_k or 1.19e6 / si.cm / si.s
-        self.medium_k = medium_k or 8e3 / si.s
-        self.large_k = large_k or 2.01e3 * si.cm ** (1 / 2) / si.s
-        self.small_r_limit = small_r_limit or 35 * si.um
-        self.medium_r_limit = medium_r_limit or 600 * si.um
 
     def __call__(self, output, radius):
         self.particulator.backend.terminal_velocity(
             values=output.data,
             radius=radius.data,
-            k1=self.small_k,
-            k2=self.medium_k,
-            k3=self.large_k,
-            r1=self.small_r_limit,
-            r2=self.medium_r_limit,
         )

PySDM/formulae.py

@@ -48,6 +48,8 @@ def __init__(  # pylint: disable=too-many-locals
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_meteoric_water_line_excess: str = "Null",
         isotope_ratio_evolution: str = "Null",
+        isotope_diffusivity_ratios: str = "Null",
+        isotope_relaxation_timescale: str = "Null",
         optical_albedo: str = "Null",
         optical_depth: str = "Null",
         particle_shape_and_density: str = "LiquidSpheres",
@@ -79,8 +81,11 @@ def __init__(  # pylint: disable=too-many-locals
         )
         self.isotope_meteoric_water_line_excess = isotope_meteoric_water_line_excess
         self.isotope_ratio_evolution = isotope_ratio_evolution
+        self.isotope_diffusivity_ratios = isotope_diffusivity_ratios
+        self.isotope_relaxation_timescale = isotope_relaxation_timescale
         self.particle_shape_and_density = particle_shape_and_density
         self.air_dynamic_viscosity = air_dynamic_viscosity
+        self.terminal_velocity = terminal_velocity
 
         self._components = tuple(
             i
@@ -298,6 +303,8 @@ def _pick(value: str, choices: dict, constants: namedtuple):
         for name, cls in choices.items():
             if name == value:
                 obj = cls(constants)
+                break
+            name = value
     else:
         parent_classes = []
         for name in value.split("+"):
@@ -317,7 +324,7 @@ def __init__(self, const):
 
     if obj is None:
         raise ValueError(
-            f"Unknown setting: '{name}'; choices are: {tuple(choices.keys())}"
+            f"Unknown setting: {name}; choices are: {', '.join(choices.keys())}"
         )
 
     obj.__name__ = value  # pylint: disable=attribute-defined-outside-init

PySDM/physics/__init__.py

@@ -31,6 +31,8 @@
     isotope_equilibrium_fractionation_factors,
     isotope_meteoric_water_line_excess,
     isotope_ratio_evolution,
+    isotope_diffusivity_ratios,
+    isotope_relaxation_timescale,
     latent_heat,
     optical_albedo,
     optical_depth,
@@ -42,5 +44,6 @@
     trivia,
     ventilation,
     air_dynamic_viscosity,
+    terminal_velocity,
 )
 from .constants import convert_to, in_unit, si

PySDM/physics/constants.py

@@ -43,6 +43,7 @@ def in_unit(value, unit):
 ONE_HALF = 1 / 2
 TWO_THIRDS = 2 / 3
 ONE_AND_A_HALF = 3 / 2
+TWO_AND_A_HALF = 5 / 2
 NaN = np.nan
 
 default_random_seed = (

PySDM/physics/constants_defaults.py

@@ -310,15 +310,15 @@
 """ [Bohren 1987](https://doi.org/10.1119/1.15109) """
 asymmetry_g = 0.85  # forward scattering from cloud droplets
 
-""" TODO #1266 """
-diffussion_thermics_D_G11_A = 1e-5 * si.m**2 / si.s
-diffussion_thermics_D_G11_B = 0.015 / si.K
-diffussion_thermics_D_G11_C = -1.9
+""" [Grabowski et al. 2011](https://doi.org/10.1016/j.atmosres.2010.10.020) """
+diffusion_thermics_D_G11_A = 1e-5 * si.m**2 / si.s
+diffusion_thermics_D_G11_B = 0.015 / si.K
+diffusion_thermics_D_G11_C = -1.9
 
-diffussion_thermics_K_G11_A = 1.5e-11 * si.W / si.m / si.K**4
-diffussion_thermics_K_G11_B = -4.8e-8 * si.W / si.m / si.K**3
-diffussion_thermics_K_G11_C = 1e-4 * si.W / si.m / si.K**2
-diffussion_thermics_K_G11_D = -3.9e-4 * si.W / si.m / si.K
+diffusion_thermics_K_G11_A = 1.5e-11 * si.W / si.m / si.K**4
+diffusion_thermics_K_G11_B = -4.8e-8 * si.W / si.m / si.K**3
+diffusion_thermics_K_G11_C = 1e-4 * si.W / si.m / si.K**2
+diffusion_thermics_K_G11_D = -3.9e-4 * si.W / si.m / si.K
 
 """
 [Pruppacher & Rasmussen (1979))](https://doi.org/10.1175/1520-0469(1979)036<1255:AWTIOT>2.0.CO;2)
@@ -343,11 +343,32 @@
 FROESSLING_1938_B = 0.276
 
 
+""" fit coefficients from [Hellmann & Harvey 2020](https://doi.org/10.1029/2020GL089999) """
+HELLMANN_HARVEY_T_UNIT = 100 * si.K
+HELLMANN_HARVEY_EQ6_COEFF0 = 0.98258
+HELLMANN_HARVEY_EQ6_COEFF1 = -0.02546
+HELLMANN_HARVEY_EQ6_COEFF2 = 0.02421
+HELLMANN_HARVEY_EQ7_COEFF0 = 0.98284
+HELLMANN_HARVEY_EQ7_COEFF1 = 0.003517
+HELLMANN_HARVEY_EQ7_COEFF2 = -0.001996
+HELLMANN_HARVEY_EQ8_COEFF0 = 0.96671
+HELLMANN_HARVEY_EQ8_COEFF1 = 0.007406
+HELLMANN_HARVEY_EQ8_COEFF2 = -0.004861
+
+
+""" terminal velocity formulation from Rogers & Yau (equations: 8.5, 8.6, 8.8) """
+ROGERS_YAU_TERM_VEL_SMALL_K = 1.19e6 / si.cm / si.s
+ROGERS_YAU_TERM_VEL_MEDIUM_K = 8e3 / si.s
+ROGERS_YAU_TERM_VEL_LARGE_K = 2.01e3 * si.cm**0.5 / si.s
+ROGERS_YAU_TERM_VEL_SMALL_R_LIMIT = 35 * si.um
+ROGERS_YAU_TERM_VEL_MEDIUM_R_LIMIT = 600 * si.um
+
+
 def compute_derived_values(c: dict):
     """
     computes derived quantities such as molar mass ratios, etc.
 
-    water molar mass is computed from from molecular masses and VSMOW isotope abundances
+    water molar mass is computed from molecular masses and VSMOW isotope abundances
     (and neglecting molecular binding energies)
     for discussion, see:
     - caption of Table 2.1 in [Gat 2010](https://doi.org/10.1142/p027)

PySDM/physics/diffusion_thermics/grabowski_et_al_2011.py

@@ -10,16 +10,16 @@ def __init__(self, _):
     @staticmethod
     def D(const, T, p):  # pylint: disable=unused-argument
         """eq (10)"""
-        return const.diffussion_thermics_D_G11_A * (
-            const.diffussion_thermics_D_G11_B * T + const.diffussion_thermics_D_G11_C
+        return const.diffusion_thermics_D_G11_A * (
+            const.diffusion_thermics_D_G11_B * T + const.diffusion_thermics_D_G11_C
         )
 
     @staticmethod
     def K(const, T, p):  # pylint: disable=unused-argument
         """eq (12)"""
         return (
-            const.diffussion_thermics_K_G11_A * T**3
-            + const.diffussion_thermics_K_G11_B * T**2
-            + const.diffussion_thermics_K_G11_C * T
-            + const.diffussion_thermics_K_G11_D
+            const.diffusion_thermics_K_G11_A * T**3
+            + const.diffusion_thermics_K_G11_B * T**2
+            + const.diffusion_thermics_K_G11_C * T
+            + const.diffusion_thermics_K_G11_D
         )

PySDM/physics/isotope_diffusivity_ratios/__init__.py

@@ -0,0 +1,6 @@
+""" heavy-to-light diffusivity ratios for water molecules """
+
+from PySDM.impl.null_physics_class import Null
+from .grahams_law import GrahamsLaw
+from .stewart_1975 import Stewart1975
+from .hellmann_and_harvey_2020 import HellmannAndHarvey2020

PySDM/physics/isotope_diffusivity_ratios/grahams_law.py

@@ -0,0 +1,14 @@
+""" [Graham's law](https://en.wikipedia.org/wiki/Graham%27s_law)
+    see also eq. (21) in [Horita et al. 2008](https://doi.org/10.1080/10256010801887174)
+"""
+
+
+class GrahamsLaw:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def ratio_2H(const, temperature):  # pylint: disable=unused-argument
+        return (
+            (2 * const.M_1H + const.M_16O) / (const.M_2H + const.M_1H + const.M_16O)
+        ) ** const.ONE_HALF

PySDM/physics/isotope_diffusivity_ratios/hellmann_and_harvey_2020.py

@@ -0,0 +1,39 @@
+"""
+Ratios of diffusivity in air of heavy vs. light isotope using fits provided in
+[Hellmann & Harvey 2020](https://doi.org/10.1029/2020GL089999)
+"""
+
+
+class HellmannAndHarvey2020:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def ratio_2H(const, temperature):
+        return (
+            const.HELLMANN_HARVEY_EQ6_COEFF0
+            + const.HELLMANN_HARVEY_EQ6_COEFF1
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT)
+            + const.HELLMANN_HARVEY_EQ6_COEFF2
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT) ** (const.TWO_AND_A_HALF)
+        )
+
+    @staticmethod
+    def ratio_17O(const, temperature):
+        return (
+            const.HELLMANN_HARVEY_EQ7_COEFF0
+            + const.HELLMANN_HARVEY_EQ7_COEFF1
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT) ** (const.ONE_HALF)
+            + const.HELLMANN_HARVEY_EQ7_COEFF2
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT) ** (const.TWO_AND_A_HALF)
+        )
+
+    @staticmethod
+    def ratio_18O(const, temperature):
+        return (
+            const.HELLMANN_HARVEY_EQ8_COEFF0
+            + const.HELLMANN_HARVEY_EQ8_COEFF1
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT) ** (const.ONE_HALF)
+            + const.HELLMANN_HARVEY_EQ8_COEFF2
+            / (temperature / const.HELLMANN_HARVEY_T_UNIT) ** (const.THREE)
+        )

PySDM/physics/isotope_diffusivity_ratios/stewart_1975.py

@@ -0,0 +1,40 @@
+"""
+Temperature-independent ratio of vapour diffusivity in air for heavy vs. light isotopes
+assuming same collision diameters for different isotopic water molecules, see:
+eq. (8) in [Stewart 1975](https://doi.org/10.1029/JC080i009p01133),
+eq. (1) in [Merlivat 1978](https://doi.org/10.1063/1.436884),
+eq. (6) in [Cappa et al. 2003](https://doi.org/10.1029/2003JD003597),
+eq. (22) in [Horita et al. 2008](https://doi.org/10.1080/10256010801887174),
+eq. (3) in [Hellmann and Harvey 2020](https://doi.org/10.1029/2020GL089999).
+
+All functions return constants, so there is a potential overhead in computing them on each call,
+but this variant is provided for historical reference only, hence leaving like that.
+"""
+
+
+class Stewart1975:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def ratio_2H(const, temperature):  # pylint: disable=unused-argument
+        return (
+            (
+                (2 * const.M_1H + const.M_16O)
+                * (const.Md + const.M_2H + const.M_1H + const.M_16O)
+            )
+            / (
+                (const.M_2H + const.M_1H + const.M_16O)
+                * (const.Md + (2 * const.M_1H + const.M_16O))
+            )
+        ) ** const.ONE_HALF
+
+    @staticmethod
+    def ratio_18O(const, temperature):  # pylint: disable=unused-argument
+        return (
+            ((2 * const.M_1H + const.M_16O) * (const.Md + 2 * const.M_1H + const.M_18O))
+            / (
+                (2 * const.M_1H + const.M_18O)
+                * (const.Md + (2 * const.M_1H + const.M_16O))
+            )
+        ) ** const.ONE_HALF