Add Spontaneous Breakup

PySDM/backends/numba/impl/_algorithmic_methods.py

@@ -173,6 +173,34 @@ def breakup(n, idx, length, attributes, gamma, n_fragment, healthy, is_first_in_
                                             attributes.data, gamma.data, n_fragment.data, healthy.data,
                                             is_first_in_pair.indicator.data)
 
+    ## TODO: Ben implement spontaneous breakup
+    @staticmethod
+    @numba.njit(**conf.JIT_FLAGS)
+    def spontaneous_breakup_body(n, length, attributes, gamma, rand, n_fragment, healthy):
+        for i in numba.prange(length):
+            gamma[i] = np.ceil(gamma[i] - rand[i])
+            # no spontaneous breakup occurs
+            if gamma[i] == 0:
+                continue
+            # breakup does occur
+            if n_fragment[i] > 0:
+                n[i] = n[i] * int(n_fragment[i])
+                # what is this?
+                # EJ: this changes the attributes and multiplicities when a coalescence occurs
+                for a in range(0, len(attributes)):
+                #     attributes[a, k] += gamma[i] * attributes[a, j]
+                     attributes[a, i] /= n_fragment[i]
+            # should not be an issue unless n_frag < 1
+            if n[i] == 0:
+                healthy[0] = 0
+
+    @staticmethod
+    def spontaneous_breakup(n, length, attributes, gamma, u01, n_fragment, healthy):
+        AlgorithmicMethods.spontaneous_breakup_body(n.data, length,
+                                            attributes.data, gamma.data, u01.data, 
+                                            n_fragment.data, healthy.data,
+                                            )
+
     # Emily: SLAMS fragmentation function
     @numba.njit(**{**conf.JIT_FLAGS})
     def slams_fragmentation_body(n_fragment, probs, rand):
@@ -270,6 +298,18 @@ def linear_collection_efficiency(params, output, radii, is_first_in_pair, unit):
         return AlgorithmicMethods.linear_collection_efficiency_body(
             params, output.data, radii.data, is_first_in_pair.indicator.data, radii.idx.data, len(is_first_in_pair), unit)
 
+    @staticmethod
+    @numba.njit(**conf.JIT_FLAGS)
+    def exponential_body(output, scaling, radii, idx, length, unit):
+        output[:] = 0
+        for i in numba.prange(length):
+            output[i] = np.exp(scaling * radii[idx[i]] / unit)
+
+    @staticmethod
+    def exponential(output, scaling, radii, unit):
+        return AlgorithmicMethods.exponential_body(
+            output.data, scaling, radii.data, radii.idx.data, len(radii), unit)
+
     @staticmethod
     @numba.njit(**conf.JIT_FLAGS)
     def interpolation_body(output, radius, factor, b, c):
@@ -328,7 +368,6 @@ def normalize(prob, cell_id, cell_idx, cell_start, norm_factor, dt, dv):
         return AlgorithmicMethods.normalize_body(
             prob.data, cell_id.data, cell_idx.data, cell_start.data, norm_factor.data, dt, dv)
 
-
     @staticmethod
     @numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}})
     def remove_zero_n_or_flagged(multiplicity, idx, length) -> int:

PySDM/dynamics/__init__.py

@@ -9,4 +9,5 @@
 from .eulerian_advection import EulerianAdvection
 from .ambient_thermodynamics import AmbientThermodynamics
 from .aqueous_chemistry import AqueousChemistry
-from .breakup import Breakup
+from .breakup import Breakup
+from .spontaneous_breakup import SpontaneousBreakup

PySDM/dynamics/spontaneous_breakup.py

