Fixes to Straub fragmentation function and sampling methods, addition of LowList1982 fragmentation function and coalescence efficiency (#1023)

Co-authored-by: Sylwester Arabas <[email protected]>
Co-authored-by: Oleksii Bulenok <[email protected]>
Co-authored-by: de Jong <[email protected]>
Co-authored-by: Oleksii Bulenok <[email protected]>

Author: 4 people

PySDM/backends/impl_numba/methods/collisions_methods.py

@@ -521,23 +521,24 @@ def __init__(self):
             param_names=(
                 "n_fragment",
                 "frag_size",
-                "v_max",
                 "x_plus_y",
                 "vmin",
                 "nfmax",
                 "nfmax_is_not_none",
             ),
             name_iter="i",
             body="""
-            frag_size[i] = min(frag_size[i], v_max[i]);
-            frag_size[i] = max(frag_size[i], vmin);
+            frag_size[i] = min(frag_size[i], x_plus_y[i])
 
             if (nfmax_is_not_none) {
                 if (x_plus_y[i] / frag_size[i] > nfmax) {
                     frag_size[i] = x_plus_y[i] / nfmax;
                 }
             }
-            if (frag_size[i] == 0.0) {
+            else if (frag_size[i] < vmin) {
+                frag_size[i] = x_plus_y[i];
+            }
+            else if (frag_size[i] == 0.0) {
                 frag_size[i] = x_plus_y[i];
             }
             n_fragment[i] = x_plus_y[i] / frag_size[i];
@@ -549,10 +550,8 @@ def __init__(self):
                 param_names=("mu", "sigma", "frag_size", "rand"),
                 name_iter="i",
                 body=f"""
-                frag_size[i] = {self.formulae.fragmentation_function.frag_size.c_inline(
-                    mu="mu",
-                    sigma="sigma",
-                    rand="rand[i]"
+                frag_size[i] = mu + sigma * {self.formulae.trivia.erfinv_approx.c_inline(
+                    c="rand[i]"
                 )};
                 """.replace(
                     "real_type", self._get_c_type()
@@ -610,6 +609,20 @@ def __init__(self):
                 Nrt[i] = Nr1[i] + Nr2[i] + Nr3[i] + Nr4[i];
         """
 
+        self.__straub_mass_remainder = """
+                Nr1[i] = Nr1[i] * exp(3 * mu1 + 9 * pow(sigma1, 2.0) / 2);
+                Nr2[i] = Nr2[i] * (pow(mu2, 3.0) + 3 * mu2 * pow(sigma2, 2.0));
+                Nr3[i] = Nr3[i] * (pow(mu3, 3.0) + 3 * mu3 * pow(sigma3, 2.0));
+                Nr4[i] = v_max[i] * 6.0 / 3.141592654 + pow(ds[i], 3.0) - Nr1[i] - Nr2[i] - Nr3[i];
+                if Nr4[i] <= 0.0 {
+                    d34[i] = 0;
+                    Nr4[i] = 0;
+                }
+                else {
+                    d34[i] = exp(log(Nr4[i]) / 3);
+                }
+        """
+
         if self.formulae.fragmentation_function.__name__ == "Straub2010Nf":
             self.__straub_fragmentation_body = trtc.For(
                 param_names=(
@@ -624,41 +637,44 @@ def __init__(self):
                     "Nr3",
                     "Nr4",
                     "Nrt",
+                    "d34",
                 ),
                 name_iter="i",
                 body=f"""
                 {self.__straub_Nr_body}
+                auto sigma1 = {self.formulae.fragmentation_function.params_sigma1.c_inline(CW="CW[i]")};
+                auto mu1 = {self.formulae.fragmentation_function.params_mu1.c_inline(sigma1="sigma1")};
+                auto sigma2 = {self.formulae.fragmentation_function.params_sigma2.c_inline(CW="CW[i]")};
+                auto mu2 = {self.formulae.fragmentation_function.params_mu2.c_inline(ds="ds[i]")};
+                auto sigma3 = {self.formulae.fragmentation_function.params_sigma3.c_inline(CW="CW[i]")};
+                auto mu3 = {self.formulae.fragmentation_function.params_mu3.c_inline(ds="ds[i]")};
+                {self.__straub_mass_remainder}
+                Nrt[i] = Nr1[i] + Nr2[i] + Nr3[i] + Nr4[i]
 
                 if (rand[i] < Nr1[i] / Nrt[i]) {{
-                    auto sigma1 = {self.formulae.fragmentation_function.sigma1.c_inline(CW="CW[i]")};
-                    frag_size[i] = {self.formulae.fragmentation_function.p1.c_inline(
-                        sigma1="sigma1",
-                        rand="rand[i] * Nrt[i] / Nr1[i]"
+                    auto X = rand[i] * Nrt[i] / Nr1[i];
+                    auto lnarg = mu1 + sqrt(2.0) * sigma1 * {self.formulae.trivia.erfinv_approx.c_inline(
+                        c="X"
                     )};
+                    frag_size[i] = exp(lnarg);
                 }}
                 else if (rand[i] < (Nr2[i] + Nr1[i]) / Nrt[i]) {{
-                    frag_size[i] = {self.formulae.fragmentation_function.p2.c_inline(
-                        CW="CW[i]",
-                        rand="(rand[i] * Nrt[i] - Nr1[i]) / (Nr2[i] - Nr1[i])"
+                    auto X = (rand[i] * Nrt[i] - Nr1[i]) / Nr2[i];
+                    frag_size[i] = mu2 + sqrt(2.0) * sigma2 * {self.formulae.trivia.erfinv_approx.c_inline(
+                        c="X"
                     )};
                 }}
                 else if (rand[i] < (Nr3[i] + Nr2[i] + Nr1[i]) / Nrt[i]) {{
-                    frag_size[i] = {self.formulae.fragmentation_function.p3.c_inline(
-                        CW="CW[i]",
-                        ds="ds[i]",
-                        rand="(rand[i] * Nrt[i] - Nr2[i]) / (Nr3[i] - Nr2[i])"
+                    auto X = (rand[i] * Nrt[i] - Nr1[i] - Nr2[i]) / Nr3[i];
+                    frag_size[i] = mu3 + sqrt(2.0) * sigma3 * {self.formulae.trivia.erfinv_approx.c_inline(
+                        c="X"
                     )};
                 }}
                 else {{
-                    frag_size[i] = {self.formulae.fragmentation_function.p4.c_inline(
-                        CW="CW[i]",
-                        ds="ds[i]",
-                        v_max="v_max[i]",
-                        Nr1="Nr1[i]",
-                        Nr2="Nr2[i]",
-                        Nr3="Nr3[i]"
-                    )};
+                    frag_size[i] = d34[i];
                 }}
+
+                frag_size[i] = pow(frag_size[i], 3.0) * 3.141592654 / 6.0
                 """.replace(
                     "real_type", self._get_c_type()
                 ),
@@ -910,7 +926,6 @@ def exp_fragmentation(
         n_fragment,
         scale,
         frag_size,
-        v_max,
         x_plus_y,
         rand,
         vmin,
@@ -932,7 +947,6 @@ def exp_fragmentation(
             args=(
                 n_fragment.data,
                 frag_size.data,
-                v_max.data,
                 x_plus_y.data,
                 self._get_floating_point(vmin),
                 self._get_floating_point(nfmax if nfmax else -1),
@@ -941,7 +955,7 @@ def exp_fragmentation(
         )
 
     def gauss_fragmentation(
-        self, *, n_fragment, mu, sigma, frag_size, v_max, x_plus_y, rand, vmin, nfmax
+        self, *, n_fragment, mu, sigma, frag_size, x_plus_y, rand, vmin, nfmax
     ):
         self.__gauss_fragmentation_body.launch_n(
             n=len(frag_size),
@@ -958,7 +972,6 @@ def gauss_fragmentation(
             args=(
                 n_fragment.data,
                 frag_size.data,
-                v_max.data,
                 x_plus_y.data,
                 self._get_floating_point(vmin),
                 self._get_floating_point(nfmax if nfmax else -1),
@@ -967,7 +980,7 @@ def gauss_fragmentation(
         )
 
     def slams_fragmentation(
-        self, n_fragment, frag_size, v_max, x_plus_y, probs, rand, vmin, nfmax
+        self, n_fragment, frag_size, x_plus_y, probs, rand, vmin, nfmax
     ):  # pylint: disable=too-many-arguments
         self.__slams_fragmentation_body.launch_n(
             n=(len(n_fragment)),
@@ -985,7 +998,6 @@ def slams_fragmentation(
             args=(
                 n_fragment.data,
                 frag_size.data,
-                v_max.data,
                 x_plus_y.data,
                 self._get_floating_point(vmin),
                 self._get_floating_point(nfmax if nfmax else -1),
@@ -999,7 +1011,6 @@ def feingold1988_fragmentation(
         n_fragment,
         scale,
         frag_size,
-        v_max,
         x_plus_y,
         rand,
         fragtol,
@@ -1022,7 +1033,6 @@ def feingold1988_fragmentation(
             args=(
                 n_fragment.data,
                 frag_size.data,
-                v_max.data,
                 x_plus_y.data,
                 self._get_floating_point(vmin),
                 self._get_floating_point(nfmax if nfmax else -1),
@@ -1049,6 +1059,7 @@ def straub_fragmentation(
         Nr3,
         Nr4,
         Nrt,
+        d34,
     ):
         self.__straub_fragmentation_body.launch_n(
             n=len(frag_size),
@@ -1064,6 +1075,7 @@ def straub_fragmentation(
                 Nr3.data,
                 Nr4.data,
                 Nrt.data,
+                d34.data,
             ),
         )
 
@@ -1072,7 +1084,6 @@ def straub_fragmentation(
             args=(
                 n_fragment.data,
                 frag_size.data,
-                v_max.data,
                 x_plus_y.data,
                 self._get_floating_point(vmin),
                 self._get_floating_point(nfmax if nfmax else -1),

PySDM/backends/impl_thrust_rtc/methods/collisions_methods.py

@@ -190,7 +190,7 @@ def to_numba(name, args, iter_var, body):
         f"""
 def make(self):
     import numpy as np
-    from numpy import floor, ceil, exp, log, power, sqrt, arctanh
+    from numpy import floor, ceil, exp, log, power, sqrt, arctanh, sinh, arcsinh
     import numba
 
     @numba.njit(parallel=False, {JIT_OPTS})

PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py

@@ -6,5 +6,6 @@
 from .exponential import ExponFrag
 from .feingold1988 import Feingold1988Frag
 from .gaussian import Gaussian
+from .lowlist82 import LowList1982Nf
 from .slams import SLAMS
 from .straub2010 import Straub2010Nf

PySDM/dynamics/collisions/breakup_fragmentations/__init__.py

@@ -10,28 +10,22 @@ def __init__(self, scale, vmin=0.0, nfmax=None):
         self.scale = scale
         self.vmin = vmin
         self.nfmax = nfmax
-        self.max_size = None
         self.sum_of_volumes = None
 
     def register(self, builder):
         self.particulator = builder.particulator
-        self.max_size = self.particulator.PairwiseStorage.empty(
-            self.particulator.n_sd // 2, dtype=float
-        )
         self.sum_of_volumes = self.particulator.PairwiseStorage.empty(
             self.particulator.n_sd // 2, dtype=float
         )
 
     def __call__(self, nf, frag_size, u01, is_first_in_pair):
-        self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair)
         self.sum_of_volumes.sum(
             self.particulator.attributes["volume"], is_first_in_pair
         )
         self.particulator.backend.exp_fragmentation(
             n_fragment=nf,
             scale=self.scale,
             frag_size=frag_size,
-            v_max=self.max_size,
             x_plus_y=self.sum_of_volumes,
             rand=u01,
             vmin=self.vmin,

PySDM/dynamics/collisions/breakup_fragmentations/exponential.py

@@ -11,29 +11,23 @@ def __init__(self, scale, fragtol=1e-3, vmin=0.0, nfmax=None):
         self.fragtol = fragtol
         self.vmin = vmin
         self.nfmax = nfmax
-        self.max_size = None
         self.sum_of_volumes = None
 
     def register(self, builder):
         self.particulator = builder.particulator
         builder.request_attribute("volume")
-        self.max_size = self.particulator.PairwiseStorage.empty(
-            self.particulator.n_sd // 2, dtype=float
-        )
         self.sum_of_volumes = self.particulator.PairwiseStorage.empty(
             self.particulator.n_sd // 2, dtype=float
         )
 
     def __call__(self, nf, frag_size, u01, is_first_in_pair):
-        self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair)
         self.sum_of_volumes.sum(
             self.particulator.attributes["volume"], is_first_in_pair
         )
         self.particulator.backend.feingold1988_fragmentation(
             n_fragment=nf,
             scale=self.scale,
             frag_size=frag_size,
-            v_max=self.max_size,
             x_plus_y=self.sum_of_volumes,
             rand=u01,
             fragtol=self.fragtol,

PySDM/dynamics/collisions/breakup_fragmentations/feingold1988.py

@@ -10,20 +10,15 @@ def __init__(self, mu, sigma, vmin=0.0, nfmax=None):
         self.sigma = sigma
         self.vmin = vmin
         self.nfmax = nfmax
-        self.max_size = None
         self.sum_of_volumes = None
 
     def register(self, builder):
         self.particulator = builder.particulator
-        self.max_size = self.particulator.PairwiseStorage.empty(
-            self.particulator.n_sd // 2, dtype=float
-        )
         self.sum_of_volumes = self.particulator.PairwiseStorage.empty(
             self.particulator.n_sd // 2, dtype=float
         )
 
     def __call__(self, nf, frag_size, u01, is_first_in_pair):
-        self.max_size.max(self.particulator.attributes["volume"], is_first_in_pair)
         self.sum_of_volumes.sum(
             self.particulator.attributes["volume"], is_first_in_pair
         )
@@ -32,7 +27,6 @@ def __call__(self, nf, frag_size, u01, is_first_in_pair):
             mu=self.mu,
             sigma=self.sigma,
             frag_size=frag_size,
-            v_max=self.max_size,
             x_plus_y=self.sum_of_volumes,
             rand=u01,
             vmin=self.vmin,