Merge pull request #731 from slayoo/MAC

mass and heat accommodation coefficients alterable from within constants. closes #730

Author: slayoo

PySDM/backends/impl_numba/methods/condensation_methods.py

@@ -201,7 +201,8 @@ def calculate_ml_old(volume, multiplicity, cell_idx):
     @staticmethod
     def make_calculate_ml_new(
         jit_flags, dx_dt, volume_of_x, x, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
-        phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol, const
+        phys_lambdaK, phys_lambdaD, phys_dk_D, phys_dk_K, within_tolerance, max_iters, RH_rtol,
+        const
     ):
         @numba.njit(**jit_flags)
         def minfun(x_new, x_old, timestep, kappa, f_org, rd3, temperature, RH, lv, pvs, D, K):
@@ -238,8 +239,8 @@ def calculate_ml_new(
                 sgm = phys_sigma(T, v[drop], vdry[drop], f_org[drop])
                 RH_eq = phys_RH_eq(r_old, T, kappa[drop], rd3, sgm)
                 if not within_tolerance(np.abs(RH - RH_eq), RH, RH_rtol):
-                    Dr = phys_DK(DTp, r_old, lambdaD)
-                    Kr = phys_DK(const.K0, r_old, lambdaK)
+                    Dr = phys_dk_D(DTp, r_old, lambdaD)
+                    Kr = phys_dk_K(const.K0, r_old, lambdaK)
                     args = (x_old, timestep, kappa[drop], f_org[drop], rd3, T, RH, lv, pvs, Dr, Kr)
                     r_dr_dt_old = phys_r_dr_dt(RH_eq, T, RH, lv, pvs, Dr, Kr)
                     dx_old = timestep * dx_dt(x_old, r_dr_dt_old)
@@ -326,7 +327,8 @@ def make_condensation_solver(self,
             phys_dthd_dt=self.formulae.state_variable_triplet.dthd_dt,
             phys_lambdaK=self.formulae.diffusion_kinetics.lambdaK,
             phys_lambdaD=self.formulae.diffusion_kinetics.lambdaD,
-            phys_DK=self.formulae.diffusion_kinetics.DK,
+            phys_dk_D=self.formulae.diffusion_kinetics.D,
+            phys_dk_K=self.formulae.diffusion_kinetics.K,
             phys_D=self.formulae.diffusion_thermics.D,
             within_tolerance=self.formulae.trivia.within_tolerance,
             dx_dt=self.formulae.condensation_coordinate.dx_dt,
@@ -346,16 +348,17 @@ def make_condensation_solver(self,
     @lru_cache()
     def make_condensation_solver_impl(
         fastmath, phys_pvs_C, phys_lv, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
-        phys_T, phys_p, phys_pv, phys_dthd_dt, phys_lambdaK, phys_lambdaD, phys_DK, phys_D,
-        within_tolerance, dx_dt, volume, x, timestep, dt_range, adaptive,
+        phys_T, phys_p, phys_pv, phys_dthd_dt, phys_lambdaK, phys_lambdaD, phys_dk_D, phys_dk_K,
+        phys_D, within_tolerance, dx_dt, volume, x, timestep, dt_range, adaptive,
         fuse, multiplier, RH_rtol, max_iters, const
     ):
         jit_flags = {**conf.JIT_FLAGS, **{'parallel': False, 'cache': False, 'fastmath': fastmath}}
 
         calculate_ml_old = CondensationMethods.make_calculate_ml_old(jit_flags, const)
         calculate_ml_new = CondensationMethods.make_calculate_ml_new(
             jit_flags, dx_dt, volume, x, phys_r_dr_dt, phys_RH_eq, phys_sigma, radius,
-            phys_lambdaK, phys_lambdaD, phys_DK, within_tolerance, max_iters, RH_rtol, const
+            phys_lambdaK, phys_lambdaD, phys_dk_D, phys_dk_K, within_tolerance, max_iters, RH_rtol,
+            const
         )
         step_impl = CondensationMethods.make_step_impl(
             jit_flags, phys_pvs_C, phys_lv, calculate_ml_old, calculate_ml_new,

PySDM/backends/impl_numba/test_helpers/bdf.py

@@ -75,7 +75,8 @@ def _make_solve(formulae):
     phys_dthd_dt = formulae.state_variable_triplet.dthd_dt
     phys_lambdaD = formulae.diffusion_kinetics.lambdaD
     phys_lambdaK = formulae.diffusion_kinetics.lambdaK
-    phys_DK = formulae.diffusion_kinetics.DK
+    phys_diff_kin_D = formulae.diffusion_kinetics.D
+    phys_diff_kin_K = formulae.diffusion_kinetics.K
     phys_D = formulae.diffusion_thermics.D
 
     @numba.njit(**{**JIT_FLAGS, **{'parallel': False}})
@@ -91,8 +92,8 @@ def _impl(dy_dt, x, T, p, n, RH, kappa, f_org, dry_volume, thd,
         for i, x_i in enumerate(x):
             v = volume(x_i)
             r = phys_radius(v)
-            Dr = phys_DK(DTp, r, lambdaD)
-            Kr = phys_DK(K0, r, lambdaK)
+            Dr = phys_diff_kin_D(DTp, r, lambdaD)
+            Kr = phys_diff_kin_K(K0, r, lambdaK)
             sgm = sigma(T, v, dry_volume[i], f_org[i])
             dy_dt[idx_x + i] = dx_dt(
                 x_i,

PySDM/backends/impl_thrust_rtc/methods/condensation_methods.py

@@ -109,11 +109,11 @@ def args(arg):
             if ( ! {phys.trivia.within_tolerance.c_inline(
                 return_type='bool', error_estimate="abs(_RH - RH_eq)", value="_RH", rtol="RH_rtol"
             )}) {{
-                Dr = {phys.diffusion_kinetics.DK.c_inline(
-                    DK="_DTp", r="r_old", lmbd="_lambdaD"
+                Dr = {phys.diffusion_kinetics.D.c_inline(
+                    D="_DTp", r="r_old", lmbd="_lambdaD"
                 )};
-                Kr = {phys.diffusion_kinetics.DK.c_inline(
-                    DK="const.K0", r="r_old", lmbd="_lambdaK"
+                Kr = {phys.diffusion_kinetics.K.c_inline(
+                    K="const.K0", r="r_old", lmbd="_lambdaK"
                 )};
                 r_dr_dt_old = {phys.drop_growth.r_dr_dt.c_inline(
                     RH_eq="RH_eq", T="_T", RH="_RH", lv="_lv", pvs="_pvs", D="Dr", K="Kr"

PySDM/physics/constants_defaults.py

@@ -24,6 +24,10 @@
 
 K0 = 2.4e-2 * si.joules / si.metres / si.seconds / si.kelvins
 
+# mass and heat accommodation coefficients
+MAC = 1.
+HAC = 1.
+
 p1000 = 1000 * si.hectopascals
 c_pd = 1005 * si.joule / si.kilogram / si.kelvin
 c_pv = 1850 * si.joule / si.kilogram / si.kelvin

PySDM/physics/diffusion_kinetics/fuchs_sutugin.py

@@ -18,5 +18,13 @@ def lambdaK(const, T, p):
         return (4. / 5) * const.K0 * T / p / np.sqrt(2 * const.Rd * T)
 
     @staticmethod
-    def DK(_, DK, r, lmbd):
-        return DK * (1 + lmbd/r) / (1 + 1.71 * lmbd/r + 1.33 * lmbd/r * lmbd/r)
+    def D(const, D, r, lmbd):
+        return D * (1 + lmbd/r) / (
+                1 + (4./3/const.MAC + 0.377) * lmbd/r + (4./3/const.MAC) * lmbd/r * lmbd/r
+        )
+
+    @staticmethod
+    def K(const, K, r, lmbd):
+        return K * (1 + lmbd/r) / (
+                1 + (4./3/const.HAC + 0.377) * lmbd/r + (4./3/const.HAC) * lmbd/r * lmbd/r
+        )

PySDM/physics/diffusion_kinetics/neglect.py

@@ -16,5 +16,9 @@ def lambdaK(_, T, p):  # pylint: disable=unused-argument
         return -1
 
     @staticmethod
-    def DK(_, DK, r, lmbd):  # pylint: disable=unused-argument
-        return DK
+    def D(_, D, r, lmbd):  # pylint: disable=unused-argument
+        return D
+
+    @staticmethod
+    def K(_, K, r, lmbd):  # pylint: disable=unused-argument
+        return K

tests/unit_tests/physics/test_accommodation_coefficients.py

@@ -0,0 +1,29 @@
+# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring
+import pytest
+import numpy as np
+from PySDM import Formulae
+
+@pytest.mark.parametrize("constants", (
+    {"MAC": 1., "HAC": 1.},
+    {"MAC": 1., "HAC": .1},
+    {"MAC": .1, "HAC": 1.},
+    {"MAC": .1, "HAC": .1},
+))
+def test_accomodation_coefficients(constants):
+    # arrange
+    formulae = Formulae(constants=constants)
+    D = 1
+    K = 2
+    r = 3
+    lmbd = 4
+
+    # act
+    D_dk = formulae.diffusion_kinetics.D(D, r, lmbd)
+    K_dk = formulae.diffusion_kinetics.K(K, r, lmbd)
+
+    # assert
+    Kn = lmbd / r
+    xx_D = 4 / 3 / constants['MAC']
+    np.testing.assert_almost_equal(D_dk, D * (1 + Kn) / (1 + (xx_D + .377) * Kn + xx_D * Kn**2))
+    xx_K = 4 / 3 / constants['HAC']
+    np.testing.assert_almost_equal(K_dk, K * (1 + Kn) / (1 + (xx_K + .377) * Kn + xx_K * Kn**2))