Add the Long 1974 piecewise collision kernel

PySDM/backends/impl_numba/storage.py

@@ -1,3 +1,4 @@
+# pylint: disable=too-many-public-methods
 """
 CPU Numpy-based implementation of Storage class
 """
@@ -181,6 +182,13 @@ def divide_if_not_zero(self, divisor):
         impl.divide_if_not_zero(self.data, divisor.data)
         return self
 
+    def where(self, condition, true_value, false_value):
+        impl.where(self.data, condition.data, true_value.data, false_value.data)
+        return self
+
+    def isless(self, comparison, value):
+        impl.isless(self.data, comparison.data, value)
+
     def sum(self, arg_a, arg_b):
         impl.sum_out_of_place(self.data, arg_a.data, arg_b.data)
         return self

PySDM/backends/impl_numba/storage_impl.py

@@ -42,6 +42,24 @@ def divide_if_not_zero(output, divisor):
             output[i] /= divisor[i]
 
 
+@numba.njit(**conf.JIT_FLAGS)
+def where(output, condition, true_value, false_value):
+    for i in numba.prange(output.shape[0]):  # pylint: disable=not-an-iterable
+        if condition[i]:
+            output[i] = true_value[i]
+        else:
+            output[i] = false_value[i]
+
+
+@numba.njit(**conf.JIT_FLAGS)
+def isless(output, comparison, value):
+    for i in numba.prange(output.shape[0]):  # pylint: disable=not-an-iterable
+        if comparison[i] < value:
+            output[i] = True
+        else:
+            output[i] = False
+
+
 @numba.njit(**conf.JIT_FLAGS)
 def floor(output):
     output[:] = np.floor(output)

PySDM/backends/impl_thrust_rtc/storage.py

@@ -306,6 +306,47 @@ def divide_if_not_zero(output, divisor):
                 n=(output.shape[0]), args=(Impl.thrust((output, divisor)))
             )
 
