Jensen & Nugent 2017 GCCN example (incl. Grabowski et al. diffusion t…

…hermics) (#1267)

Co-authored-by: Sanket Bhiogade <sbhiogade@agh.edu.pl>
Co-authored-by: Sylwester Arabas <sylwester.arabas@agh.edu.pl>

PySDM/initialisation/sampling/spectral_sampling.py

@@ -55,7 +55,8 @@ def _sample(self, grid, spectrum):
         diff = abs(1 - np.sum(y_float) / spectrum.norm_factor)
         if diff > self.error_threshold:
             raise ValueError(
-                f"{100*diff}% error in total real-droplet number due to sampling"
+                f"{diff * 100:.3g}% error in total real-droplet number due to sampling "
+                f"({len(x)} samples)"
             )
 
         return x, y_float

PySDM/physics/constants_defaults.py

@@ -321,6 +321,16 @@
 """ [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.15 / si.K
+diffussion_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
+
 
 def compute_derived_values(c: dict):
     c["eps"] = c["Mv"] / c["Md"]

PySDM/physics/diffusion_kinetics/__init__.py

@@ -5,3 +5,5 @@
 from .fuchs_sutugin import FuchsSutugin
 from .lowe_et_al_2019 import LoweEtAl2019
 from .neglect import Neglect
+from .grabowski_et_al_2011 import GrabowskiEtAl2011
+from .jensen_and_nugent_2017 import JensenAndNugent2017

PySDM/physics/diffusion_kinetics/grabowski_et_al_2011.py

@@ -0,0 +1,9 @@
+"""
+as in [Grabowski et al. (2011)](https://doi.org/10.1016/j.atmosres.2010.10.020)
+"""
+
+from .pruppacher_and_klett_2005 import PruppacherKlett
+
+
+class GrabowskiEtAl2011(PruppacherKlett):
+    pass

PySDM/physics/diffusion_kinetics/jensen_and_nugent_2017.py

@@ -0,0 +1,14 @@
+"""
+as in [Jensen and Nugent 2017](https://doi.org/10.1175/JAS-D-15-0370.1)
+(which refers to [Grabowski et al. (2011)](https://doi.org/10.1016/j.atmosres.2010.10.020)
+but uses different gas constant, typo?)
+"""
+
+import numpy as np
+from .grabowski_et_al_2011 import GrabowskiEtAl2011
+
+
+class JensenAndNugent2017(GrabowskiEtAl2011):
+    @staticmethod
+    def lambdaD(const, D, T):
+        return D / np.sqrt(2 * const.Rd * T)

PySDM/physics/diffusion_kinetics/pruppacher_and_klett_2005.py

@@ -19,3 +19,11 @@ def D(const, D, r, lmbd):
         return D / (
             (r / (r + const.dv_pk05)) + 2 * np.sqrt(const.PI) * lmbd / r / const.MAC
         )
+
+    @staticmethod
+    def lambdaK(_, T, p):  # pylint: disable=unused-argument
+        return -1
+
+    @staticmethod
+    def K(const, K, r, lmbd):  # pylint: disable=unused-argument
+        return K  # TODO #1266

PySDM/physics/diffusion_thermics/__init__.py

@@ -5,3 +5,4 @@
 from .lowe_et_al_2019 import LoweEtAl2019
 from .neglect import Neglect
 from .tracy_welch_porter import TracyWelchPorter
+from .grabowski_et_al_2011 import GrabowskiEtAl2011

PySDM/physics/diffusion_thermics/grabowski_et_al_2011.py

@@ -0,0 +1,23 @@
+"""
+as in [Grabowski et al. (2011)](https://doi.org/10.1016/j.atmosres.2010.10.020)
+"""
+
+
+class GrabowskiEtAl2011:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def D(const, T, p):  # pylint: disable=unused-argument
+        return const.diffussion_thermics_D_G11_A * (
+            const.diffussion_thermics_D_G11_B * T + const.diffussion_thermics_D_G11_C
+        )
+
+    @staticmethod
+    def K(const, T, p):  # pylint: disable=unused-argument
+        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
+        )

PySDM/physics/impl/fake_unit_registry.py

@@ -26,6 +26,7 @@ def __init__(self, si):
                 "litre",
                 "hour",
                 "newton",
+                "watt",
             ):
                 self.__setattr__(prefix + unit, _fake(si.__getattr__(prefix + unit)))
                 self.__setattr__(
@@ -48,5 +49,6 @@ def __init__(self, si):
                 "h",
                 "bar",
                 "N",
+                "W",
             ):
                 self.__setattr__(prefix + unit, _fake(si.__getattr__(prefix + unit)))

PySDM/physics/trivia.py

@@ -125,3 +125,7 @@ def isotopic_enrichment_to_delta_SMOW(E, delta_0_SMOW):
     @staticmethod
     def mixing_ratio_to_specific_content(mixing_ratio):
         return mixing_ratio / (1 + mixing_ratio)
+
+    @staticmethod
+    def dn_dlogr(r, dn_dr):
+        return np.log(10) * r * dn_dr

PySDM/products/condensation/peak_supersaturation.py

@@ -19,6 +19,10 @@ def __init__(self, unit="dimensionless", name=None):
     def register(self, builder):
         super().register(builder)
         self.particulator.observers.append(self)
+
+        assert (
+            "Condensation" in self.particulator.dynamics
+        ), "It seems the Condensation dynamic was not added when building particulator"
         self.condensation = self.particulator.dynamics["Condensation"]
         self.RH_max = np.full_like(self.buffer, np.nan)