@@ -0,0 +1,143 @@
+"""
+Single particle spontaneous breakup into uniformly sized fragments.
+
+Created at 10.05.21 by jb-mackay & edejong
+"""
+
+import numpy as np
+from PySDM.physics import si
+from PySDM.dynamics.impl.random_generator_optimizer import RandomGeneratorOptimizer
+from PySDM.dynamics.impl.random_generator_optimizer_nopair import RandomGeneratorOptimizerNoPair
+import warnings
+
+# checks that timestep is within this range
+default_dt_coal_range = (.1 * si.second, 100 * si.second)
+
+# TODO: define some a rate function
+class SpontaneousBreakup:
+
+    def __init__(self,
+                 breakup_rate, # rate function
+                 fragmentation, # frag function 
+                 seed=None,
+                 croupier=None,
+                 optimized_random=False,
+                 substeps: int = 1,
+                 adaptive: bool = False,
+                 dt_coal_range=default_dt_coal_range
+                 ):
+        assert substeps == 1 or adaptive is False
+
+        self.core = None # .core is the parent for all attributes, dynamics, etc
+        self.enable = True # on/off for the dynamic
+
+        self.breakup_rate = breakup_rate
+        self.fragmentation = fragmentation
+        self.seed = seed
+
+        assert dt_coal_range[0] > 0
+        self.croupier = croupier
+        self.optimized_random = optimized_random
+        self.__substeps = substeps
+        self.adaptive = adaptive
+        self.stats_n_substep = None
+        self.stats_dt_min = None
+        self.dt_coal_range = tuple(dt_coal_range)
+
+        # temporary storage for outputs (may need add some (and remove))
+        self.n_fragment = None
+        self.breakup_rate_temp = None
+        self.norm_factor_temp = None # renormalizes the system
+        self.prob = None
+        self.dt_left = None
+
+        # simulation diagnositcs (TODO)
+
+
+    # registering the process with the builder (initialization)
+    def register(self, builder):
+        self.core = builder.core
+
+        # outputs a random number (will need one for IF breakup? and for fragmentation number)
+        self.rnd_opt = RandomGeneratorOptimizerNoPair(optimized_random=self.optimized_random,
+                                                dt_min=self.dt_coal_range[0],
+                                                seed=self.seed) #self.core.formulae.seed+1)
+        self.rnd_opt_frag = RandomGeneratorOptimizerNoPair(optimized_random=self.optimized_random,
+                                                dt_min=self.dt_coal_range[0],
+                                                seed=self.seed) #self.core.formulae.seed+1)
+
+        self.optimised_random = None
+
+        if self.dt_coal_range[1] > self.core.dt:
+            self.dt_coal_range = (self.dt_coal_range[0], self.core.dt)
+        assert self.dt_coal_range[0] <= self.dt_coal_range[1]
+
+        # init temporary storage; swap to Storage (see Condensation)
+        self.n_fragment = self.core.Storage.empty(self.core.n_sd, dtype=float)
+        self.breakup_rate_temp = self.core.Storage.empty(self.core.n_sd, dtype=float)
+        self.norm_factor_temp = self.core.Storage.empty(self.core.mesh.n_cell, dtype=float)
+        self.prob = self.core.Storage.empty(self.core.n_sd, dtype=float) # does the event happen?
+        self.dt_left = self.core.Storage.empty(self.core.mesh.n_cell, dtype=float)
+
+        self.stats_n_substep = self.core.Storage.empty(self.core.mesh.n_cell, dtype=int)
+        self.stats_n_substep[:] = 0 if self.adaptive else self.__substeps
+        self.stats_dt_min = self.core.Storage.empty(self.core.mesh.n_cell, dtype=float)
+        self.stats_dt_min[:] = np.nan
+
+        # registering helper processes
+        self.rnd_opt.register(builder)
+        self.rnd_opt_frag.register(builder)
+        self.breakup_rate.register(builder)
+        self.fragmentation.register(builder)
+
+        if self.croupier is None:
+            self.croupier = self.core.backend.default_croupier
+
+        # Diagnostics
+        # self.collision_rate = self.core.Storage.from_ndarray(np.zeros(self.core.mesh.n_cell, dtype=int))
+        # self.collision_rate_deficit = self.core.Storage.from_ndarray(np.zeros(self.core.mesh.n_cell, dtype=int))
+
+    def __call__(self):
+        if self.enable:
+            if not self.adaptive:
+                for _ in range(self.__substeps):
+                    self.step()
+            else:
+                self.dt_left[:] = self.core.dt
+
+                # work on all particles
+                while self.core.particles.get_working_length() != 0:
+                    self.core.particles.cell_idx.sort_by_key(self.dt_left)
+                    self.step() # do breakup
+
+                self.core.particles.reset_working_length()
+                self.core.particles.reset_cell_idx()
+            self.rnd_opt.reset()
+            self.rnd_opt_frag.reset()
+
+    def step(self):
+        # (1) Make the random numbers for fragmentation
+        rand = self.rnd_opt.get_random_arrays() # does event occur
+        rand_frag = self.rnd_opt_frag.get_random_arrays() # for # fragments
+
+        # (2) Compute the probability of spontaneous breakup
+        self.compute_probability(self.prob)
+
+        # (3) Compute the number of fragments
+        self.compute_n_fragment(self.n_fragment, rand_frag)
+
+        # (4) Perform the collisional-breakup step: 
+        self.core.particles.spontaneous_breakup(gamma=self.prob, u01=rand, n_fragment=self.n_fragment)
+        if self.adaptive:
+            self.core.particles.cut_working_length(self.core.particles.adaptive_sdm_end(self.dt_left))
+
+    # (2) Does spont breakup occur?
+    def compute_probability(self, prob):
+        self.breakup_rate(self.breakup_rate_temp)
+        prob[:] = self.breakup_rate_temp
+        # breakup_rate * dt = expected number of breakups in time dt
+        prob[:] *= self.core.dt
+
+    # (3) Compute n_fragment
+    def compute_n_fragment(self, n_fragment, u01):
+        self.fragmentation(n_fragment, u01)

PySDM/physics/__init__.py