+        __where_body = trtc.For(
+            param_names=("output", "condition", "t_value", "f_value"),
+            name_iter="i",
+            body="""
+            if (condition[i]) {
+                output[i] = t_value[i]
+            }
+            else {
+                output[i] = f_value[i]
+            }
+            """,
+        )
+
+        @staticmethod
+        @nice_thrust(**NICE_THRUST_FLAGS)
+        def where(output, condition, t_value, f_value):
+            Impl.__where_body.launch_n(
+                n=(output.shape[0]),
+                args=(Impl.thrust((output, condition, t_value, f_value))),
+            )
+
+        __isless_body = trtc.For(
+            param_names=("output", "comparison", "value"),
+            name_iter="i",
+            body="""
+            if (comparison[i] < value) {
+                output[i] = true
+            }
+            else {
+                output[i] = false
+            }
+            """,
+        )
+
+        @staticmethod
+        @nice_thrust(**NICE_THRUST_FLAGS)
+        def isless(output, comparison, value):
+            Impl.__isless_body.launch_n(
+                n=(output.shape[0]), args=(Impl.thrust((output, comparison, value)))
+            )
+
         __exp_body = trtc.For(
             ("output",),
             "i",
@@ -550,6 +591,14 @@ def divide_if_not_zero(self, divisor):
             Impl.divide_if_not_zero(self, divisor)
             return self
 
+        def where(self, condition, t_value, f_value):
+            Impl.where(self, condition, t_value, f_value)
+            return self
+
+        def isless(self, comparison, value):
+            Impl.isless(self, comparison, value)
+            return self
+
         def fill(self, other):
             if isinstance(other, Storage):
                 trtc.Copy(other.data, self.data)

PySDM/dynamics/collisions/collision_kernels/__init__.py

@@ -12,3 +12,4 @@
 from .hydrodynamic import Hydrodynamic
 from .linear import Linear
 from .simple_geometric import SimpleGeometric
+from .long1974 import Long1974

PySDM/dynamics/collisions/collision_kernels/long1974.py

@@ -0,0 +1,74 @@
+"""
+piecewise kernel from Eq (11) in
+[Long 1974](https://doi.org/10.1175/1520-0469%281974%29031%3C1040%3ASTTDCE%3E2.0.CO%3B2)
+
+Default parameters are:
+    lin_coeff = 5.78e3 / s
+    sq_coeff = 9.44e15 / m^3 / s
+    r_thresh = 5e-5 m
+"""
+
+
+class Long1974:
+    def __init__(self, lin_coeff=5.78e3, sq_coeff=9.44e15, r_thres=5e-5):
+        self.lc = lin_coeff
+        self.sc = sq_coeff
+        self.rt = r_thres
+        self.particulator = None
+        self.largeR = None
+        self.arrays = {}
+
+    def register(self, builder):
+        self.particulator = builder.particulator
+        builder.request_attribute("volume")
+        builder.request_attribute("radius")
+        for key in (
+            "r_lg",
+            "v_lg",
+            "v_sm",
+            "v_ratio",
+            "tmp",
+            "tmp1",
+        ):
+            self.arrays[key] = self.particulator.PairwiseStorage.empty(
+                self.particulator.n_sd // 2, dtype=float
+            )
+            self.arrays["condition"] = self.particulator.PairwiseStorage.empty(
+                self.particulator.n_sd // 2, dtype=bool
+            )
+
+    def __call__(self, output, is_first_in_pair):
+        # get smaller and larger radii, volume
+        self.arrays["r_lg"].max(
+            self.particulator.attributes["radius"], is_first_in_pair
+        )
+        self.arrays["v_lg"].max(
+            self.particulator.attributes["volume"], is_first_in_pair
+        )
+        self.arrays["v_sm"].min(
+            self.particulator.attributes["volume"], is_first_in_pair
+        )
+
+        # compute volume ratio
+        self.arrays["v_ratio"].fill(self.arrays["v_sm"])
+        self.arrays["v_ratio"].divide_if_not_zero(self.arrays["v_lg"])
+
+        # compute small radius limit
+        self.arrays["tmp1"].fill(self.arrays["v_ratio"])
+        self.arrays["tmp1"] **= 2.0
+        self.arrays["tmp1"] += 1.0
+        self.arrays["tmp"].fill(self.arrays["v_lg"])
+        self.arrays["tmp"] **= 2.0
+        self.arrays["tmp1"] *= self.arrays["tmp"]
+        self.arrays["tmp1"] *= self.sc
+
+        # compute large radius (linear) limit
+        self.arrays["v_ratio"] += 1.0
+        self.arrays["v_ratio"] *= self.arrays["v_lg"]
+        self.arrays["v_ratio"] *= self.lc
+
+        # apply piecewise
+        self.arrays["condition"].isless(self.arrays["r_lg"], self.rt)
+        output.where(
+            self.arrays["condition"], self.arrays["tmp1"], self.arrays["v_ratio"]
+        )

docs/bibliography.json

@@ -1012,5 +1012,14 @@
     ],
     "label": "Hayes 2004",
     "title": "An Introduction to Isotopic Calculations"
+  },
+  "https://doi.org/10.1175/1520-0469%281974%29031%3C1040%3ASTTDCE%3E2.0.CO%3B2" : {
+    "usages": [
+      "PySDM/dynamics/collisions/collision_kernels/long1974.py",
+      "examples/PySDM_examples/Long_1974/__init__.py",
+      "examples/PySDM_examples/Long_1974/fig10_11_13_14.ipynb"
+    ],
+    "label": "Long 1974 (J. Atmos. Sci. 31)",
+    "title": "Solutions to the Droplet Collection Equation for Polynomial Kernels"
   }
 }

examples/PySDM_examples/Long_1974/__init__.py

@@ -0,0 +1,8 @@
+# pylint: disable=invalid-name
+"""
+piecewise collision-coalescence example from
+[Long 1974 (JAS)](https://doi.org/10.1175/1520-0469%281974%29031%3C1040%3ASTTDCE%3E2.0.CO%3B2)
+"""
+
+from .settings import Settings
+from .simulation import Simulation