ice diffusion capacity: refactor for spherical & added columnar (+ several speedups in unrelated codes which started timeouting due to longer compilation time) (#1643)

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

Author: tluettm

.github/workflows/tests.yml

@@ -170,6 +170,7 @@ jobs:
         run: |
           python -m pip uninstall -y jupyterlab-server # https://github.com/pypa/pip/issues/6275
           python -m pip install -e . -e ./examples
+          python -m pip install pytest-timeout
           python -m pip install -r tests/devops_tests/requirements.txt
       - if: steps.cache.outputs.cache-hit != 'true'
         uses: actions/cache/save@v4
@@ -245,7 +246,7 @@ jobs:
         run: python -m pytest --basetemp=/tmp/pytest ${{ env.pytest_options }} tests/examples_tests/test_run* -k "not Rozanski_and_Sonntag_1982" --suite ${{ matrix.test-suite }}
 
       - if: ( ! startsWith(matrix.platform, 'macos-13') )
-        run: python -m pytest ${{ env.pytest_options }} --basetemp=/tmp/pytest tests/examples_tests/test_run* --suite ${{ matrix.test-suite }}
+        run: python -m pytest ${{ env.pytest_options }} --basetemp=/tmp/pytest --timeout_method=thread --timeout=${{ startsWith(matrix.platform, 'windows-') && 1000 || 900 }} tests/examples_tests/test_run* --suite ${{ matrix.test-suite }}
 
       - if: ( ! startsWith(matrix.platform, 'windows-') ) && matrix.test-suite == env.anim_test-suite && matrix.python-version == env.anim_python-version
         run: |

PySDM/backends/impl_numba/methods/deposition_methods.py

@@ -16,7 +16,7 @@ def _deposition(self):
         formulae = self.formulae_flattened
 
         @numba.njit(**{**self.default_jit_flags, **{"parallel": False}})
-        def body(  # pylint: disable=too-many-arguments
+        def body(  # pylint: disable=too-many-arguments,unused-argument
             *,
             multiplicity,
             signed_water_mass,
@@ -47,22 +47,15 @@ def body(  # pylint: disable=too-many-arguments
                         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)
+                    capacity = formulae.diffusion_ice_capacity__capacity(ice_mass)
 
-                    mass_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],
-                        )
-                    )
+                    mass_ventilation_factor = 1  # TODO #1655
                     heat_ventilation_factor = mass_ventilation_factor  # TODO #1588
 
                     Dv_const = formulae.diffusion_thermics__D(temperature, pressure)
@@ -95,6 +88,7 @@ def body(  # pylint: disable=too-many-arguments
                         D=diffusion_coefficient * mass_ventilation_factor,
                         pvs=pvs_ice,
                     )
