add Mason's approximated drop growth formula in the very form given i…

…n the 1951 paper (without the -1 term) + unit test (#1320)

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

@@ -58,7 +58,7 @@ jobs:
   pdoc:
     strategy:
       matrix:
-        platform: [ubuntu-latest, macos-latest, windows-latest]
+        platform: [ubuntu-latest, macos-12, windows-latest]
     runs-on: ${{ matrix.platform }}
     steps:
       - uses: actions/checkout@v4.1.1

PySDM/formulae.py

@@ -34,7 +34,7 @@ def __init__(  # pylint: disable=too-many-locals
         saturation_vapour_pressure: str = "FlatauWalkoCotton",
         latent_heat: str = "Kirchhoff",
         hygroscopicity: str = "KappaKoehlerLeadingTerms",
-        drop_growth: str = "MaxwellMason",
+        drop_growth: str = "Mason1971",
         surface_tension: str = "Constant",
         diffusion_kinetics: str = "FuchsSutugin",
         diffusion_thermics: str = "Neglect",

PySDM/physics/drop_growth/__init__.py

@@ -2,4 +2,5 @@
 Formulation of the coupled heat-moisture diffusion problem
 """
 
-from .maxwell_mason import MaxwellMason
+from .mason_1951 import Mason1951
+from .mason_1971 import Mason1971

PySDM/physics/drop_growth/mason_1951.py

@@ -0,0 +1,23 @@
+"""
+single-equation approximation of the vapour and heat diffusion problem
+as derived in [Mason 1951](https://doi.org/10.1088/0370-1301/64/9/307)
+(note that Mason's work did not include the ventilation factor, it
+is included here for code maintainability - to reduce duplication at the
+calling scopes; it is not ignored and is used as a multiplicative factor
+in the same way as in `PySDM.physics.drop_growth.mason_1971.Mason1971`)
+"""
+
+
+class Mason1951:  # pylint: disable=too-few-public-methods
+    def __init__(self, _):
+        pass
+
+    # pylint: disable=too-many-arguments
+    @staticmethod
+    def r_dr_dt(const, RH_eq, T, RH, lv, pvs, D, K, ventilation_factor):
+        return (
+            ventilation_factor
+            * (RH - RH_eq)
+            / const.rho_w
+            / (const.Rv * T / D / pvs + lv**2 / K / T**2 / const.Rv)
+        )

PySDM/physics/drop_growth/maxwell_mason.py → PySDM/physics/drop_growth/mason_1971.py

@@ -1,9 +1,11 @@
 """
 single-equation approximation of the vapour and heat diffusion problem
+as given in eq. 3.11 in [Mason 1971](https://archive.org/details/physicsofclouds0000maso)
+(see also discussion of the ventilation effect on page 125)
 """
 
 
-class MaxwellMason:  # pylint: disable=too-few-public-methods
+class Mason1971:  # pylint: disable=too-few-public-methods
     def __init__(self, _):
         pass
 
@@ -13,8 +15,6 @@ def r_dr_dt(const, RH_eq, T, RH, lv, pvs, D, K, ventilation_factor):
         return (
             ventilation_factor
             * (RH - RH_eq)
-            / (
-                const.rho_w * const.Rv * T / D / pvs
-                + const.rho_w * lv / K / T * (lv / const.Rv / T - 1)
-            )
+            / const.rho_w
+            / (const.Rv * T / D / pvs + lv / K / T * (lv / T / const.Rv - 1))
         )

tests/unit_tests/physics/test_drop_growth.py

@@ -0,0 +1,76 @@
+""" tests for drop growth formulae """
+
+import pytest
+from matplotlib import pyplot
+import numpy as np
+
+from PySDM.formulae import _choices, Formulae
+from PySDM.physics import drop_growth, constants_defaults, si, in_unit
+from PySDM.physics.dimensional_analysis import DimensionalAnalysis
+
+
+class TestDropGrowth:
+    @staticmethod
+    @pytest.mark.parametrize("paper", _choices(drop_growth))
+    def test_unit(paper):
+        """checks dimensionality of the returned value"""
+        with DimensionalAnalysis():
+            # arrange
+            formulae = Formulae(drop_growth=paper)
+
+            # act
+            r_dr_dt = formulae.drop_growth.r_dr_dt(
+                RH_eq=1,
+                T=constants_defaults.T_tri,
+                RH=1.05,
+                lv=constants_defaults.l_tri,
+                pvs=constants_defaults.p_tri,
+                D=constants_defaults.D0,
+                K=constants_defaults.K0,
+                ventilation_factor=1,
+            )
+
+            # assert
+            assert r_dr_dt.check("[area]/[time]")
+
+    @staticmethod
+    def test_mason_1971_vs_1951_difference_vs_temperature(plot=False):
+        """checks the relative difference between Mason's 1951 and 1971 formulae
+        for a range of temperatures"""
+        # arrange
+        temperatures = constants_defaults.T0 + np.linspace(-10, 40) * si.K
+        papers = ("Mason1951", "Mason1971")
+
+        # act
+        formulae = {paper: Formulae(drop_growth=paper) for paper in papers}
+        r_dr_dt = {
+            paper: formulae[paper].drop_growth.r_dr_dt(
+                RH_eq=1,
+                T=temperatures,
+                RH=1.05,
+                lv=constants_defaults.l_tri,
+                pvs=constants_defaults.p_tri,
+                D=constants_defaults.D0,
+                K=constants_defaults.K0,
+                ventilation_factor=1,
+            )
+            for paper in papers
+        }
+        relative_error = r_dr_dt["Mason1971"] / r_dr_dt["Mason1951"] - 1
+
+        # plot
+        pyplot.plot(temperatures, in_unit(relative_error, constants_defaults.PER_CENT))
+        pyplot.title("")
+        pyplot.xlabel("temperature [K]")
+        pyplot.ylabel("r dr/dt relative difference (1971 vs. 1951) [%]")
+        pyplot.grid()
+        pyplot.legend()
+        if plot:
+            pyplot.show()
+        else:
+            pyplot.clf()
+
+        # assert
+        (relative_error < 0.03).all()
+        (relative_error > 0.02).all()
+        (np.diff(relative_error) < 0).all()