@@ -18,3 +18,5 @@
 import PySDM.physics.hydrostatics
 import PySDM.physics.breakup_fragmentations
 import PySDM.physics.coalescence_efficiencies
+import PySDM.physics.spontaneous_breakup_fragmentations
+import PySDM.physics.spontaneous_breakup_rates

PySDM/physics/spontaneous_breakup_fragmentations/__init__.py

@@ -0,0 +1,5 @@
+"""
+Single droplet breakup fragmentation functions.
+"""
+from .always_n import SpontaneousAlwaysN
+from .expo_cube import SpontaneousExpoCube

PySDM/physics/spontaneous_breakup_fragmentations/always_n.py

@@ -0,0 +1,24 @@
+import numpy as np
+
+
+class SpontaneousAlwaysN:
+    """Breakup into `n` droplets"""
+
+    def __init__(self, n):
+        self.core = None
+        self.N = n
+        self.N_vec = None
+        self.zeros = None
+
+
+    def __call__(self, output, u01):
+        output *= self.zeros
+        output += self.N_vec
+
+    def register(self, builder):
+        self.core = builder.core
+        N_vec_tmp = np.tile([self.N], self.core.n_sd)
+        zeros_tmp = np.tile([0], self.core.n_sd)
+        self.N_vec = self.core.PairwiseStorage.from_ndarray(N_vec_tmp)
+        self.zeros = self.core.PairwiseStorage.from_ndarray(zeros_tmp)
+

PySDM/physics/spontaneous_breakup_fragmentations/expo_cube.py

@@ -0,0 +1,29 @@
+'''
+Linear fragmentation function in volume.
+
+For a fragmentation function of the form:
+    Q(D0, D) dD = A * (D0)**3 * exp(-B * D) dD
+where D0 is the parent droplet size, Q(D) dD is the number of drops of size D
+formed from the breakup. Then it can be shown that:
+    N_frag = int (Q(D0, D) dD) = A * (D0)**3 / B
+and furthermore for consistency with mass conservation:
+    D0**3 = int(D**3 Q(D0, D) dD) => 6*A = B**4
+
+Therefore this function takes a single argument C = A/B = 1/6 * B**3, and 
+translates to a fragmentation function which is cubic in the parent size D0.
+    N_frag = C * (D0)**3
+'''
+class SpontaneousExpoCube:
+
+    def __init__(self, C):
+        self.core = None
+        self.C = C  
+
+    def __call__(self, output, u01):
+        output[:] = self.core.particles['volume']
+        output *= self.C
+
+    def register(self, builder):
+        self.core = builder.core
+        builder.request_attribute('volume')
+

PySDM/physics/spontaneous_breakup_rates/__init__.py

@@ -0,0 +1,5 @@
+"""
+Single droplet spontaneous breakup rates.
+"""
+from .constant import Constant
+from .exponential import Expon

PySDM/physics/spontaneous_breakup_rates/constant.py

@@ -0,0 +1,12 @@
+class Constant:
+    '''Constant rate of breakup for particles of any size'''
+    def __init__(self, a):
+        self.a = a
+        self.core = None
+
+    def __call__(self, output):
+        output *= 0
+        output += self.a
+
+    def register(self, builder):
+        self.core = builder.core

PySDM/physics/spontaneous_breakup_rates/exponential.py

@@ -0,0 +1,16 @@
+from PySDM.physics import constants as const
+
+class Expon:
+    '''Exponential breakup rate by particle diameter: a * r**b'''
+    def __init__(self, a, b):
+        self.a = a
+        self.b = b
+        self.core = None
+
+    def __call__(self, output):
+        self.core.backend.exponential(output, self.b, self.core.particles['radius'], const.si.um)
+        output *= self.a
+
+    def register(self, builder):
+        self.core = builder.core
+        builder.request_attribute('radius')

PySDM/state/particles.py

@@ -164,6 +164,23 @@ def breakup(self, gamma, n_fragment, is_first_in_pair):
             if isinstance(attr, ExtensiveAttribute):
                 attr.mark_updated()
 
+    ## TODO Ben: def spontaneous_breakup(self, gamma=prob, u01, n_fragment) gamma = prob of breakup occuring
+    def spontaneous_breakup(self, gamma, u01, n_fragment):
+        self.core.bck.spontaneous_breakup(n=self['n'],
+                              length=self.SD_num,
+                              attributes=self.get_extensive_attrs(),
+                              gamma=gamma,
+                              u01=u01,
+                              n_fragment=n_fragment,
+                              healthy=self.__healthy_memory,
+                             )
+        self.healthy = bool(self.__healthy_memory)
+        self.core.particles.sanitize()
+        self.attributes['n'].mark_updated()
+        for attr in self.attributes.values():
+            if isinstance(attr, ExtensiveAttribute):
+                attr.mark_updated()
+
     def adaptive_sdm_end(self, dt_left):
         return self.core.bck.adaptive_sdm_end(dt_left, self.core.particles.cell_start)