collision dynamic: kernel choice move to formulae and setup for collisions including ice

PySDM/dynamics/collisions/collision.py

@@ -18,6 +18,7 @@
 from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb
 from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN
 from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc
+from PySDM.dynamics.collisions import collision_kernels
 from PySDM.dynamics.impl.random_generator_optimizer import RandomGeneratorOptimizer
 from PySDM.dynamics.impl.random_generator_optimizer_nopair import (
     RandomGeneratorOptimizerNoPair,
@@ -42,7 +43,6 @@ class Collision:  # pylint: disable=too-many-instance-attributes
     def __init__(
         self,
         *,
-        collision_kernel,
         coalescence_efficiency,
         breakup_efficiency,
         fragmentation_function,
@@ -63,7 +63,6 @@ def __init__(
         self.warn_overflows = warn_overflows
         self.max_multiplicity = DEFAULTS.max_multiplicity
 
-        self.collision_kernel = collision_kernel
         self.compute_coalescence_efficiency = coalescence_efficiency
         self.compute_breakup_efficiency = breakup_efficiency
         self.compute_number_of_fragments = fragmentation_function
@@ -98,6 +97,8 @@ def __init__(
         self.breakup_rate = None
         self.breakup_rate_deficit = None
 
+        self.collision_kernel = None
+
     def register(self, builder):
         self.particulator = builder.particulator
         rnd_args = {
@@ -136,6 +137,10 @@ def register(self, builder):
         self.stats_dt_min.fill(np.nan)
 
         self.rnd_opt_coll.register(builder)
+        self.collision_kernel = getattr(
+            collision_kernels,
+            self.particulator.formulae.collision_kernel_liquid_liquid.__name__,
+        )()
         self.collision_kernel.register(builder)
 
         if self.croupier is None:
@@ -295,7 +300,6 @@ class Coalescence(Collision):
     def __init__(
         self,
         *,
-        collision_kernel,
         coalescence_efficiency=ConstEc(Ec=1),
         croupier=None,
         optimized_random=False,
@@ -306,7 +310,6 @@ def __init__(
         breakup_efficiency = ConstEb(Eb=0)
         fragmentation_function = AlwaysN(n=1)
         super().__init__(
-            collision_kernel=collision_kernel,
             coalescence_efficiency=coalescence_efficiency,
             breakup_efficiency=breakup_efficiency,
             fragmentation_function=fragmentation_function,
@@ -324,7 +327,6 @@ class Breakup(Collision):
     def __init__(
         self,
         *,
-        collision_kernel,
         fragmentation_function,
         croupier=None,
         optimized_random=False,
@@ -336,7 +338,6 @@ def __init__(
         coalescence_efficiency = ConstEc(Ec=0.0)
         breakup_efficiency = ConstEb(Eb=1.0)
         super().__init__(
-            collision_kernel=collision_kernel,
             coalescence_efficiency=coalescence_efficiency,
             breakup_efficiency=breakup_efficiency,
             fragmentation_function=fragmentation_function,

PySDM/dynamics/collisions/collision_kernels/__init__.py

@@ -13,3 +13,4 @@
 from .linear import Linear
 from .simple_geometric import SimpleGeometric
 from .long1974 import Long1974
+from .neglect import Neglect

PySDM/dynamics/collisions/collision_kernels/constantK.py

@@ -4,12 +4,11 @@
 
 
 class ConstantK:
-    def __init__(self, a):
-        self.a = a
+    def __init__(self):
         self.particulator = None
 
     def __call__(self, output, is_first_in_pair):
-        output.fill(self.a)
+        output.fill(self.particulator.formulae.constants.CONSTANTK_a)
 
     def register(self, builder):
         self.particulator = builder.particulator

PySDM/dynamics/collisions/collision_kernels/geometric.py

@@ -7,15 +7,13 @@
 
 
 class Geometric(Gravitational):
-    def __init__(self, collection_efficiency=1.0, x="volume"):
-        super().__init__()
-        self.collection_efficiency = collection_efficiency
-        self.x = x
-
     def __call__(self, output, is_first_in_pair):
         output.sum(self.particulator.attributes["radius"], is_first_in_pair)
         output **= 2
-        output *= const.PI * self.collection_efficiency
+        output *= (
+            const.PI
+            * self.particulator.formulae.constants.GEOMETRIC_collection_efficiency
+        )
         self.pair_tmp.distance(
             self.particulator.attributes["relative fall velocity"], is_first_in_pair
         )

PySDM/dynamics/collisions/collision_kernels/golovin.py

@@ -4,43 +4,15 @@
 [Golovin 1963](http://mi.mathnet.ru/dan27630))
 """
 
-import numpy as np
-from scipy import special
-
 
 class Golovin:
-    def __init__(self, b):
-        self.b = b
+    def __init__(self):
         self.particulator = None
 
     def __call__(self, output, is_first_in_pair):
         output.sum(self.particulator.attributes["volume"], is_first_in_pair)
-        output *= self.b
+        output *= self.particulator.formulae.constants.GOLOVIN_b
 
     def register(self, builder):
         self.particulator = builder.particulator
         builder.request_attribute("volume")
-
-    def analytic_solution(self, x, t, x_0, N_0):
-        tau = 1 - np.exp(-N_0 * self.b * x_0 * t)
-
-        if isinstance(x, np.ndarray):
-            func = np.vectorize(
-                lambda i: Golovin.analytic_solution_helper(x[int(i)], tau, x_0)
-            )
-            result = np.fromfunction(func, x.shape, dtype=float)
-            return result
-
-        return Golovin.analytic_solution_helper(x, tau, x_0)
-
-    @staticmethod
-    def analytic_solution_helper(x, tau, x_0):
-        sqrt_tau = np.sqrt(tau)
-        result = float(
-            (1 - tau)
-            * 1
-            / (x * np.sqrt(tau))
-            * special.ive(1, 2 * x / x_0 * sqrt_tau)  # pylint: disable=no-member
-            * np.exp(-(1 + tau - 2 * sqrt_tau) * x / x_0)
-        )
-        return result

PySDM/dynamics/collisions/collision_kernels/linear.py

@@ -4,15 +4,13 @@
 
 
 class Linear:
-    def __init__(self, a, b):
-        self.a = a
-        self.b = b
+    def __init__(self):
         self.particulator = None
 
     def __call__(self, output, is_first_in_pair):
         output.sum_pair(self.particulator.attributes["volume"], is_first_in_pair)
-        output *= self.b
-        output += self.a
+        output *= self.particulator.formulae.constants.LINEAR_b
+        output += self.particulator.formulae.constants.LINEAR_a
 
     def register(self, builder):
         self.particulator = builder.particulator

PySDM/dynamics/collisions/collision_kernels/neglect.py

@@ -0,0 +1,14 @@
+"""
+kernel for disabling collision
+"""
+
+
+class Neglect:
+    def __init__(self):
+        self.particulator = None
+
+    def __call__(self, output, is_first_in_pair):
+        output.fill(0)
+
+    def register(self, builder):
+        self.particulator = builder.particulator

PySDM/dynamics/collisions/collision_kernels/simple_geometric.py

@@ -4,10 +4,9 @@
 
 
 class SimpleGeometric:
-    def __init__(self, C):
+    def __init__(self):
         self.particulator = None
         self.pair_tmp = None
-        self.C = C
 
     def register(self, builder):
         self.particulator = builder.particulator
@@ -18,7 +17,7 @@ def register(self, builder):
         )
 
     def __call__(self, output, is_first_in_pair):
-        output[:] = self.C
+        output[:] = self.particulator.formulae.constants.SIMPLE_GEOMETRIC_c
         self.pair_tmp.sum(self.particulator.attributes["radius"], is_first_in_pair)
         self.pair_tmp **= 2
         output *= self.pair_tmp

PySDM/formulae.py

@@ -55,6 +55,7 @@ def __init__(  # pylint: disable=too-many-locals,too-many-arguments
         freezing_temperature_spectrum: str = "Null",
         heterogeneous_ice_nucleation_rate: str = "Null",
         homogeneous_ice_nucleation_rate: str = "Null",
+        collision_kernel: str = "Neglect+Neglect+Neglect",
         fragmentation_function: str = "AlwaysN",
         isotope_equilibrium_fractionation_factors: str = "Null",
         isotope_kinetic_fractionation_factors: str = "Null",
@@ -96,6 +97,9 @@ def __init__(  # pylint: disable=too-many-locals,too-many-arguments
         self.freezing_temperature_spectrum = freezing_temperature_spectrum
         self.heterogeneous_ice_nucleation_rate = heterogeneous_ice_nucleation_rate
         self.homogeneous_ice_nucleation_rate = homogeneous_ice_nucleation_rate
+        self.collision_kernel_liquid_liquid = collision_kernel.split("+")[0]
+        self.collision_kernel_liquid_ice = collision_kernel.split("+")[1]
+        self.collision_kernel_ice_ice = collision_kernel.split("+")[2]
         self.fragmentation_function = fragmentation_function
         self.isotope_equilibrium_fractionation_factors = (
             isotope_equilibrium_fractionation_factors

PySDM/physics/__init__.py

@@ -54,5 +54,6 @@
     terminal_velocity,
     terminal_velocity_ice,
     bulk_phase_partitioning,
+    collision_kernel,
 )
 from .constants import convert_to, in_unit, si

PySDM/physics/collision_kernel/__init__.py

@@ -0,0 +1,13 @@
+"""
+Collision kernels physics
+"""
+
+from .neglect import Neglect
+from .constantK import ConstantK
+from .electric import Electric
+from .geometric import Geometric
+from .golovin import Golovin
+from .hydrodynamic import Hydrodynamic
+from .linear import Linear
+from .simple_geometric import SimpleGeometric
+from .long1974 import Long1974

PySDM/physics/collision_kernel/constantK.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.ConstantK`
+"""
+
+
+class ConstantK:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/electric.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Electric`
+"""
+
+
+class Electric:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/geometric.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Geometric`
+"""
+
+
+class Geometric:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/golovin.py

@@ -0,0 +1,36 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.golovin`
+"""
+
+import numpy as np
+from scipy import special
+
+
+class Golovin:
+    def __init__(self, _):
+        pass
+
+
+@staticmethod
+def analytic_solution(b, x, t, x_0, N_0):
+
+    def analytic_solution_helper(x, tau, x_0):
+        sqrt_tau = np.sqrt(tau)
+        result = float(
+            (1 - tau)
+            * 1
+            / (x * np.sqrt(tau))
+            * special.ive(1, 2 * x / x_0 * sqrt_tau)  # pylint: disable=no-member
+            * np.exp(-(1 + tau - 2 * sqrt_tau) * x / x_0)
+        )
+        return result
+
+    tau = 1 - np.exp(-N_0 * b * x_0 * t)
+
+    if isinstance(x, np.ndarray):
+        func = np.vectorize(lambda i: analytic_solution_helper(x[int(i)], tau, x_0))
+        result = np.fromfunction(func, x.shape, dtype=float)
+    else:
+        result = analytic_solution_helper(x, tau, x_0)
+
+    return result

PySDM/physics/collision_kernel/hydrodynamic.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Hydrodynamic`
+"""
+
+
+class Hydrodynamic:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/linear.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Linear`
+"""
+
+
+class Linear:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/long1974.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.Long1974`
+"""
+
+
+class Long1974:
+    def __init__(self, _):
+        pass

PySDM/physics/collision_kernel/neglect.py

@@ -0,0 +1,12 @@
+"""
+disable liquid liquid collisions
+"""
+
+
+class Neglect:
+    def __init__(self, _):
+        pass
+
+    @staticmethod
+    def collision_kernel(const, output, is_first_in_pair):
+        pass

PySDM/physics/collision_kernel/simple_geometric.py

@@ -0,0 +1,8 @@
+"""
+Formulae supporting `PySDM.dynamics.collisions.collision_kernels.SimpleGeometric`
+"""
+
+
+class SimpleGeometric:
+    def __init__(self, _):
+        pass

PySDM/physics/constants_defaults.py

@@ -370,6 +370,22 @@
 J_HOM = np.nan
 """ constant ice nucleation rates """
 
+GOLOVIN_b = 5e3 * si.s
+""" Golovin collision kernel constant """
+
+CONSTANTK_a = np.nan
+""" CONSTANTK collision kernel constant """
+
+LINEAR_a = np.nan
+LINEAR_b = np.nan
+""" LINEAR collision kernel constants """
+
+GEOMETRIC_collection_efficiency = 1
+""" GEOMETRIC collision kernel constants """
+
+SIMPLE_GEOMETRIC_c = 1
+""" SIMPLE_GEOMETRIC collision kernel constants """
+
 STRAUB_E_D1 = 0.04 * si.cm
 """ [Straub et al. 2010](https://doi.org/10.1175/2009JAS3175.1) """
 STRAUB_MU2 = 0.095 * si.cm