+
                     howell_factor_x_diffcoef_x_rhovsice_x_icess = (
                         formulae.drop_growth__r_dr_dt(
                             RH_eq=1,

PySDM/physics/constants_defaults.py

@@ -90,22 +90,24 @@
 N_A = sci.N_A / si.mole
 """ Avogadro constant (value from SciPy) """
 
-# mass and heat accommodation coefficients for condensation
 MAC = 1.0
 """ mass accommodation coefficient of unity as recommended in
 [Laaksonen et al. 2005](https://doi.org/10.5194/acp-5-461-2005) """
 HAC = 1.0
 """ thermal accommodation coefficient of unity as recommended in
 [Laaksonen et al. 2005](https://doi.org/10.5194/acp-5-461-2005) """
 
-# mass and heat accommodation coefficients for vapour deposition on ice
 MAC_ice = 0.5
 """ mass accommodation coefficient for vapour deposition as recommended in
 [Kaercher & Lohmann 2002](https://doi.org/10.1029/2001JD000470) """
-HAC_ice = 0.7
+HAC_ice = 1.0
 """ thermal accommodation coefficient for vapour deposition as recommended in
 [Pruppacher & Klett](https://doi.org/10.1007/978-0-306-48100-0) """
 
+C_cunn = 0.7
+""" Cunningham correction factor as used in
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) """
+
 ARM_C1 = 6.1094 * si.hectopascal
 """ [August](https://doi.org/10.1002/andp.18280890511) Roche Magnus formula coefficients
 (values from [Alduchov & Eskridge 1996](https://doi.org/10.1175%2F1520-0450%281996%29035%3C0601%3AIMFAOS%3E2.0.CO%3B2))
@@ -527,6 +529,28 @@
 CRAIG_1961_INTERCEPT_COEFF = 10 * PER_MILLE
 """ 〃 """
 
+capacity_columnar_ice_B1 = 0.3
+""" eq. A11 & A12 in [Spichtinger et al. 2023](https://doi.org/10.5194/acp-23-2035-2023)  """
+capacity_columnar_ice_B2 = 0.43
+""" 〃 """
+capacity_columnar_ice_A1 = 0.015755 * si.m / si.kg**capacity_columnar_ice_B1
+""" 〃 """
+capacity_columnar_ice_A2 = 0.33565 * si.m / si.kg**capacity_columnar_ice_B2
+""" 〃 """
+
+columnar_ice_mass_transition = 2.146e-13 * si.kg
+""" tab. 1 in [Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009) """
+columnar_ice_length_beta_1 = 3.0
+""" 〃 """
+columnar_ice_length_beta_2 = 2.2
+""" 〃 """
+columnar_ice_length_alpha_1 = 526.1 * si.kg / si.m**columnar_ice_length_beta_1
+""" 〃 """
+columnar_ice_length_alpha_2 = 0.04142 * si.kg / si.m**columnar_ice_length_beta_2
+""" 〃 """
+columnar_bulk_ice_density = 0.81e3 * si.kg / si.m**3
+""" 〃 """
+
 asymmetry_g = 0.85  # forward scattering from cloud droplets
 """ [Bohren 1987](https://doi.org/10.1119/1.15109) """
 

PySDM/physics/diffusion_ice_capacity/__init__.py

@@ -3,3 +3,4 @@
 """
 
 from .spherical import Spherical
+from .columnar import Columnar

PySDM/physics/diffusion_ice_capacity/columnar.py

@@ -0,0 +1,29 @@
+"""
+capacity for columnar ice crystals approximated as prolate ellipsoids
+eq. A11 & A12 in Spichtinger et al. 2023 (https://doi.org/10.5194/acp-23-2035-2023)
+"""
+
+import numpy as np
+
+
+class Columnar:  # pylint: disable=too-few-public-methods
+
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def capacity(const, mass):
+        return (
+            const.capacity_columnar_ice_A1 * mass**const.capacity_columnar_ice_B1
+            + const.capacity_columnar_ice_A2 * mass**const.capacity_columnar_ice_B2
+        )
+
+    @staticmethod
+    def reference_capacity(
+        const, polar_diameter, eccentricity
+    ):  # pylint: disable=unused-argument
+        return (
+            polar_diameter
+            * eccentricity
+            / np.log((1 + eccentricity) / (1 - eccentricity))
+        )

PySDM/physics/diffusion_ice_capacity/spherical.py

@@ -2,12 +2,14 @@
 capacity for approximation of ice crystals as spheres
 """
 
+import numpy as np
+
 
 class Spherical:  # pylint: disable=too-few-public-methods
 
     def __init__(self, _):
         pass
 
     @staticmethod
-    def capacity(const, diameter):
-        return diameter / const.TWO
+    def capacity(const, mass):
+        return np.power(mass / const.PI_4_3 / const.rho_i, const.ONE_THIRD)

PySDM/physics/diffusion_ice_kinetics/standard.py

@@ -23,16 +23,15 @@ def lambdaK(const, T, p):
     @staticmethod
     def D(const, D, r, lmbd, T):
         return D / (
-            r / (r + lmbd)
+            r / (r + lmbd * const.C_cunn)
             + 4.0 * D / const.MAC_ice / np.sqrt(8.0 * const.Rv * T / const.PI) / r
         )
 
     @staticmethod
     def K(const, K, r, lmbd, T, rho):  # pylint: disable=too-many-arguments
         return K / (
             r / (r + lmbd)
-            + 4.0
-            * K
+            + K
             / const.HAC_ice
             / np.sqrt(8.0 * const.Rd * T / const.PI)
             / const.c_pd

PySDM/physics/particle_shape_and_density/__init__.py

@@ -5,3 +5,4 @@
 from .liquid_spheres import LiquidSpheres
 from .mixed_phase_spheres import MixedPhaseSpheres
 from .porous_spheroids import PorousSpheroid
+from .columnar_ice import ColumnarIce

PySDM/physics/particle_shape_and_density/columnar_ice.py

@@ -0,0 +1,46 @@
+"""
+for columnar ice crystals as in
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)
+"""
+
+import numpy as np
+from .porous_spheroids import PorousSpheroid
+
+
+class ColumnarIce(PorousSpheroid):
+
+    @staticmethod
+    def polar_radius_empirical_parametrisation(const, mass):
+        """Sec. 3.1.2 in [Spichtinger & Gierens 2009].
+        Note that the notation in the paper is incorrect.
+        """
+        return np.where(
+            mass < const.columnar_ice_mass_transition,
+            (mass / const.columnar_ice_length_alpha_1)
+            ** (1 / const.columnar_ice_length_beta_1)
+            / 2,
+            (mass / const.columnar_ice_length_alpha_2)
+            ** (1 / const.columnar_ice_length_beta_2)
+            / 2,
+        )
+
+    @staticmethod
+    def aspect_ratio_empirical_parametrisation(const, mass):
+        """Eq. 17 in [Spichtinger & Gierens 2009]"""
+        return np.where(
+            mass < const.columnar_ice_mass_transition,
+            1,
+            np.sqrt(
+                np.sqrt(27)
+                * const.columnar_bulk_ice_density
+                / 8.0
+                / const.columnar_ice_length_alpha_2
+                ** (3 / const.columnar_ice_length_beta_2)
+            )
+            * mass
+            ** (
+                (3 - const.columnar_ice_length_beta_2)
+                / 2
+                / const.columnar_ice_length_beta_2
+            ),
+        )

PySDM/physics/particle_shape_and_density/porous_spheroids.py

@@ -1,10 +1,33 @@
 """
-for mixed-phase microphysics as in
+for porous spheriods  as in
 [Shima et al. 2020](https://doi.org/10.5194/gmd-13-4107-2020)
+and prolate spheriods as in
+[Spichtinger & Gierens 2009](https://doi.org/10.5194/acp-9-685-2009)
 """
 
+import numpy as np
+
 
 class PorousSpheroid:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
     @staticmethod
     def supports_mixed_phase(_=None):
         return True
+
+    @staticmethod
+    def aspect_ratio(polar_radius, equatorial_radius):
+        return polar_radius / equatorial_radius
+
+    @staticmethod
+    def equatorial_radius(polar_radius, aspect_ratio):
+        return polar_radius / aspect_ratio
+
+    @staticmethod
+    def polar_radius(equatorial_radius, aspect_ratio):
+        return equatorial_radius * aspect_ratio
+
+    @staticmethod
+    def eccentricity(aspect_ratio):
+        return np.sqrt(1 - aspect_ratio**-2.0)