diff --git a/PySDM/backends/numba/impl/_algorithmic_methods.py b/PySDM/backends/numba/impl/_algorithmic_methods.py index f6cbdf583f..9c53c6f875 100644 --- a/PySDM/backends/numba/impl/_algorithmic_methods.py +++ b/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: diff --git a/PySDM/dynamics/__init__.py b/PySDM/dynamics/__init__.py index c4652229e0..956212e77f 100644 --- a/PySDM/dynamics/__init__.py +++ b/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 \ No newline at end of file +from .breakup import Breakup +from .spontaneous_breakup import SpontaneousBreakup \ No newline at end of file diff --git a/PySDM/dynamics/spontaneous_breakup.py b/PySDM/dynamics/spontaneous_breakup.py new file mode 100644 index 0000000000..1115ba8e08 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/PySDM/physics/__init__.py b/PySDM/physics/__init__.py index 078bd7eccb..69b351cb97 100644 --- a/PySDM/physics/__init__.py +++ b/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 \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_fragmentations/__init__.py b/PySDM/physics/spontaneous_breakup_fragmentations/__init__.py new file mode 100644 index 0000000000..47ab773f12 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_fragmentations/always_n.py b/PySDM/physics/spontaneous_breakup_fragmentations/always_n.py new file mode 100644 index 0000000000..7808a31f2d --- /dev/null +++ b/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) + \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_fragmentations/expo_cube.py b/PySDM/physics/spontaneous_breakup_fragmentations/expo_cube.py new file mode 100644 index 0000000000..7b906c616c --- /dev/null +++ b/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') + \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_rates/__init__.py b/PySDM/physics/spontaneous_breakup_rates/__init__.py new file mode 100644 index 0000000000..7581bd755d --- /dev/null +++ b/PySDM/physics/spontaneous_breakup_rates/__init__.py @@ -0,0 +1,5 @@ +""" +Single droplet spontaneous breakup rates. +""" +from .constant import Constant +from .exponential import Expon \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_rates/constant.py b/PySDM/physics/spontaneous_breakup_rates/constant.py new file mode 100644 index 0000000000..95cf959dc1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/PySDM/physics/spontaneous_breakup_rates/exponential.py b/PySDM/physics/spontaneous_breakup_rates/exponential.py new file mode 100644 index 0000000000..300abc91ec --- /dev/null +++ b/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') \ No newline at end of file diff --git a/PySDM/state/particles.py b/PySDM/state/particles.py index 911c054889..0439b4e8cf 100644 --- a/PySDM/state/particles.py +++ b/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) diff --git a/PySDM_tests/breakup_tests/spontaneous_breakup_test.ipynb b/PySDM_tests/breakup_tests/spontaneous_breakup_test.ipynb new file mode 100644 index 0000000000..1dbde9d58a --- /dev/null +++ b/PySDM_tests/breakup_tests/spontaneous_breakup_test.ipynb @@ -0,0 +1,515 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 2, + "source": [ + "import sys\n", + "import numpy as np\n", + "\n", + "from PySDM.backends import CPU, GPU\n", + "from PySDM.builder import Builder\n", + "from PySDM.environments import Box\n", + "from PySDM.dynamics import Coalescence\n", + "from PySDM.dynamics import SpontaneousBreakup\n", + "from PySDM.initialisation.spectral_sampling import ConstantMultiplicity\n", + "\n", + "from PySDM.products.state import ParticlesVolumeSpectrum, ParticlesConcentration, ParticleMeanRadius\n", + "from PySDM.products.stats.timers import WallTime\n", + "\n", + "from matplotlib import pyplot" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## AlwaysN fragments, Golovin kernel" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 81, + "source": [ + "from PySDM.initialisation.spectra import Exponential\n", + "from PySDM.physics.spontaneous_breakup_rates import Constant\n", + "from PySDM.physics.spontaneous_breakup_fragmentations import SpontaneousAlwaysN\n", + "from PySDM.physics.constants import si\n", + "from PySDM.physics.formulae import Formulae\n", + "#from pystrict import strict\n", + "\n", + "#@strict\n", + "class Settings:\n", + "\n", + " def __init__(self):\n", + " self.formulae = Formulae()\n", + " self.n_sd = 2**8\n", + " self.n_part = 100 / si.cm**3\n", + " self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres)\n", + " self.dv = 1 * si.m**3\n", + " self.norm_factor = self.n_part * self.dv\n", + " self.rho = 1000 * si.kilogram / si.metre**3\n", + " self.dt = 1 * si.seconds\n", + " self.adaptive = False\n", + " self.seed = 44\n", + " self._steps = [0, 1, 2]\n", + " self.breakup_rate = Constant(1 / si.second)\n", + " self.fragmentation = SpontaneousAlwaysN(n=3)\n", + " self.spectrum = Exponential(norm_factor=self.norm_factor, scale=self.X0)\n", + " self.radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(100 * si.um), num=128, endpoint=True)\n", + " self.radius_range = [0 * si.um, 1e6 * si.um]\n", + "\n", + " @property\n", + " def output_steps(self):\n", + " return [int(step/self.dt) for step in self._steps]" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 85, + "source": [ + "settings = Settings()\n", + "backend = CPU\n", + "\n", + "builder = Builder(n_sd=settings.n_sd, backend=backend, formulae=settings.formulae)\n", + "builder.set_environment(Box(dv=settings.dv, dt=settings.dt))\n", + "attributes = {}\n", + "attributes['volume'], attributes['n'] = ConstantMultiplicity(settings.spectrum).sample(settings.n_sd)\n", + "breakup = SpontaneousBreakup(settings.breakup_rate, settings.fragmentation, adaptive=settings.adaptive)\n", + "builder.add_dynamic(breakup)\n", + "\n", + "products = [ParticlesVolumeSpectrum(), WallTime(), ParticleMeanRadius(), ParticlesConcentration(radius_range = settings.radius_range)]\n", + "core = builder.build(attributes, products)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 86, + "source": [ + "for step in settings.output_steps:\n", + " #print(step, flush=True)\n", + " core.run(step - core.n_steps)\n", + " pyplot.step(x=settings.radius_bins_edges[:-1] / si.micrometres, \n", + " y=core.products['dv/dlnr'].get(settings.radius_bins_edges) * settings.rho,\n", + " where='post', label=\"t = {step*settings.dt}s\")\n", + " print(core.products['radius_m1'].get(), core.products['n_a_cm3'].get())\n", + " \n", + "pyplot.xscale(\"log\")\n", + "pyplot.xlabel(\"radius (um)\")\n", + "pyplot.ylabel(\"dm/dlnr\")\n", + "pyplot.legend()" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[27.26101585] [99.997952]\n", + "[18.95139778] [296.087686]\n", + "[13.13134485] [884.356888]\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 86 + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsr0lEQVR4nO3deZhU1bnv8e8LNqDg2IoohMEJ0wy22sKJ05WoRPEoxmBA8YitRj3EjE/yiCfOV07gmEdv4qxROEa5tpp4QUFMjPGiRhFUDIOSCwgCgcggyCDQdL/3j6pd7C6qu6u6a9rVv8/z8HTVrl27Vonst9da73qXuTsiIiLpalfoBoiISLQocIiISEYUOEREJCMKHCIikhEFDhERyYgCh4iIZGSfQjcgHw499FDv3bt3oZshIhIp77///np3Pyz5eJsIHL1792bu3LmFboaISKSY2YpUxzVUJSIiGVHgEBGRjChwiIhIRtrEHIdI1NXW1rJq1Sp27NhR6KZICerUqRM9evSgrKwsrfMVOEQiYNWqVey///707t0bMyt0c6SEuDsbNmxg1apV9OnTJ633aKhKJAJ27NhBeXm5goZknZlRXl6eUW9WPQ6Rlpg7Cea/sOf5gBFQVZ3Tj1TQkFzJ9P8t9TgkEp7/+/NUz6ymemY1z//9+UI3JxY01s6PPV47v2EQKUGbNm3ioYceyuo1Z8yYwfHHH88vfvGLJs/7z//8z1Z/1uTJk1m+fDkt2X9o3rx5zJgxI/H8jTfe4K9//Wvi+SOPPMJTTz3V6jamq3fv3qxfvz7l38nIkSMZMGAACxcuzGkbFDgkEmYsm8HijYtZvHExM5bNaP4N+dBtAFRPj/2MoA1bd7J03VaWrtvKhq07mzw3F4Hj8ccf595772X8+PFNnteawLF69WquvfZaVq5cyVtvvcUNN9yQ8TWaCxw33HADV155ZYvb2FKp/k5qamoYPXo0zz77bE4/W4FDIqPvIX3pe0jfQjejZGz6qpYdtXXsqK1j01e1TZ47btw4li5dSmVlJT//+c+z8/mbNtG1a9fE8zVr1nDmmWdSWVlJ//79efPNNxk3bhxfffUVlZWVjB49GoCnn36aQYMGUVlZyfXXX09dXR0AXbp04Sc/+Qn9+vXj7LPPZt26dXTv3p3x48fzxBNP8Oyzz/Lwww8D8Jvf/IaKigoGDhzIqFGjANi2bRtXX301gwYN4sQTT2Tq1Kns2rWL2267jZqaGiorK5k4cSKPPPII9913H5WVlbz55pvccccd/OpXvwLgrLPO4qabbmLQoEEcd9xxvPnmmwBs376d7373u1RUVPDtb3+bwYMHM3fuXOrq6rjqqqvo378/AwYM4L777tvrv9OGDRsYOnQo/fr149prr030mhr7O+nWrRubNm1qcI1U37c1NMchEjF3vrSQRf/4stXX2VFbl3h8dNcu3PvdykbPnTBhAgsWLGDevHkpXz/jjDPYsmXLXsd/9atfcc4556R8T11dHe3a7fnddcqUKXzrW9/iF7/4BXV1dWzfvp0zzjiDBx54IPG5H3/8MTU1Nbz99tuUlZUxduxYnnnmGa688kq2bdtGVVUV9913H3fddRd33nkn//Ef/8Htt9/O1VdfTZ8+ffj+97/Pww8/zIQJE/j000/p2LFj4iY7fvx4vvnNb/Lkk0+yadMmBg0axDnnnMNdd93F3LlzeeCBBwD46quv6NKlCz/72c8A+POf/9zge+3evZv33nuPGTNmcOedd/Laa6/x0EMPcfDBB7No0SIWLFhAZWXsv/W8efNYvXo1CxYsANjrhg9w5513cvrpp3Pbbbcxffp0nnjiiSb/Ttq1a5cIpoFU37c1FDhEpNWC36zTtW3bNpYuXcqRRx6ZOHbKKadw9dVXU1tby8UXX5y4uYb9+c9/5v333+eUU04BYjfxoNfSrl07Ro4cCcAVV1zBJZdcwpFHHsnjjz/O5MmTOeOMM7jiiisAGDhwIKNHj+biiy/m4osvBuCPf/wj06ZNS/QeduzYwWeffZbR9wK45JJLADj55JNZvnw5AG+99RY/+tGPAOjfvz8DBw4E4KijjmLZsmX84Ac/4IILLmDo0KF7XW/WrFn84Q9/AOCCCy7g4IMPbvLzu3fvzsKFC6mvr08E5lTftzUUOEQi5vYL+2XlOkvXbc3KdSCzHsfs2bMZNmwYQ4cOpVu3bonjZ555JrNmzWL69OlcddVV/PSnP91r7sDdGTNmDL/85S+bbVM4U+iqq65q8Nr06dOZNWsWL730EuPHj2f+/Pm4O7///e/p27fhcOjs2bOb/aywjh07AtC+fXt2797d5LkHH3wwH330Ea+++iqPPPIIzz33HE8++WRGn5fszDPPZMuWLfTs2ZOPPvqI8vLylN93n31afvvXHIeINGv//fdPGRgCb775JvPmzdvrT6phqsGDB/OPf/yDd999l6VLlyaOr1ixgsMPP5zvfe97XHvttXzwwQcAlJWVUVsbm4M5++yzeeGFF/j8888B2LhxIytWxAq41tfX88ILsey2KVOmcPrpp6dsa319PStXrmTIkCFMnDiRzZs3s3XrVr71rW9x//33J+YQPvzww5Tfvbn/FqmcdtppPPfccwAsWrSI+fNjGXnr16+nvr6e73znO9x9992J7xx25plnMmXKFABeeeUVvvjiiybbMWPGDA4//HBWrlxJeXl5o9+3NdTjkLx5/u/PN8iIGnbUMC497tICtkjSVV5ezmmnnUb//v05//zzueeee1p1vY4dO9KrV6/ETRBi2Ur33HMPZWVldOnSJZHiet111zFw4EBOOukknnnmGe6++26GDh1KfX09ZWVlPPjgg/Tq1YvOnTvz3nvvcffdd9O1a1dqampSfnZdXR1XXHEFmzdvxt354Q9/yEEHHcStt97Kj3/8YwYOHEh9fT19+vTh5ZdfZsiQIUyYMIHKykpuvvlmLrzwQkaMGMHUqVO5//770/q+Y8eOZcyYMVRUVHD88cfTr18/DjzwQFavXk11dTX19fUAiZ7UI488AsQytm6//XYuu+wy+vXrx6mnnkrPnj2Bxv9OvvjiC44++uhEj6ux79sa1pK85qipqqpy7cdReNUzq1m8cTF9D+mb+DnpvElpvzcs3fflzKQLYj+rpzd8nCMff/wxX//617N6zeShqqMP65LV6zdn2LBh3HjjjQwbNiwr1+vSpUurf5POlbq6Ompra+nUqRNLly7lnHPOYfHixXTo0CHrn/Vf//VfbNiwgYkTJ2b0vlT/j5nZ++5elXyuehySV0GwSA4E0vbccMMN3HLLLbz99tvNruWIuu3btzNkyBBqa2txdx566KGcBI1Ro0axYsWKRNpxrihwiEhBXHTRRVx00UVZu16x9jYgNh+Rj1GPXC/8C2hyXEREMqLAISIiGVHgEBGRjGiOQyRbwqXW81BmPR82bN3ZoI7VQfuWUd6lYwFbJMVAPQ6RbAlKrZdQmfWgEOKXmzcx+fFHmy2GmAmVVW8ZlVWXSAvvkVE0+2QUWrcBkS2z3phOZe05eJ/d1PzuiaxeV2XVW0dl1SWSgj0ygOLaJ0Oybty4cXy2/FMuHHKqyqqrrLrmOKR1Ir+gL4rzEq+M27P7YCscGSqrvvPQCrjk3kbPnTBhAh989Dde+stfU64wV1n1PVRWXaTUhbeAhdSBI9X+4tKAyqrvobLqIm1Bc3MSQXDpNiArv+m32vkTEg/DWU+ZZjz9I6lW1UGtaJLKqu+hsuoiElOk+4unu/1rJvuLp7L//vuzrYmSHiqr3rRSK6uuwCGSC2vnx6rmTrogNtSVQ53K2tOprH2T52Syv3gq5eXlnDzoXzj/zEFZmRxvrKz6CSecwIknnkhNTU1iaCcoqz569GgqKioSZdUHDhzIueeey5o1awASZdX79+/P66+/zm233Zbys4My4wMGDODEE09sUFa9traWgQMH0q9fP2699VYAhgwZwqJFi6isrKSmpoYLL7yQF198MTE5no6xY8eybt06KioquOWWWxqUVT/rrLOorKzkiiuuaFBWPSitfvvttzNr1iz69evHH/7wh5Rl1cN/J42VVU/+vq2hoSqRbAvPgQRDW0Uw6d5ccGnOfY/EhlCyVX59v/32S/QcAMaMGcOYMWP2Om/ixIkNSoSPHDkyMZeR7N57G5/gD5SVlfHWW2/tdXzffffl0Ucf3ev4IYccwpw5cxoc+9vf/pZ4fMYZZyQev/HGG4nHhx56aGKOo1OnTjz99NMNyqr36tWLDh06pOxlhNOGy8vL+eMf/5jyuwQ9kbDPP/+czp07J5439n1bQ4FDJNuqqvcEimCvDtmLyqqrrHpKZnYe8GugPfBbd5+Q9HpH4CngZGADMNLdl5vZucAEoAOwC/i5u78ef8/JwGRgX2AG8CNvC7tRiZQYlVXPvnyVVc9Z4DCz9sCDwLnAKmCOmU1z90Wh064BvnD3Y8xsFDARGAmsBy5093+YWX/gVaB7/D0PA98DZhMLHOcBr+Tqe0jLhbeKHXZUdnZ5E5HCy2WPYxCwxN2XAZjZs8BwIBw4hgN3xB+/ADxgZubuH4bOWQjsG++dHAIc4O7vxq/5FHAxChxFKbyyXERKRy4DR3dgZej5KmBwY+e4+24z2wyUE+txBL4DfODuO82se/w64Wt2R4pW30P6Nn+SiERKUU+Om1k/YsNXey+nbP691wHXAYn0NSluyUNblx53aWEbFKwYDxb/iQiQ23Ucq4GvhZ73iB9LeY6Z7QMcSGySHDPrAbwIXOnuS0Pn92jmmgC4+2PuXuXuVYcddlgrv4rkQzC0VTQFE8NBo42XGdm0aRNPP/l4Vq+psuotU+pl1ecAx5pZHzPrAIwCpiWdMw0IErdHAK+7u5vZQcB0YJy7vx2c7O5rgC/N7F8strrlSmBqDr+D5FnfQ/oW1/BWsGK8CNZhZCq8WnxHbV3zb2jCpk2beGZydgOHyqq3TkmWVXf33cCNxDKiPgaec/eFZnaXmQU5eE8A5Wa2BPgpMC5+/EbgGOA2M5sX/xPUXx4L/BZYAixFE+MiKQWrxSG2+O+gfcuaPD8caJJLk6isusqqh+V0jsPdZxBLmQ0fuy30eAew10C2u98N3N3INecC/bPbUpHomPjeRD7Z+AkAXyX1JPYNrQ4PXmvsWPi9fQ44lmsrfsKO2jo6lbWPlSYJXVdl1VVWPayoJ8dFMlV0E+wR06msPUcf1oWl6zJbTKey6nuorLpISDoL+sLnLN64OO/zFclrR0oxcNw06KbE4+QbfLg3ELzW2LFMg0NTVFZ9D5VVFwlJJ+spfOPue0jfgqwYb3SCfe6kvFWsLTUqq66y6mHqcUhG0ulBBNvJBooitRbS2+1PUgqXVR/+rxdwzz33tOp6jZVVv+eeeygrK6NLly6JFNegrPpJJ53EM888kyirXl9fT1lZGQ8++CC9evVKlFW/++676dq1KzU1NSk/OygzvnnzZty9QVn1H//4xwwcOJD6+nr69OnDyy+/zJAhQ5gwYQKVlZXcfPPNXHjhhYwYMYKpU6dy//33p/V9x44dy5gxY6ioqOD4449vUFa9urqa+vp6gAZl1SGWsXX77bdz2WWX0a9fP0499dSUZdXPP//8xN9JY2XVk79vayhwSEHlfWgrnYV8RbrwL7zbH9BsllS2qaz6Hm29rLqGqqSgimFoay9FuvAvnF7b0g2ZiklQVr25BYClYPv27Zx++umccMIJfPvb385pWfUXX3yRyy67LOvXDlOPQwquKIe2goV/qQTDXQXojaST9RT0TILU2mKlsurZF/my6iIlKdwDSbXTXxEMbYWDRrrDWeGFgiLNUeAQyUR4d79AY8Eky9y9QYppU4KeSToO2rcssdgv3/MmUhwyreGlwCHSWqmCSZZ16tSJbVs202n/A7N+7fIuHSnv0jHr15VocHc2bNhAp06d0n6PAoc0KZOsp/Akt2RXjx49mPr2fA7bbz21GzqybsvOlOftWr/ntV3r9w4G4dfSOS/d60q0derUiR49ejR/YpwChzQpyHoKFtUNO2pYysnrcDZUUWRGBYLUWii69NpMlJWV8fzi2I275vpK7nj0nZTnhV+rub5yr9fDr6VzXrrXlbZFgUOalU7W06XHXZpxeY+89FDCqbVBem0QSJKtnR9bVR7hAFMoU2Z/xtR5sa1xhld25/LB2jytlClwSEHktYeSnFqbKnCEJ7WLbP1GcFNetOZLKo44oNDNSSloX0CBo7QpcEhBtKSHklN5mOBuqXDQGF7ZvdHzFq35kpGPvlOwAFOsQU2yT4FDJAIqjjiAmuu/0ejr4YDSXIARaS0FDpEScPngnhoekrxR4JDIivymTeGML4jNqxTpcJlImIocSkrP//15qmdWN9gUqdiksz9IUQuXeV87v/FsL5Eio8AhKYXXbxTVuowkjW7aFBVBxpfSfyVCNFQljUpevxEpW9bCtnWw9nPdlEWyTIFDStO2dbBrW9GtyUhXeEFdMa/fkLZJgUNKV4fOcFUje2oUufDaDaXXSrFR4BApUs2t3RApFE2Oi7Qx4RXmIi2hHodIG6IV5pINChySoL038iRY+FeAKrzZWGGuSriiwCEJqfbeSKVo996APWm4u7bFJsezIBxQIQur1MNBI98ZX1lYrV7ISrjhoAUKXIWiwCENpLN2o2j33oCGQaPzYY2fl8HeG+GAGnyPVpc3SS71ni/hoBWsWm9BmZNCpQeHs82C4KXAkX+aHJecG3bUsGZ7MVnVoXPsxrh/t9SvDxixJ1ik+Vt/EFBzHfimzP4s7YnrRWu+bNkEd8RXqwfZZlrbUjjqcUjOae+N9KW790b4NU1wS74pcIgUmXTWb6iMuhSShqpERCQj6nFIImuouRRcaYVwNlMB0nCzJZhT0fxC25bTHoeZnWdmi81siZmNS/F6RzOrib8+28x6x4+Xm9lfzGyrmT2Q9J434tecF//TNZffoS2ISgn1Rm1ZuydLate2QrcmtfDeGxEtvDi8srtqZwmQwx6HmbUHHgTOBVYBc8xsmrsvCp12DfCFux9jZqOAicBIYAdwK9A//ifZaHefm6u2t0WRLqEepOCW0XwabiEVKgU3S5LnVcLrKaRtyWWPYxCwxN2Xufsu4FlgeNI5w4H/jj9+ATjbzMzdt7n7W8QCiEjzOnTek2LaWBpuEfvnlh0sXLNZ9aMkEnIZOLoDK0PPV8WPpTzH3XcDm4HyNK49KT5MdauZWaoTzOw6M5trZnPXrVuXeetF8mj91p1s31WnYSCJhChOjo9299Vmtj/we+DfgKeST3L3x4DHAKqqqjy/TRTJ3H4d2quMukRCLnscq4GvhZ73iB9LeY6Z7QMcCGxo6qLuvjr+cwswhdiQmIiI5EkuexxzgGPNrA+xADEKuDzpnGnAGOAdYATwurs32juIB5eD3H29mZUB/wq8lovGl7pMKuEWUkmkCoezqURKQM4Ch7vvNrMbgVeB9sCT7r7QzO4C5rr7NOAJ4HdmtgTYSCy4AGBmy4EDgA5mdjEwFFgBvBoPGu2JBY3Hc/UdSlm6lXALLTlVOFylNheVcLMunHYbwRRckVRyOsfh7jOAGUnHbgs93gGkLGLk7r0buezJ2WpfWxeVFNxwOxsEjnQr4RZSE3WxwiXCf7arjv06tM9ny0RaTCVHJNqaq4RbxML7WuzXoT2HdulY4BaJpCeKWVUiJSNR0HDSgYVuikja1OOQorR442KqZ1YnNk4SkeKhHkcbE4UspfBEfTFP3Iu0Vc0GjvjK7B7uvrK5c6X4RaGgYdFt/CQiDTQbONzdzWwGoCT0EhGVbKqUghTcSRdAbRGn4ZYAlVCXxqQ7VPWBmZ3i7nNy2hqR5kSlEm4TgjTcYIvYYqStaaUp6QaOwcBoM1sBbAOMWGdkYM5aJtKYDp3hqukwszj3DW9OuvuKF5K2ppWmpBs4vpXTVogUWLgEC8Qm6HM5z5LOvuIZCe8wCLFV6o0sPCxW4QWRwyu7K3AVsbTScd19BbGy6LWAh/5IRDz/9+eV3tqEIGkAYqnADVaoR0F4h8G18xsGkYgIemKL1nypTaKKXFqBw8x+APwT+BMwPf7n5Ry2S7IsCtlUhRYkDRRrmnKzgh0GI1xMMdiaVopbukNVPwL6unuTJc+luEU6myoYilEmlUjBpRs4VhLbnU+kMIKhmCO6RjqTCijqbCqRdKQbOJYBb5jZdGBncNDd781Jq0RS6TYAunUtdCtaJJxJVczZVCLpSDdwfBb/0yH+R0QylPVMKpECSStwuPuduW6IZF9UdvkTkWhpMnCY2Us0kXbr7hdlvUWSNVHZ5a9JwaT42vmRzhYSKSXN9Th+lZdWSM5EOpMKGgaNASNg/axCtygjUSgvIpKpJgOHu//ffDVEpFHB+gSAmdEKHFEoLyKSqeaGqubT9FCValVJ5KQqL5JLOZkUD5cYGTAiu9fOA5UXibbmhqr+Nf7z+/Gfv4v/vAKVHJGICs/9RLYES7jESASF91sHFDgiprmhqhUAZnauu58YeukmM/sAGJfLxknLRGGXv0IL5n6qI1phF4h8soDmfKIr3XUcZmanufvb8Senov3Ki1bk61KFh2GUTSVSdNINHNcAT5rZgfHnm4Crc9IiyYpIZ1OFM6mCbKoIUXkRKXXNTY5/A3jX3d8HTggCh7urbpXkVjiTKmJUXkRKXXM9jiuBB83s78BMYKa7r819s0SiTeVFpJQ1Nzn+7wBmdjxwPjA53uv4C7FA8ra71+W8lSIiUjTSrVX1CfAJcJ+Z7QsMAS4F7gWqctc8kdbJ95qNnIp40oDmfkpH2plRZnawmQ0Evg6sBSa5u4KGFLV8bgk7ZfZnjHz0nQbrE7IqvHYjgkkD4bUbmvuJtrR6HGb2P4GriO3LUR8/7MA3c9MskezJ15qNvJQXSU4aiNje4slzP9pbPJrSTcf9LnC0u+/KZWNEok6T4tIWpBs4FgAHAZ/nrinSEuEx/GFHDePS4y4tcItEpNSlGzh+CXxoZgtouHWs9uMosPAYPhC9wBGe8IWsjtsH/10iX3Yl4nuSqLR86Uk3cPw3MBGYz545jmaZ2XnAr4H2wG/dfULS6x2Bp4CTgQ3ASHdfbmblwAvAKcBkd78x9J6TgcnAvsAM4Efu3qYLLkb6xhi+IWaxaF84eyrXmVQ5vzEm70kSMSotX3rSDRzb3f03mVzYzNoDDwLnAquAOWY2zd0XhU67BvjC3Y8xs1HEgtNIYAdwK9A//ifsYeB7wGxigeM84JVM2iZFJpjwnXRB1i556XGX5q33VZBJ8YjR3E9pSTdwvGlmvwSm0XCo6oMm3jMIWOLuywDM7FlgOBAOHMOBO+KPXwAeMDNz923AW2Z2TPiCZnYEcIC7vxt//hRwMQocUmC6MUpbkm7gCEqqD47/NJpPx+0OrAw9XxV6/17nuPtuM9sMlAPrm7jmqqRrpvwVz8yuA64D6NlTtf5FRLKluSKHP40/fJlYoLDQy0U9r+DujwGPAVRVVRV1WyV6cr4KOuKrxEGT4qWsuZXj+8f/nAz8O3AEcCRwPXBSM+9dDXwt9LxH/FjKc8xsH+BAYpPkTV2zRzPXFMm5nK+CjvgqcdCkeClrrsjhnQBmNgs4yd23xJ/fATQ3UzcHONbM+hC7uY8CLk86ZxowBngHGAG83lSGlLuvMbMvzexfiE2OXwnc30w7RHIi5/MaEZ8QB839lKp0a1UdDoRXje+KH2uUu+8GbgReBT4GnnP3hWZ2l5kF6z+eAMrNbAnwU0Jb0ZrZcmJFFK8ys1VmVhF/aSzwW2AJsBRNjEspmTspll0W4f3EF635Mrc1u6Tg0p0cfwp4z8xejD+/mNhaiia5+wxiKbPhY7eFHu8gVmU31Xt7N3J8Lnun6IqUhoiv2QgPSWmIqnSlW1Z9vJm9ApwRP1Tt7h/mrlkixenzL3eyfttOtudywjfCQ1SXD+7J5YOVxVjq0u1xBGs2mlq3IVIwizcupnpmNYs3Ls7pSvr123ayfefu3Pw2HZ4MFyliaQcOKZx0Cxnm6+ZZbMIlRfoe0jfnJUb267gPNdVZnvAND0tFcIhK2hYFjghIp5Bhvm+eGQmvSRgwAqqa2BejBb9157O8SM5UVTf930WkiChwRERzPYiivnmG1yRA4zfIIv2tO5jXGPnoO2z33ezXUf9spG3TvwDJj3R6EEX6W3cwr4HFhqkO7dyx0E0SKSgFjiIVntdoa3MWxSiY16ieqdIZYeE1Gyor0nakuwBQ8iw8r1F0cxbpKoHFbNK44ZXdE8FCazbaFvU4iljfQ/oy6bxJiedBDyQykhezhXf6i4CgSJ/mNVLTmo22Sz0Oya1gMVsW5y7Cace5FBTp07yGSEP6NSqigjmQtjb/ke+044ojDmA/jd2LNKDAEVHhoBHJ+Y8WKuq0Y5E2QoEjwpLnQERE8kFzHCIikhH1OIpAeM0GEM2hp0zKiohIpClwFIHwfEWuM4VyJt2yIkUqXFZkeYcvObRzR7oUulEiRUpDVUUimK+IdIZUtwGRLQmeKCsCbN+5m/Xbdha4RSLFS4FDJG6/jvtQc/03tNhPpBkKHCIikhEFDhERyYj65AWQvKOfiEiUKHAUQPKOfpGj1FuRNk2Bo0AinT1Vgqm3IpI+BQ5pmYim3ULDHf2279zN+kI3SCRiNDlehBZvXBztoay184t+Ayel3oq0nP7VFJnwZHkkJ84HjNjzONjAKYK279zNorXaDlUkFQWOIhP5suFV1ZGb80h2aOeOrAd6x7dD/dMXhW6RSHFR4BBJ0vWAjnQ9oCOTzvsGAH+aWeAGiRQZzXFIYcydVPTzICKSmnocBRZMgkc6PbclgpTeiM6DpCqFH+khRpEMKHAUUOQnwqFhBlWmKbrdBkD19Ny0KwNBVdxMMqxSlcJX4JC2QoGjgCI/Ed6SDKrWBJocCCbCg8eZCErhV8+MdjKASKYUOKTlMs2gKsJU3WAiXETSl9PAYWbnAb8G2gO/dfcJSa93BJ4CTgY2ACPdfXn8tZuBa4A64Ifu/mr8+HJgS/z4bnevyuV3kCwqgVRdEclh4DCz9sCDwLnAKmCOmU1z90Wh064BvnD3Y8xsFDARGGlmFcAooB9wJPCamR3n7nXx9w1x90hVighPpgZj4yIiUZTLdNxBwBJ3X+buu4BngeFJ5wwH/jv++AXgbDOz+PFn3X2nu38KLIlfL7LCFXH7HtI3mpPhSqEVEXI7VNUdWBl6vgoY3Ng57r7bzDYD5fHj7ya9t3v8sQN/NDMHHnX3x3LQ9pwIJlMjK+IptEFV3O07d6tGlUgrRPFfz+nuvtrMugJ/MrNP3H1W8klmdh1wHUDPnj3z3caEYIgqssNT4b03gqBRBCm0LREOGiqlLtJyuQwcq4GvhZ73iB9Ldc4qM9sHOJDYJHmj73X34OfnZvYisSGsvQJHvCfyGEBVVZVn4fu0SDhoRHJ4KtzLiGBPY8rsz5g6L/a/3XaPBQ0VLhRpnVwGjjnAsWbWh9hNfxRwedI504AxwDvACOB1d3czmwZMMbN7iU2OHwu8Z2adgXbuviX+eChwVw6/Q4ukmgiP3BBV0NOIeC9j6rzVLFoTq3KrnoZIduQscMTnLG4EXiWWjvukuy80s7uAue4+DXgC+J2ZLQE2EgsuxM97DlgE7Aa+7+51ZnY48GJs/px9gCnuXnQl6MK9jJLoaUSslwF7ehpB0Ki5/htUz1RPQyQbcjrH4e4zgBlJx24LPd4BpFw67e7jgfFJx5YBJ2S/pdmRPJ8R2V4GlFRPY3hl9+bfICJpi+LkeNHSfEbhpeppiEh2KXBkQcn0NCLeywD1NETyQYEjC0qqpxHBXkYy9TREckuBI0si2dMIi2BPI5xqO7yyO5cPLtx6HZG2RIEjQ+FUW23eU1jBsFRAgUMkPxQ4MhSuOQX537xn8cbFVM+sbvlK9HDm1IARka9Wq8V8IvmnwNEChSodEp4/afF8SjCfEYh44BCR/FPgiJCs7RhYBDvvZUrzGSLFQ4GjCRPfm8gnGz8BNJ9RaJrPECkeChxpCOY0FDgKS/MZIsVBgaMJNw26CYDqmZoHKHaL1nzJyEffSSz+E5HcyeUOgCUvnOEUSWvnl8SOfsMruyeChVaMi+SeehwtlJUMp0IKrxCP+Irxywf31JyHSB4pcLRQ1jKcCqWqOrKpuBqWEiksDVWlKfLDUnMnaVhKRLJCPY40RHZYKrxKfMVbsZ+9To/UsFRymXQNS4kUngJHGiI3LBUEjHCwCAJGBIanwov9Zn+6EYDBfQ5R70KkSChwlIqmehcRCBZh4R5GEDCi0MtQAUxpKxQ4oi7ivYuwqO/eV+gCmCL5osARVakCRoSDBZTGsFShCmCK5JMCR5SU0HBUEDDCwSIqw1JBr0JBQtoqBY4oCW/xGqGAkaqybTAkFZVgEQhn1EUmu04kyxQ4ilW4dwF7UmgjuMVrY5VtoziPEbkMO5EcUOAoVuHeRYQW7YV7F0BirkIrvEVKh1aOF7OgdxGhjZfCvYtFa75sEEREpDQocBRCUP5j0gWxxyUmGIJSL0OkNClwFEIwDLV2fsN5jKYE50dcuEChiEST5jgKJTz8lDwRDnvmN6BhbakI1ZkCEgEiuSChChSKRJcCR6Gtnd9wTUYgvEdGkZZAT54IBxqUOg8HhiDlNipptyLSOAWOfAn3KoLeRHJPIiLBIRBevBcI9yQUKERKkwJHNqQaagobMKJhem0QNIqkJ5FqgV6qld3JorZ4T0SyQ4EjU+EgEdz8w0EhWXhCO8+L9xpbU9HYKu5AlFd2i0juKXBkKggSgaDH0FhQmHRBftqVQrjSbDgwNLaKO1kUV3aLSO4pHbclguGmtfPT2461gKm0qdZUVBxxQCKYJKfGKl1WRJqT0x6HmZ0H/BpoD/zW3Sckvd4ReAo4GdgAjHT35fHXbgauAeqAH7r7q+lcM6teGbf3DT/VxHby87AWptI2NSmdrnCGU/Ac0kuNVbqsiDQmZ4HDzNoDDwLnAquAOWY2zd0XhU67BvjC3Y8xs1HARGCkmVUAo4B+wJHAa2Z2XPw9zV0z6/65ZQfrt+6MP+vJ1i7fZGnd2Uzddfyek94H3n9nr3kEOB64BYDhdd25nPSCQlOT0ukK3/zTTY3VXIaINCeXPY5BwBJ3XwZgZs8Cw4HwTX44cEf88QvAA2Zm8ePPuvtO4FMzWxK/HmlcM2vu3P1vLNr1JbP/secmPvvTjfEgMT9xLDD7042JG35jrzWXqRR+bzYnpZUaG5NqL43GjlXPrE48bm7vjeRrhN+fU/ZP2LUt9rhDZ5hZvefY5Kq9z6/dlt55eeK2m/p657uPGfXutGtnAIljyYJzqifv0+C9ktqR7Q7jf137p6xfN5eBozuwMvR8FTC4sXPcfbeZbQbK48ffTXpv8Ctzc9cEwMyuA64D6NmzdTfM8E08VepqIN3XlKlUGKn20mjuGMSCQVN7bxR0j47Oh+39OHwsWYfO6Z2XJ2Xt21FLPQDt2hll7WPTrsGxZOFzwu+V/DJ3z82FzUYA57n7tfHn/wYMdvcbQ+csiJ+zKv58KbFAcAfwrrs/HT/+BPBK/G1NXjOVqqoqnzt3bja/nohIyTOz9919ry5pLrOqVgNfCz3vET+W8hwz2wc4kNgkeWPvTeeaIiKSQ7kMHHOAY82sj5l1IDbZPS3pnGnAmPjjEcDrHusCTQNGmVlHM+sDHAu8l+Y1RUQkh3I2xxGfs7gReJVY6uyT7r7QzO4C5rr7NOAJ4Hfxye+NxAIB8fOeIzbpvRv4vrvXAaS6Zq6+g4iI7C1ncxzFRHMcIiKZK8Qch4iIlCAFDhERyYgCh4iIZESBQ0REMtImJsfNbB2wCdjcgrcfCqzPaoOkKQfSsr+nYlfM36sQbcvHZ2b7M7J1vdZep6Xvb8m9rJe771VioE0EDgAze8zdr2vB++amyiqQ3Gjp31OxK+bvVYi25eMzs/0Z2bpea69TDPeytjRU9VKhGyBpKdW/p2L+XoVoWz4+M9ufka3rtfY6Bf9/qc30OFpKPQ4RKQXqceTXY4VugIhIFmTtXqYeh4iIZEQ9DhERyYgCh4iIZESBQ0REMqLAkSEzO8rMnjCzFwrdFhGRljCzi83scTOrMbOhmb5fgQMwsyfN7PP4Vrbh4+eZ2WIzW2Jm4wDcfZm7X1OYloqIpJbhfez/uPv3gBuAkZl+lgJHzGTgvPABM2sPPAicD1QAl5lZRf6bJiKSlslkfh+7Jf56RhQ4AHefRWwHwrBBwJJ4D2MX8CwwPO+NExFJQyb3MYuZCLzi7h9k+lkKHI3rDqwMPV8FdDezcjN7BDjRzG4uTNNERNKS8j4G/AA4BxhhZjdketGc7Tleqtx9A7FxQRGRSHL33wC/aen71eNo3Grga6HnPeLHRESiIif3MQWOxs0BjjWzPmbWARgFTCtwm0REMpGT+5gCB2Bm/xt4B+hrZqvM7Bp33w3cCLwKfAw85+4LC9lOEZHG5PM+piKHIiKSEfU4REQkIwocIiKSEQUOERHJiAKHiIhkRIFDREQyosAhIiIZUeAQEZGMKHCItJKZnWVmL8cfXxTseZCF675gZkdl6VqvmdnB2biWiAKHSArxstMZ//tw92nuPiELn98PaO/uy1p7rbjfAWOzdC1p4xQ4ROLMrHd8p7SngAXA18zsYTOba2YLzezO0LnnmdknZvYBcEno+FVm9kD88WQzGxF6bWv85xFmNsvM5pnZAjM7I0VzRgNTk98bfzzCzCaHPuNhM3vXzJbFez9PmtnHwTlx04DLWvUfSCROgUOkoWOBh9y9n7uvAH7h7lXAQOB/mNlAM+sEPA5cCJwMdMvwMy4HXnX3SuAEYF6Kc04D3k/zegcD3wB+QixA3Af0AwaYWSWAu38BdDSz8gzbKrIXBQ6Rhla4+7uh59+N9yo+JHYzrgCOBz519//nsWJvT2f4GXOAajO7Axjg7ltSnHMEsC7N670Ub8d84J/uPt/d64GFQO/QeZ8DR2bYVpG9KHCINLQteGBmfYCfAWe7+0BgOtApg2vtJv5vLD5f0gESW3yeSWxfhMlmdmWK936V9FnhaqTJbdgZ/1kfehw8D2/W1il+XZFWUeAQadwBxALJZjM7HDg/fvwToLeZHR1/3tjcwXJiQ1kAFwFlAGbWi1jP4HHgt8BJKd77MXBM6Pk/zezr8QD07Uy/iJkZsSG15Zm+VySZto4VaYS7f2RmHxILFCuBt+PHd5jZdcB0M9sOvAnsn+ISjwNTzewjYCZ7ejNnAT83s1pgK5CqxzE9ft5r8efjgJeJDV/NBbpk+HVOBt6N788g0iraj0OkCJnZvsBfgNPcvS4L1/s1MM3d/9zqxkmbp6EqkSLk7l8BtwPds3TJBQoaki3qcYiISEbU4xARkYwocIiISEYUOEREJCMKHCIikhEFDhERycj/B1cf8oirle8AAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Exponential rate, AlwaysN fragments" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 5, + "source": [ + "from PySDM.initialisation.spectra import Exponential\n", + "from PySDM.physics.spontaneous_breakup_rates import Expon\n", + "from PySDM.physics.spontaneous_breakup_fragmentations import SpontaneousAlwaysN\n", + "from PySDM.physics.constants import si\n", + "from PySDM.physics.formulae import Formulae\n", + "#from pystrict import strict\n", + "\n", + "#@strict\n", + "class Settings:\n", + "\n", + " def __init__(self):\n", + " self.formulae = Formulae()\n", + " self.n_sd = 2**8\n", + " self.n_part = 100 / si.cm**3\n", + " self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres)\n", + " self.dv = 1 * si.m**3\n", + " self.norm_factor = self.n_part * self.dv\n", + " self.rho = 1000 * si.kilogram / si.metre**3\n", + " self.dt = 1 * si.seconds\n", + " self.adaptive = False\n", + " self.seed = 44\n", + " self._steps = [0, 1, 2]\n", + " self.breakup_rate = Expon(1.58e-4 / si.s, 0.773 * si.mm)\n", + " self.fragmentation = SpontaneousAlwaysN(n=3)\n", + " self.spectrum = Exponential(norm_factor=self.norm_factor, scale=self.X0)\n", + " self.radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(100 * si.um), num=128, endpoint=True)\n", + " self.radius_range = [0 * si.um, 1e6 * si.um]\n", + "\n", + " @property\n", + " def output_steps(self):\n", + " return [int(step/self.dt) for step in self._steps]" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 6, + "source": [ + "settings = Settings()\n", + "backend = CPU\n", + "\n", + "builder = Builder(n_sd=settings.n_sd, backend=backend, formulae=settings.formulae)\n", + "builder.set_environment(Box(dv=settings.dv, dt=settings.dt))\n", + "attributes = {}\n", + "attributes['volume'], attributes['n'] = ConstantMultiplicity(settings.spectrum).sample(settings.n_sd)\n", + "breakup = SpontaneousBreakup(settings.breakup_rate, settings.fragmentation, adaptive=settings.adaptive)\n", + "builder.add_dynamic(breakup)\n", + "\n", + "products = [ParticlesVolumeSpectrum(), WallTime(), ParticleMeanRadius(), ParticlesConcentration(radius_range = settings.radius_range)]\n", + "core = builder.build(attributes, products)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 7, + "source": [ + "for step in settings.output_steps:\n", + " #print(step, flush=True)\n", + " core.run(step - core.n_steps)\n", + " pyplot.step(x=settings.radius_bins_edges[:-1] / si.micrometres, \n", + " y=core.products['dv/dlnr'].get(settings.radius_bins_edges) * settings.rho,\n", + " where='post', label=\"t = {step*settings.dt}s\")\n", + " print(core.products['radius_m1'].get(), core.products['n_a_cm3'].get())\n", + " \n", + "pyplot.xscale(\"log\")\n", + "pyplot.xlabel(\"radius (um)\")\n", + "pyplot.ylabel(\"dm/dlnr\")\n", + "pyplot.legend()" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[27.26101585] [99.997952]\n", + "[23.15512423] [199.995904]\n", + "[17.18684981] [499.98976]\n" + ] + }, + { + "output_type": "execute_result", + "data": { + "text/plain": [ + "" + ] + }, + "metadata": {}, + "execution_count": 7 + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApD0lEQVR4nO3deZxU9Znv8c9DS4OKuLRREERQEW0W27GF8apcjYoEo6jRQMRRUYNeY5JJJpmQqKhcSSTXGyfGBTEC0YQRIfGKoZWMmlzUUbRVvCyKAwSkEZRFZJOtee4fdao5XV3dXdVdp2v7vl+vflWdU79z6lcs9fRve37m7oiIiKSqXbYrICIi+UWBQ0RE0qLAISIiaVHgEBGRtChwiIhIWhQ4REQkLQdkuwJt4cgjj/SePXtmuxoiInnlnXfe2eDuX0k8XxSBo2fPnlRXV2e7GiIiecXMViU7r64qERFJiwKHiIikRYFDRETSUhRjHCL5bs+ePdTU1LBz585sV0UKUMeOHenevTvt27dPqbwCh0geqKmp4ZBDDqFnz56YWbarIwXE3dm4cSM1NTX06tUrpWvUVSWSB3bu3ElZWZmChmScmVFWVpZWa1YtDpHGVE+FhbNiz/tfCZWjs1odBQ2JSrr/ttTiEGnMwlmwbiGseg3+/M8w9eLYT/XUbNeszW3evJlHHnkko/esqqri5JNP5vbbb2+y3M9//vNWv9e0adNYuXIlLdl/aMGCBVRVVdUd/+1vf+M///M/644nTZrEk08+2eo6pqpnz55s2LAh6d/JiBEj6N+/P4sXL460DgocIk3p0h++/m9w3Nmx43UL97dCikgUgePxxx/nV7/6FRMmTGiyXGsCx5o1a7jppptYvXo1r732Grfcckva92gucNxyyy1ce+21La5jSyX7O5kxYwajRo3i6aefjvS9FThEmlM5GkbPif106Z/t2mTF2LFjWb58ORUVFfz4xz/OyD03b97MUUcdVXe8du1aBg8eTEVFBf369ePVV19l7NixfPnll1RUVDBq1CgAfv/73zNw4EAqKiq4+eabqa2tBaBTp0784Ac/oG/fvpx//vmsX7+ebt26MWHCBJ544gmefvppHn30UQAefPBBysvLGTBgACNHjgRg+/bt3HDDDQwcOJDTTjuN5557jt27dzNu3DhmzJhBRUUFEydOZNKkSTzwwANUVFTw6quvcvfdd3P//fcDcO655/KTn/yEgQMHctJJJ/Hqq68CsGPHDr75zW9SXl7O5ZdfzqBBg6iurqa2tpbrr7+efv360b9/fx544IEGf04bN25kyJAh9O3bl5tuuqmu1dTY30mXLl3YvHlzvXsk+7ytoTEOkTxzz/OLWfLJlozes/yYztx1Sd9GX7/vvvtYtGgRCxYsSPr6Oeecw9atWxucv//++7nggguSXlNbW0u7dvt/d50+fToXXXQRt99+O7W1tezYsYNzzjmHhx56qO59P/jgA2bMmMHrr79O+/btufXWW/nDH/7Atddey/bt26msrOSBBx5g/Pjx3HPPPfzsZz/jrrvu4oYbbqBXr1585zvf4dFHH+W+++7j73//Ox06dKj7kp0wYQJf/epXmTJlCps3b2bgwIFccMEFjB8/nurqah566CEAvvzySzp16sSPfvQjAF5++eV6n2vv3r289dZbVFVVcc899/DSSy/xyCOPcPjhh7NkyRIWLVpERUUFEGvNrFmzhkWLFgE0+MIHuOeeezj77LMZN24cc+bM4Yknnmjy76Rdu3Z1wTQu2edtDQUOEWm1+G/Wqdq+fTvLly/nmGOOqTt3xhlncMMNN7Bnzx4uu+yyui/XsJdffpl33nmHM844A4h9icdbLe3atWPEiBEAXHPNNVxxxRUcc8wxPP7440ybNo1zzjmHa665BoABAwYwatQoLrvsMi677DIA/vKXvzB79uy61sPOnTv5+OOP0/pcAFdccQUAp59+OitXrgTgtdde4/vf/z4A/fr1Y8CAAQAcf/zxrFixgu9+97tcfPHFDBkypMH95s2bx5/+9CcALr74Yg4//PAm379bt24sXryYffv21QXmZJ+3NRQ4RPJMUy2DbEmnxTF//nyGDRvGkCFD6NKlS935wYMHM2/ePObMmcP111/PD3/4wwZjB+7Oddddxy9+8Ytm6xSeKXT99dfXe23OnDnMmzeP559/ngkTJrBw4ULcnT/+8Y/06dOnQX3T0aFDBwBKSkrYu3dvk2UPP/xw3n//febOncukSZN45plnmDJlSlrvl2jw4MFs3bqVHj168P7771NWVpb08x5wQMu//jXGISLNOuSQQ5IGhrhXX32VBQsWNPhJ1k01aNAgPvnkE958802WL19ed37VqlUcffTRfPvb3+amm27i3XffBaB9+/bs2bMHgPPPP59Zs2bx2WefAbBp0yZWrYolcN23bx+zZsUmLkyfPp2zzz47aV337dvH6tWrOe+885g4cSJffPEF27Zt46KLLuI3v/lN3RjCe++9l/SzN/dnkcxZZ53FM888A8CSJUtYuHAhABs2bGDfvn184xvf4N577637zGGDBw9m+vTpALzwwgt8/vnnTdajqqqKo48+mtWrV1NWVtbo520NBQ4RaVZZWRlnnXUW/fr1y8jgeIcOHTjuuOPqvgQhNlvp1FNP5bTTTmPGjBl1XTtjxoyp62opLy/n3nvvZciQIQwYMIALL7yQtWvXAnDwwQfz1ltv0a9fP1555RXGjRuX9L1ra2u55ppr6N+/P6eddhrf+973OOyww7jzzjvZs2cPAwYMoG/fvtx5550AnHfeeSxZsoSKigpmzJjBJZdcwrPPPls3OJ6KW2+9lfXr11NeXs4dd9xB3759OfTQQ1mzZg3nnnsuFRUVXHPNNXUtqUmTJjFp0iQA7rrrLubNm0ffvn3505/+RI8ePYDG/04+//xzTjjhhLoWV2OftzWsJfOa801lZaVrPw5J29SLY4+j5zR9rg188MEHnHLKKW36nlEbNmwYt912G8OGDcvI/Tp16tTq36SjUltby549e+jYsSPLly/nggsuYOnSpZSWlmb8vX75y1+yceNGJk6cmNZ1yf6Nmdk77l6ZWFZjHJIXZn40k6oVsbn0w44fxlUnXZWZG+fY6vBicsstt3DHHXfw+uuvN7uWI9/t2LGD8847jz179uDuPPLII5EEjZEjR7Jq1aq6acdRUeCQvFC1ooqlm5bWHWcscMRXh8cpcLSZSy+9lEsvvTRj98vV1gbExiPaotcj6oV/cQockjf6HNGn+UItUaSL+kRaSoPjIiKSFgUOERFJiwKHiIikRYFDRJqltOpKqx6mwCEizVJadaVVD1PgEJFmKa260qqHaTquSL55YWz9tSeZ0KU/fO2+Rl9WWnWlVQ+LNHCY2VDg10AJ8Ft3vy/h9Q7Ak8DpwEZghLuvNLMLgfuAUmA38GN3fyW45m9AV+DL4DZD3P2zKD+HFIl1C/enFOl/ZXbrkmeUVn0/pVVvBTMrAR4GLgRqgLfNbLa7LwkVuxH43N1PNLORwERgBLABuMTdPzGzfsBcoFvoulHuruRTkjnhQJHp3+YzrYmWQbYorfp+SqveOgOBZe6+wt13A08DwxPKDAd+FzyfBZxvZubu77n7J8H5xcCBQetEJBqJ28OuW5j7AaQNKa260qqHRRk4ugGrQ8c11G811Cvj7nuBL4CyhDLfAN51912hc1PNbIGZ3WnhXytEMqH/lbHg0aW/uqwCSquutOphkaVVN7MrgaHuflNw/E/AIHe/LVRmUVCmJjheHpTZEBz3BWYTG8dYHpzr5u5rzOwQ4I/A7929wSRqMxsDjAHo0aPH6fHfSiQ/jX6xfvLBqUOnZubG6aZJV1r1jFFadaVVT2YNcGzouHtwLlmZGjM7ADiU2CA5ZtYdeBa4Nh40ANx9TfC41cymE+sSaxA43H0yMBli+3Fk6DOJSIYorbrSqifzNtDbzHoRCxAjgasTyswGrgPeAK4EXnF3N7PDgDnAWHd/PV44CC6HufsGM2sPfB14KcLPICIRUVr1zGurtOqRjXEEYxa3EZsR9QHwjLsvNrPxZhb/1/IEUGZmy4AfAmOD87cBJwLjgrGMBWZ2FNABmGtm/w9YQCwgPR7VZxARkYYiXcfh7lVAVcK5caHnO4EGO/K4+73AvY3c9vRM1lFERNKjlCMiIpIWpRyR4qH9xUUyQi0OKR7x/cXXLdwfQCQlSquutOphChzSZmZ+NJPRL46u+5n50cy2r0R8YZ+kRWnVlVY9TIFD2kzViiqWbloKwNJNS6laUdXMFekLB6esBKYCpbTqSqsepjEOaVN9jujD1KFTG6wEz5RwcAK46qQGk/by3sS3JvLhpg8zes+TjziZnwz8SaOvK6260qqHKXBIwelzRJ/mC0lGKa36fkqrLiI5p6mWQTo2btvF5i9jWWcPO7B9q+6ltOr7Ka26iNQf1LdPmUnuprZIx+Yv97BzTy0799TWBZDGKK260qqHKXCINKPeoD67qbLtWa5R5nRsX0LH9iXNllNadaVVD4ssrXouqays9LZIMCZNiw+IhwfHU02Pnmpa9Qbljhy8f83GuoX7p+LGn8cfm0iTXq/e02IZpqde37b/nqJIq758ff3fOk/4SqeM3r85SquutOoiuSm+6K+xjZm0WVPWKK260qpLEZr50cx6azGGHT8sN6e/JmtRKN1I1imteublfVp1KXxtsaBPRHKPAoe0SnxBXybWTmjVd9OKYTxSsiPdf1vqqpKcUQyrvluqY8eObNy4kbKysnrrE0Ray93ZuHEjHTt2TPkaBQ7JKVr1nVz37t2pqalh/fr1Gbvn+q276h3v3tAhY/eW/NKxY0e6d++ecnkFDpE80L59e3r16pXRe9792Bv1jmfcXNGgzPT5H/PcgjV1x8MrunH1oB4ZrYfkH41xSGHaum7/3hvSYs8tWMOStVsAWLJ2S70gIsVLgUMK0/b1sHu71mlkQHnXzsy4+UzKu3bOdlUkR6irSiITXucx7PjMrA5OS+nBcH3jK8JFpGUUOCQyibOkMi0xMGkWlkjbUOCQSEU5S0rTd0WyQ4FDsircali6aWnagabYp++GZz01NeMp1XIiqVDgkIxKNxDEWw19juhDnyP6MOz4YVlJXZLRvFvVU/dn5IXY4HxEubHCs56ARgNCquVEUqFZVZJR4e6jeCBoTjxtydShU7PW3ZTRvFvxjLwQewwHkQiUd+2c0oynVMuJNCfSwGFmQ81sqZktM7OxSV7vYGYzgtfnm1nP4PyFZvaOmS0MHr8auub04PwyM3vQlH8h5+RCIGiJTObdqsvIG9//Q6SARBY4zKwEeBj4GlAOfMvMyhOK3Qh87u4nAg8A8Z1HNgCXuHt/4DrgqdA1jwLfBnoHP0Oj+gySh+IL/3YXzi59IrkmyjGOgcAyd18BYGZPA8OBJaEyw4G7g+ezgIfMzNz9vVCZxcCBZtYBOALo7O5vBvd8ErgMeCHCzyGBrK/LSEV84V/pwXDwV7JdG5GCFGXg6AasDh3XAIMaK+Pue83sC6CMWIsj7hvAu+6+y8y6BfcJ37NbpisuyUW9LiNjSg9WF5FIhHJ6VpWZ9SXWfTWkBdeOAcYAdZu7S+tlcvpreBBdRPJHlIFjDXBs6Lh7cC5ZmRozOwA4FNgIYGbdgWeBa919eah8OPdvsnsC4O6TgckAlZWV2gEnx4S7unK22yvPhdduLFm7RTOqJGOiDBxvA73NrBexL/eRwNUJZWYTG/x+A7gSeMXd3cwOA+YAY9399Xhhd19rZlvM7B+B+cC1wG8i/AwSkatOuiqvZlylJL5+Y93CnOgqi6/diE/DHV7RrcnstslSqIskE1ngCMYsbgPmAiXAFHdfbGbjgWp3nw08ATxlZsuATcSCC8BtwInAODMbF5wb4u6fAbcC04ADiQ2Ka2BcckM4aORIRt54Ztu4pgJHONCEFwuKJIp0jMPdq4CqhHPjQs93Ag1+7XT3e4F7G7lnNdAvszUVyZD4+o08FQ80IxI2eRIJ08pxyZilm5bmx6wrEWmVnJ5VJdmXau4pDXaLFA8FDmlSqkkIC3KwO8do8FpyhQKHNCuewykuG9lrRYPXkjsUOKToJEuhni9SGbyOt0y0dkOiosAhRSfc/ZaRwfwcXr+RandWvAWjQCOpUOCQ/LV1XSypYdzUi2FPkOCwGfHut9EvZmCDpTxYv9GUcHDRuImkQoFD8lc4E25ctrLi5vH6jasH9dCOgJIWBQ7Jb+FMuEOnQiZaECLSJAUOSSo+gJzKvuEi4anCwyu6qQVT4LRyXJIKB41MzDrSqvKWmz7/Y0Y89kazU3CXrN2SUrkoxAfkl6zd0mQ+LCkManFIoxLXb7RUQa4qj8+kgshnU6UySyp8Pp3ZVJmkGVnFQ4FDIleQq8rDM6naYDZVc7OkNMAtbUmBQ6SlWjuTKtxqgVjwqdTgvuQ+jXFIi2jMIgPirRaIPYaDiEgOU4tD0taWYxYFP7sr3mqZenG2ayKSMgUOSVtbjlkkzu4qlgSL2i9ccpm6qiTnxWd35cwA++7t+7uYIhKfSQXZmyUl0phmWxxmZkB3d1/dBvWRLEp106aiFk9n0uX4rM+kEsmWZgOHu7uZVQHZT/spkUq2aVPOiM9Ask9jx2kkNMyoQ7rEfjKwviVbwgsF1QUmLZHqGMe7ZnaGu78daW0k6zK16K+1lm5ayugXR+9v+cRnIHU9an+hbCU0zGO5sFBQ8l+qgWMQMMrMVgHbASPWGBkQWc2kzeTazKVwS6eu5bN2WrDYLggc+ZzQMIv7d2RioaDyUkmqgeOiSGshWZXpvFStlXTW1uvTslKXSGRz/44MLDoMD9wDbRo4ku27rsDV9lIKHO6+ysxKgKNTvUbyS650UWVazm4T28iq83S2fW3Rrn3hoBWfGdaC1erZGhtJtu+6AkfbSykImNl3gbuAT4F9wWkH1FUlOS3j28RGLNVtX1u1a1+eLzpMZd91iVaqrYfvA33cfWOUlRGJQka3iW0DqUzDVVJDyaZUFwCuBr5I9+ZmNtTMlprZMjMbm+T1DmY2I3h9vpn1DM6XmdlfzWybmT2UcM3fgnsuCH6OSryviIhEJ9UWxwrgb2Y2B9gVP+nuv2rsgmBM5GHgQqAGeNvMZrv7klCxG4HP3f1EMxsJTARGADuBO4F+wU+iUe5enWLdRbKvDffviFKLxlWk4KTa4vgY+A+gFDgk9NOUgcAyd1/h7ruBp4HhCWWGA78Lns8Czjczc/ft7v4asQAikv/CmXCzMZsqA4ZXdKO8a2et/5CUZ1Xd04J7dyPWxRVXQ2w9SNIy7r7XzL4AyoANzdx7qpnVAn8E7nV3b0H9RNpWE/t3fLp1Jxu27WLJ7txdzZ04rqItYotXk4HDzJ4nNnsqKXe/NOM1at4od19jZocQCxz/BDyZWMjMxgBjAHr00CBiomLISxWfRZUPn23Dtl3s2F2r3+YlLzTX4ri/FfdeAxwbOu4enEtWpsbMDgAOBZqcueXua4LHrWY2nViXWIPA4e6TgckAlZWVapEkyOm8VGEtXGWdj/ucH1RaoqSGkheaDBzu/n9bce+3gd5m1otYgBgJXJ1QZjZwHfAGcCXwSlPdTkFwOczdN5hZe+DrwEutqGNRy4tFf4mrrDfMS+mynNrnPDy2IVIAmuuqWkjTXVWNLgAMxixuA+YCJcAUd19sZuOBanefDTwBPGVmy4BNxIJL/L1XAp2BUjO7DBgCrALmBkGjhFjQeDyFzymBXMtLlZLw2MCLqQWOnBEeBM/DAXGRZJrrqvp68Pid4PGp4PEamggoce5eBVQlnBsXer4TSPprobv3bOS2pzf3vtK4XMtLVfAqRzea0iOcd+lHu2s5qLSkLWsm0mLNdVWtAjCzC939tNBLPzGzd4EGi/ok9+VFF1URCKcXOai0hCM7dch2lURSkuoCQDOzs9z99eDgv6FtZyUqBbJYLhV16UWmHprtqoikLNXAcSMwxczi/7o3AzdEUiOR8IB4ni6WEylkzQ2Onwm86e7vAKfGA4e7p523SiQtTSyWk7ah9CLSmOZaHNcCD5vZR8CLwIvuvi76aokUrnT23MiWVqVtl4LX3OD4/wAws5OBrwHTglbHX4kFktfdvTbyWkqr5eU03AKV6p4b2aS07dKUVHNVfQh8CDxgZgcC5xGbRvsroDK66kmm5MU03Czuxd3WUtlzIy0Z2BI227SXef5IeRtYMzucWHqQA4B1wFR3/25UFZPMy/lpuNncizvfZWhL2GzK5l7mkp5Ut479n8D1xPblCG8d+9VoqiWZkJeJDDUo3nJ5viUsaCA+X6Ta4vgmcEKwr4bkibxJZFgEwt0wuTwoLpKKVAPHIuAw4LPoqiJRyPnuqSIRHhDP5UFxkVSkGjh+AbxnZouov3VsNvbjkEJSjKvERfJcqoHjd8T2A1/I/jEOkdbTKnGRvJNq4Njh7g9GWhMpXhoQF8krqSYqfNXMfmFmZ5rZP8R/Iq2ZSBqWblrK6BdH120Xmyumz/+YEY+9UW+aqUi+S7XFEU+pPih4NDQdN2cV2yrx8GyxXJs9lg+rxEXS1VySwx8GT/9MLFBY6GXt452jim2VeE5tE5tEJIPi4UkFeTgupFXi+a25FschwWMf4AzgOWLB4xLgrQjrJa2U89NwtUq8deJ/fnlKq8TzW3NJDu8BMLN5wD+4+9bg+G5Ao5nSOgU6KN5mi/3yfOqyFkHmr1THOI4GwqvGdwfnJEfkZXqRAqXFflLoUg0cTwJvmdmzwfFlwLQoKiQto/QiuUWL/aSQpZpWfYKZvQCcE5wa7e7vRVctaYmcH9cQkYKQclp1d38XeDfCuoikJLxmo6i65PI8PYsSPRaOVBcASo6a+dHMnFz4FpVhxw+rCxa51iUX+WK/8EyqPJyNFp5JpbGf/JZyi0NyU16s2cigXF6z0SaL/RJnooV3/csDiWM/8RaI5BcFjgKQN2Mbed7VkgoNiksxUFeVtJ0872rJmuqpsV398nTBn/J1FZ5IA4eZDTWzpWa2zMzGJnm9g5nNCF6fb2Y9g/NlZvZXM9tmZg8lXHO6mS0MrnnQzCzxvsUgb8c24l0to+fk3Z7YjWmzsY08DbbK11V4IuuqMrMS4GHgQqAGeNvMZrv7klCxG4HP3f1EMxtJbM+PEcBO4E6gX/AT9ijwbWA+UAUMBV6I6nPkqmIb28hlWRnbyDPqwissUY5xDASWufsKADN7GhgOhAPHcODu4Pks4CEzM3ffDrxmZieGb2hmXYHO7v5mcPwkscWIRRc4II/GNoqAvhilmEQZOLoBq0PHNexPy96gjLvvNbMvgDJgQxP3rEm4Z9Jf8cxsDDAGoEcPJVCTzIp8TUIBTCSI/xlpzUbhKdhZVe4+GZgMUFlZWRAp4JWPKndEno+qALbU1dhG4YoycKwBjg0ddw/OJStTY2YHAIcCG5u5Z/dm7lmwlI8qt0TePZXn4xqgLrxCFWXgeBvobWa9iH25jwSuTigzG7gOeAO4EnjF3RttHbj7WjPbYmb/SGxw/FrgN1FUPldpXKPAZXCDq2xZsnZL3SwzdVEVpsgCRzBmcRswFygBprj7YjMbD1S7+2zgCeApM1sGbCIWXAAws5VAZ6DUzC4DhgQzsm4llpn3QGKD4kU5MC7Z8dmWXWzYvosdUX0p5vnU23CXlLqoClekYxzuXkVsymz43LjQ851A0vwR7t6zkfPVNJyiK9ImNmzfxY5de6P5UgwvjszTLqqrB/XQbn5FoGAHx0WiclCHA5gxOsP99uHWRR62NKS4KHCI5ILK0QWzkl4KnwKHSDPi4xojHnuDHb6Xgzrov40UNyU5FGlGfFwDYt1URx7cIcs1EskuBQ6RFBzU4QBm3Hwm5V07c1RnBY648NRbKR5qc4s04sdzH2PeJ39hh3/MQaaZQok09bZ4KXCINCIcNAYfMyTb1ck5mnpbvBQ4RJpwkPVg/ug/ZrsaIjlFYxwiIpIWBQ4REUmLAoeIiKRFYxwi1F/kt7J0i9ZqiDRBgUOE0CI/gx279ja6BaWIqKtKpE58kZ9Siog0TYFDRETSosAhIiJpUZtcio4GwkVaR4FDio4GwkVaR11VUpQ0EC7ScgocIkns8I/Z4R9nuxoiOUmBQyTB4GOGcJD1UFZckUaonS6S4H9ddDNwc7arIZKz1OIQEZG0KHCIiEhaFDjywMyPZjL6xdEs3bQ021UpSBoIF0mPAkceqFpRxdJNS+lzRB+GHT8s29VJX/VUmHoxrFuY7Zo0oIFwkfRFOjhuZkOBXwMlwG/d/b6E1zsATwKnAxuBEe6+Mnjtp8CNQC3wPXefG5xfCWwNzu9198ooP0Ou6HNEH6YOnZrtarTMwlmxoNGlP/S/Mtu1qUcD4SLpiyxwmFkJ8DBwIVADvG1ms919SajYjcDn7n6imY0EJgIjzKwcGAn0BY4BXjKzk9y9NrjuPHfXgt980qU/jJ6T7VqISAZE2VU1EFjm7ivcfTfwNDA8ocxw4HfB81nA+WZmwfmn3X2Xu/8dWBbcT6TFPtuyi3fXLdZ4hkgrRRk4ugGrQ8c1wbmkZdx9L/AFUNbMtQ78xczeMbMxjb25mY0xs2ozq16/fn2rPogUhj1bTmXfzq4azxBppXxcAHi2u68xs6OA/zCzD919XmIhd58MTAaorKz0tq6k5Ibp8z/muQVrAFi7toLyroOZMfrMLNdKJL9FGTjWAMeGjrsH55KVqTGzA4BDiQ2SN3qtu8cfPzOzZ4l1YTUIHFLc4gFj/t83ATCo1xGUd+3M8IrERq+IpCvKwPE20NvMehH70h8JXJ1QZjZwHfAGcCXwiru7mc0GppvZr4gNjvcG3jKzg4F27r41eD4EGB/hZ5A89dyCNSxZu4VBvY5geEU3rh7UI9tVEikYkQUOd99rZrcBc4lNx53i7ovNbDxQ7e6zgSeAp8xsGbCJWHAhKPcMsATYC3zH3WvN7Gjg2dj4OQcA0939xag+g+SfeEtjydotlHftzIyb1S0lkmmRjnG4exVQlXBuXOj5TuCqRq6dAExIOLcCODXzNZVCEQ4a6pYSiUY+Do6LNEktDZFoKXBI3grPmNI4hkjbUeCQvBXvlopT4BBpGwocktfKu3bOdhVEio4Ch+QFdUuJ5A4FDskL6pYSyR0KHJI31C0lkhu0kZMUhCVrtzDisTfqtUpEJBpqcUjeCy/008I/kegpcEjeCbcuyrt25upBPTTmIdKGFDgkr6h1IZJ9ChyS0xKTFqp1IZJ9ChySc8JrNsL7aah1IZIbFDgkZyTbfEn7aYjkHgUOyYpwqyIusXWhYCGSmxQ4JCvC4xZxChgi+UGBQyL36dadbNi2i/GPvVEvMGjfDJH8pMAhGTV/5v+m0389W3e8rffldNq2ix27a+tWdatFIZLflHJEMqrTfz3LsbuXA3Ds7uV1QeSg0hLlmhIpEAocknGrS0+g789eY3XpCQ1eU04pkfynriqJVLz1sbr0BK36FikQChzSYonjGRALFPGWxrbel7M6eH1b78u16lukQChwSJOSBYe4QbsXArC4tH/dudWlJ7Ct9+Wx16/6F+BfIq+jiLQtBQ6pFxy29b6cQVf9S925ZMEhbnFp/7ryIlI8FDgKWLKpsUCDIBGeCRXrWtp/TsFBRBIpcBSw+Jf/6tITYo9BwEgMEkDSGVDx2VEiImGRBg4zGwr8GigBfuvu9yW83gF4Ejgd2AiMcPeVwWs/BW4EaoHvufvcVO5ZbJoag4gHjb4/e43FPz+73gyn+Ovx802dExEJiyxwmFkJ8DBwIVADvG1ms919SajYjcDn7n6imY0EJgIjzKwcGAn0BY4BXjKzk4Jrmrtnxv3y38dQvfWtuuPKQwYC1DuX6mv/+q3JDe7XnFXtazluTwmLf352g9eaGoMID1QnznAC6o7D5ZKdExEJM3eP5sZmZwJ3u/tFwfFPAdz9F6Eyc4Myb5jZAcA64CvA2HDZeLngsibvmUxlZaVXV1en/Rn++bcX8sm+9XzQoRaAU3aV1D2PO2VXSd3zVF4L3yP8enPO2lHKBTtKk76mMYi2M/rF0SzdtJQ+R/TJdlVi1i2E3dtjz0sPhi79958rPbhh+fj55sq1ke2797Jvn9OundU9AvWeh8XPH1x6QL1rJblj2n2Ff7vpP1p8vZm94+6Vieej7KrqBqwOHdcAgxor4+57zewLoCw4/2bCtfHVYs3dEwAzGwOMAejRo3VrB07ZVRJrLVxXv7UQPxeX6mvh+0l+GXb8sGxXob6Dv9LwefhcotKDUyvXRtqXtGMP+wBo185oXxJLZhE/lyhcJnyttK2CHRx398nAZIi1OFpyj2SR+l+/1fiXfUtfk/xx1UlXcdVJV2W7GiJZFWWuqjXAsaHj7sG5pGWCrqpDiQ2SN3ZtKvcUEZEIRRk43gZ6m1kvMyslNtg9O6HMbOC64PmVwCseG3SZDYw0sw5m1gvoDbyV4j1FRCRCkXVVBWMWtwFziU2dneLui81sPFDt7rOBJ4CnzGwZsIlYICAo9wywBNgLfMfdawGS3TOqzyAiIg1FNqsql7R0VpWISDFrbFaV9uMQEZG0KHCIiEhaFDhERCQtChwiIpKWohgcN7P1wGbgixZcfiSwIaMVkqYcSsv+nnJdLn+ubNStLd4z0++Rqfu19j4tvb4l32XHuXuDFANFETgAzGyyu49pwXXVyWYVSDRa+veU63L5c2Wjbm3xnpl+j0zdr7X3yYXvsmLqqno+2xWQlBTq31Muf65s1K0t3jPT75Gp+7X2Pln/t1Q0LY6WUotDRAqBWhxtS9kJRaQQZOy7TC0OERFJi1ocIiKSFgUOERFJiwKHiIikRYEjTWZ2vJk9YWazsl0XEZGWMLPLzOxxM5thZkPSvV6BAzCzKWb2mZktSjg/1MyWmtkyMxsL4O4r3P3G7NRURCS5NL/H/o+7fxu4BRiR7nspcMRMA4aGT5hZCfAw8DWgHPiWmZW3fdVERFIyjfS/x+4IXk+LAgfg7vOI7UAYNhBYFrQwdgNPA8PbvHIiIilI53vMYiYCL7j7u+m+lwJH47oBq0PHNUA3Myszs0nAaWb20+xUTUQkJUm/x4DvAhcAV5rZLeneNLI9xwuVu28k1i8oIpKX3P1B4MGWXq8WR+PWAMeGjrsH50RE8kUk32MKHI17G+htZr3MrBQYCczOcp1ERNIRyfeYAgdgZv8OvAH0MbMaM7vR3fcCtwFzgQ+AZ9x9cTbrKSLSmLb8HlOSQxERSYtaHCIikhYFDhERSYsCh4iIpEWBQ0RE0qLAISIiaVHgEBGRtChwiIhIWhQ4RFrJzM41sz8Hzy+N73mQgfvOMrPjM3Svl8zs8EzcS0SBQySJIO102v8/3H22u9+XgffvC5S4+4rW3ivwFHBrhu4lRU6BQyRgZj2DndKeBBYBx5rZo2ZWbWaLzeyeUNmhZvahmb0LXBE6f72ZPRQ8n2ZmV4Ze2xY8djWzeWa2wMwWmdk5SaozCngu8drg+ZVmNi30Ho+a2ZtmtiJo/Uwxsw/iZQKzgW+16g9IJKDAIVJfb+ARd+/r7quA2929EhgA/HczG2BmHYHHgUuA04Euab7H1cBcd68ATgUWJClzFvBOivc7HDgT+AGxAPEA0Bfob2YVAO7+OdDBzMrSrKtIAwocIvWtcvc3Q8ffDFoV7xH7Mi4HTgb+7u7/5bFkb79P8z3eBkab2d1Af3ffmqRMV2B9ivd7PqjHQuBTd1/o7vuAxUDPULnPgGPSrKtIAwocIvVtjz8xs17Aj4Dz3X0AMAfomMa99hL8HwvGS0qhbovPwcT2RZhmZtcmufbLhPcKZyNNrMOu4HFf6Hn8OLxZW8fgviKtosAh0rjOxALJF2Z2NPC14PyHQE8zOyE4bmzsYCWxriyAS4H2AGZ2HLGWwePAb4F/SHLtB8CJoeNPzeyUIABdnu4HMTMj1qW2Mt1rRRJp61iRRrj7+2b2HrFAsRp4PTi/08zGAHPMbAfwKnBIkls8DjxnZu8DL7K/NXMu8GMz2wNsA5K1OOYE5V4KjscCfybWfVUNdErz45wOvBnszyDSKtqPQyQHmdmBwF+Bs9y9NgP3+zUw291fbnXlpOipq0okB7n7l8BdQLcM3XKRgoZkilocIiKSFrU4REQkLQocIiKSFgUOERFJiwKHiIikRYFDRETS8v8BkS8/2jSgk8MAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Exponential rate, cubic fragmentation function" + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 3, + "source": [ + "from PySDM.initialisation.spectra import Exponential\n", + "from PySDM.physics.spontaneous_breakup_rates import Expon\n", + "from PySDM.physics.spontaneous_breakup_fragmentations import SpontaneousExpoCube\n", + "from PySDM.physics.constants import si\n", + "from PySDM.physics.formulae import Formulae\n", + "#from pystrict import strict\n", + "\n", + "#@strict\n", + "class Settings:\n", + "\n", + " def __init__(self):\n", + " self.formulae = Formulae()\n", + " self.n_sd = 12\n", + " self.n_part = 100 / si.cm**3\n", + " self.X0 = self.formulae.trivia.volume(radius=30.531 * si.micrometres)\n", + " self.dv = 1 * si.m**3\n", + " self.norm_factor = self.n_part * self.dv\n", + " self.rho = 1000 * si.kilogram / si.metre**3\n", + " self.dt = 1 * si.seconds\n", + " self.adaptive = False\n", + " self.seed = 44\n", + " self._steps = [0, 1, 2]\n", + " self.breakup_rate = Expon(1.58e-4 / si.s, 0.773 * si.mm)\n", + " self.fragmentation = SpontaneousExpoCube(C = 0.034 / si.mm**3)\n", + " self.spectrum = Exponential(norm_factor=self.norm_factor, scale=self.X0)\n", + " self.radius_bins_edges = np.logspace(np.log10(10 * si.um), np.log10(100 * si.um), num=128, endpoint=True)\n", + " self.radius_range = [0 * si.um, 1e6 * si.um]\n", + "\n", + " @property\n", + " def output_steps(self):\n", + " return [int(step/self.dt) for step in self._steps]" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 4, + "source": [ + "settings = Settings()\n", + "backend = CPU\n", + "\n", + "builder = Builder(n_sd=settings.n_sd, backend=backend, formulae=settings.formulae)\n", + "builder.set_environment(Box(dv=settings.dv, dt=settings.dt))\n", + "attributes = {}\n", + "attributes['volume'], attributes['n'] = ConstantMultiplicity(settings.spectrum).sample(settings.n_sd)\n", + "breakup = SpontaneousBreakup(settings.breakup_rate, settings.fragmentation, adaptive=settings.adaptive)\n", + "builder.add_dynamic(breakup)\n", + "\n", + "products = [ParticlesVolumeSpectrum(), WallTime(), ParticleMeanRadius(), ParticlesConcentration(radius_range = settings.radius_range)]\n", + "core = builder.build(attributes, products)" + ], + "outputs": [], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 7, + "source": [ + "for step in settings.output_steps:\n", + " #print(step, flush=True)\n", + " core.run(step - core.n_steps)\n", + " pyplot.step(x=settings.radius_bins_edges[:-1] / si.micrometres, \n", + " y=core.products['dv/dlnr'].get(settings.radius_bins_edges) * settings.rho,\n", + " where='post', label=\"t = {step*settings.dt}s\")\n", + " print(core.products['radius_m1'].get(), core.products['n_a_cm3'].get())\n", + " \n", + "pyplot.xscale(\"log\")\n", + "pyplot.xlabel(\"radius (um)\")\n", + "pyplot.ylabel(\"dm/dlnr\")\n", + "pyplot.legend()" + ], + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": [ + "[27.26101585] [99.997952]\n", + " 0.000158\n", + "[1.00296041 1.00427038 1.00506727 1.00567359 1.00617463 1.00660734\n", + " 1.00699151 1.00733912 1.00765802 1.00795369 1.0082301 1.00849026\n", + " 1.00873648 1.0089706 1.00919409 1.00940817 1.00961384 1.00981196\n", + " 1.01000324 1.01018831 1.01036769 1.01054186 1.01071122 1.01087614\n", + " 1.01103693 1.01119388 1.01134726 1.01149728 1.01164416 1.0117881\n", + " 1.01192925 1.01206778 1.01220384 1.01233756 1.01246906 1.01259845\n", + " 1.01272584 1.01285132 1.01297499 1.01309693 1.01321722 1.01333593\n", + " 1.01345313 1.01356888 1.01368325 1.0137963 1.01390807 1.01401862\n", + " 1.014128 1.01423625 1.01434341 1.01444953 1.01455464 1.01465878\n", + " 1.014762 1.01486431 1.01496575 1.01506636 1.01516615 1.01526517\n", + " 1.01536343 1.01546097 1.0155578 1.01565395 1.01574945 1.01584431\n", + " 1.01593855 1.01603221 1.01612528 1.0162178 1.01630978 1.01640124\n", + " 1.0164922 1.01658267 1.01667267 1.01676221 1.01685131 1.01693998\n", + " 1.01702824 1.01711611 1.01720358 1.01729068 1.01737742 1.01746381\n", + " 1.01754987 1.0176356 1.01772101 1.01780613 1.01789095 1.01797548\n", + " 1.01805975 1.01814376 1.01822751 1.01831102 1.0183943 1.01847735\n", + " 1.01856019 1.01864283 1.01872526 1.01880751 1.01888958 1.01897148\n", + " 1.01905322 1.0191348 1.01921623 1.01929753 1.01937869 1.01945973\n", + " 1.01954066 1.01962148 1.01970219 1.01978281 1.01986335 1.01994381\n", + " 1.02002419 1.02010452 1.02018478 1.020265 1.02034517 1.0204253\n", + " 1.02050541 1.02058549 1.02066556 1.02074562 1.02082568 1.02090575\n", + " 1.02098583 1.02106593 1.02114605 1.02122621 1.02130641 1.02138666\n", + " 1.02146696 1.02154732 1.02162776 1.02170827 1.02178886 1.02186954\n", + " 1.02195032 1.02203121 1.02211221 1.02219333 1.02227458 1.02235596\n", + " 1.02243749 1.02251917 1.02260101 1.02268301 1.02276519 1.02284756\n", + " 1.02293011 1.02301287 1.02309584 1.02317903 1.02326244 1.02334609\n", + " 1.02342999 1.02351414 1.02359856 1.02368325 1.02376823 1.0238535\n", + " 1.02393908 1.02402498 1.02411121 1.02419777 1.02428469 1.02437197\n", + " 1.02445963 1.02454767 1.02463612 1.02472497 1.02481426 1.02490399\n", + " 1.02499418 1.02508483 1.02517598 1.02526762 1.02535978 1.02545248\n", + " 1.02554574 1.02563956 1.02573397 1.025829 1.02592465 1.02602095\n", + " 1.02611792 1.02621559 1.02631397 1.02641309 1.02651298 1.02661366\n", + " 1.02671516 1.02681751 1.02692074 1.02702487 1.02712995 1.027236\n", + " 1.02734306 1.02745116 1.02756035 1.02767067 1.02778216 1.02789487\n", + " 1.02800883 1.02812411 1.02824076 1.02835883 1.02847838 1.02859948\n", + " 1.02872219 1.02884657 1.02897272 1.02910071 1.02923062 1.02936254\n", + " 1.02949658 1.02963284 1.02977143 1.02991247 1.0300561 1.03020246\n", + " 1.03035169 1.03050397 1.03065949 1.03081842 1.03098101 1.03114748\n", + " 1.03131809 1.03149314 1.03167295 1.03185788 1.03204832 1.03224474\n", + " 1.03244764 1.0326576 1.03287529 1.03310146 1.033337 1.03358293\n", + " 1.03384049 1.0341111 1.03439653 1.03469889 1.03502083 1.03536566\n", + " 1.03573767 1.0361425 1.03658782 1.03708446 1.03764843 1.038305\n", + " 1.03909776 1.04011279 1.04156481 1.04439057] 0.000158\n", + "[0.00015847 0.00015867 0.0001588 0.0001589 0.00015898 0.00015904\n", + " 0.0001591 0.00015916 0.00015921 0.00015926 0.0001593 0.00015934\n", + " 0.00015938 0.00015942 0.00015945 0.00015949 0.00015952 0.00015955\n", + " 0.00015958 0.00015961 0.00015964 0.00015967 0.00015969 0.00015972\n", + " 0.00015974 0.00015977 0.00015979 0.00015982 0.00015984 0.00015986\n", + " 0.00015988 0.00015991 0.00015993 0.00015995 0.00015997 0.00015999\n", + " 0.00016001 0.00016003 0.00016005 0.00016007 0.00016009 0.00016011\n", + " 0.00016013 0.00016014 0.00016016 0.00016018 0.0001602 0.00016021\n", + " 0.00016023 0.00016025 0.00016027 0.00016028 0.0001603 0.00016032\n", + " 0.00016033 0.00016035 0.00016036 0.00016038 0.0001604 0.00016041\n", + " 0.00016043 0.00016044 0.00016046 0.00016047 0.00016049 0.0001605\n", + " 0.00016052 0.00016053 0.00016055 0.00016056 0.00016058 0.00016059\n", + " 0.00016061 0.00016062 0.00016063 0.00016065 0.00016066 0.00016068\n", + " 0.00016069 0.0001607 0.00016072 0.00016073 0.00016075 0.00016076\n", + " 0.00016077 0.00016079 0.0001608 0.00016081 0.00016083 0.00016084\n", + " 0.00016085 0.00016087 0.00016088 0.00016089 0.00016091 0.00016092\n", + " 0.00016093 0.00016095 0.00016096 0.00016097 0.00016098 0.000161\n", + " 0.00016101 0.00016102 0.00016104 0.00016105 0.00016106 0.00016107\n", + " 0.00016109 0.0001611 0.00016111 0.00016113 0.00016114 0.00016115\n", + " 0.00016116 0.00016118 0.00016119 0.0001612 0.00016121 0.00016123\n", + " 0.00016124 0.00016125 0.00016127 0.00016128 0.00016129 0.0001613\n", + " 0.00016132 0.00016133 0.00016134 0.00016135 0.00016137 0.00016138\n", + " 0.00016139 0.0001614 0.00016142 0.00016143 0.00016144 0.00016146\n", + " 0.00016147 0.00016148 0.00016149 0.00016151 0.00016152 0.00016153\n", + " 0.00016155 0.00016156 0.00016157 0.00016158 0.0001616 0.00016161\n", + " 0.00016162 0.00016164 0.00016165 0.00016166 0.00016168 0.00016169\n", + " 0.0001617 0.00016172 0.00016173 0.00016174 0.00016176 0.00016177\n", + " 0.00016178 0.0001618 0.00016181 0.00016182 0.00016184 0.00016185\n", + " 0.00016186 0.00016188 0.00016189 0.00016191 0.00016192 0.00016193\n", + " 0.00016195 0.00016196 0.00016198 0.00016199 0.00016201 0.00016202\n", + " 0.00016204 0.00016205 0.00016207 0.00016208 0.0001621 0.00016211\n", + " 0.00016213 0.00016214 0.00016216 0.00016217 0.00016219 0.0001622\n", + " 0.00016222 0.00016224 0.00016225 0.00016227 0.00016229 0.0001623\n", + " 0.00016232 0.00016234 0.00016235 0.00016237 0.00016239 0.00016241\n", + " 0.00016243 0.00016244 0.00016246 0.00016248 0.0001625 0.00016252\n", + " 0.00016254 0.00016256 0.00016258 0.0001626 0.00016262 0.00016264\n", + " 0.00016266 0.00016268 0.0001627 0.00016273 0.00016275 0.00016277\n", + " 0.0001628 0.00016282 0.00016284 0.00016287 0.00016289 0.00016292\n", + " 0.00016295 0.00016298 0.000163 0.00016303 0.00016306 0.00016309\n", + " 0.00016313 0.00016316 0.00016319 0.00016323 0.00016327 0.00016331\n", + " 0.00016335 0.00016339 0.00016343 0.00016348 0.00016353 0.00016359\n", + " 0.00016365 0.00016371 0.00016378 0.00016386 0.00016395 0.00016405\n", + " 0.00016418 0.00016434 0.00016457 0.00016501] 1.0\n", + " 34000000.0\n" + ] + }, + { + "output_type": "error", + "ename": "TypingError", + "evalue": "Failed in nopython mode pipeline (step: nopython frontend)\nNo implementation of function Function() found for signature:\n \n >>> imul(array(int64, 1d, C), float64)\n \nThere are 8 candidate implementations:\n - Of which 4 did not match due to:\n Overload of function 'imul': File: : Line N/A.\n With argument(s): '(array(int64, 1d, C), float64)':\n No match.\n - Of which 2 did not match due to:\n Overload in function 'NumpyRulesInplaceArrayOperator.generic': File: numba/core/typing/npydecl.py: Line 213.\n With argument(s): '(array(int64, 1d, C), float64)':\n Rejected as the implementation raised a specific error:\n AttributeError: 'NoneType' object has no attribute 'args'\n raised from /Users/emilydejong/opt/anaconda3/envs/edjPySDM/lib/python3.9/site-packages/numba/core/typing/npydecl.py:224\n - Of which 2 did not match due to:\n Operator Overload in function 'imul': File: unknown: Line unknown.\n With argument(s): '(array(int64, 1d, C), float64)':\n No match for registered cases:\n * (int64, int64) -> int64\n * (int64, uint64) -> int64\n * (uint64, int64) -> int64\n * (uint64, uint64) -> uint64\n * (float32, float32) -> float32\n * (float64, float64) -> float64\n * (complex64, complex64) -> complex64\n * (complex128, complex128) -> complex128\n\nDuring: typing of intrinsic-call at /Users/emilydejong/Documents/PySDM/PySDM/backends/numba/impl/storage_impl.py (36)\n\nFile \"../../PySDM/backends/numba/impl/storage_impl.py\", line 36:\ndef multiply(output, multiplier):\n \n output *= multiplier\n\n ^\n", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypingError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mstep\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moutput_steps\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;31m#print(step, flush=True)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstep\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_steps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m pyplot.step(x=settings.radius_bins_edges[:-1] / si.micrometres, \n\u001b[1;32m 5\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mproducts\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'dv/dlnr'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mradius_bins_edges\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msettings\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrho\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/PySDM/PySDM/core.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, steps)\u001b[0m\n\u001b[1;32m 116\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdynamic\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdynamics\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 117\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimers\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 118\u001b[0;31m \u001b[0mdynamic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 119\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_steps\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 120\u001b[0m \u001b[0mreversed_order_so_that_environment_is_last\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreversed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mobservers\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/PySDM/PySDM/dynamics/spontaneous_breakup.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 100\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madaptive\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 101\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0m_\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__substeps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 102\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 103\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 104\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdt_left\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/PySDM/PySDM/dynamics/spontaneous_breakup.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 123\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 124\u001b[0m \u001b[0;31m# (3) Compute the number of fragments\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_n_fragment\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_fragment\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrand_frag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 126\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 127\u001b[0m \u001b[0;31m# (4) Perform the collisional-breakup step:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/PySDM/PySDM/dynamics/spontaneous_breakup.py\u001b[0m in \u001b[0;36mcompute_n_fragment\u001b[0;34m(self, n_fragment, u01)\u001b[0m\n\u001b[1;32m 139\u001b[0m \u001b[0;31m# (3) Compute n_fragment\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 140\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mcompute_n_fragment\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_fragment\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu01\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 141\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfragmentation\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_fragment\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mu01\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m~/Documents/PySDM/PySDM/physics/spontaneous_breakup_fragmentations/expo_cube.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, output, u01)\u001b[0m\n\u001b[1;32m 21\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparticles\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'volume'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 22\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 23\u001b[0;31m \u001b[0moutput\u001b[0m \u001b[0;34m*=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mC\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mregister\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbuilder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/Documents/PySDM/PySDM/backends/numba/storage.py\u001b[0m in \u001b[0;36m__imul__\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 71\u001b[0m \u001b[0mimpl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 72\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 73\u001b[0;31m \u001b[0mimpl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmultiply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 74\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 75\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/envs/edjPySDM/lib/python3.9/site-packages/numba/core/dispatcher.py\u001b[0m in \u001b[0;36m_compile_for_args\u001b[0;34m(self, *args, **kws)\u001b[0m\n\u001b[1;32m 418\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpatch_message\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 419\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 420\u001b[0;31m \u001b[0merror_rewrite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'typing'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 421\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mUnsupportedError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 422\u001b[0m \u001b[0;31m# Something unsupported is present in the user code, add help info\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m~/opt/anaconda3/envs/edjPySDM/lib/python3.9/site-packages/numba/core/dispatcher.py\u001b[0m in \u001b[0;36merror_rewrite\u001b[0;34m(e, issue_type)\u001b[0m\n\u001b[1;32m 359\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 360\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 361\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwith_traceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 362\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 363\u001b[0m \u001b[0margtypes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mTypingError\u001b[0m: Failed in nopython mode pipeline (step: nopython frontend)\nNo implementation of function Function() found for signature:\n \n >>> imul(array(int64, 1d, C), float64)\n \nThere are 8 candidate implementations:\n - Of which 4 did not match due to:\n Overload of function 'imul': File: : Line N/A.\n With argument(s): '(array(int64, 1d, C), float64)':\n No match.\n - Of which 2 did not match due to:\n Overload in function 'NumpyRulesInplaceArrayOperator.generic': File: numba/core/typing/npydecl.py: Line 213.\n With argument(s): '(array(int64, 1d, C), float64)':\n Rejected as the implementation raised a specific error:\n AttributeError: 'NoneType' object has no attribute 'args'\n raised from /Users/emilydejong/opt/anaconda3/envs/edjPySDM/lib/python3.9/site-packages/numba/core/typing/npydecl.py:224\n - Of which 2 did not match due to:\n Operator Overload in function 'imul': File: unknown: Line unknown.\n With argument(s): '(array(int64, 1d, C), float64)':\n No match for registered cases:\n * (int64, int64) -> int64\n * (int64, uint64) -> int64\n * (uint64, int64) -> int64\n * (uint64, uint64) -> uint64\n * (float32, float32) -> float32\n * (float64, float64) -> float64\n * (complex64, complex64) -> complex64\n * (complex128, complex128) -> complex128\n\nDuring: typing of intrinsic-call at /Users/emilydejong/Documents/PySDM/PySDM/backends/numba/impl/storage_impl.py (36)\n\nFile \"../../PySDM/backends/numba/impl/storage_impl.py\", line 36:\ndef multiply(output, multiplier):\n \n output *= multiplier\n\n ^\n" + ] + }, + { + "output_type": "display_data", + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUL0lEQVR4nO3df4xl5X3f8fenu4HYrhdngVjrXbY7FZtGg62QeMRiNarSuLhEtbNExd0NTo0iGizFKGnVqsKRjBOUVEWqmtoFWd4YYkDdgEXiMm1JaBwcJa7wZIcYFTMO6hQw7HYTlh9hE0eYLv32j3vu5u7dO8ydnR/33jnvlzSae859ztlzj87ezzw/znNSVUiS2ulvjPoAJEmjYwhIUosZApLUYoaAJLWYISBJLbZ11AewEhdddFHt2bNn1IchSRPlsccee7GqLh703kSFwJ49e5ifnx/1YUjSREnyraXeszlIklrMEJCkFjMEJKnFDAFJajFDQJJazBCQpBYzBCSpxSbqPgFtLofnnuPBx4+dXt5/+U6u27d7hEcktY81AY3Mg48fY+H4SQAWjp88IxAkbQxDQCM1vWMb93/sfUzv2DbqQ5FayRCQpBYzBCSpxQwBSWoxQ0CSWswQkKQWMwQkqcUMAUlqMUNAklrMEJCkFjMEJKnFnEBOE6F3sjknmpPWjiGgidA72RxgCEhrxBDQxHCSOWnt2ScgSS1mCEhSixkCktRihoAktZgdw1o1h29Kk8sQ0KoNM3xz0EPlJY2eIaA1sdzwzW5QTO/YdkZgSBot+wS0YXyovDR+DAFJarGhQiDJ1UmeSrKY5OYB75+f5P7m/bkke5r1VyV5LMkTze8f7dnmvc36xSSfSZI1+1SSpKEsGwJJtgB3AD8GTAM/mWS6r9gNwCtVdSnwq8BtzfoXgQ9V1XuA64F7e7b5LPAzwN7m5+pVfA5J0jkYpiZwBbBYVU9X1evAfcD+vjL7gbub1w8A70+Sqvp6Vf2fZv2TwFuaWsMOYFtVfa2qCrgHuGa1H0aStDLDhMBO4Pme5aPNuoFlquoU8CpwYV+Zfwz8cVV9pyl/dJl9ApDkxiTzSeZPnDgxxOFKkoa1IR3DSS6j00T0sZVuW1WHqmqmqmYuvvjitT84SWqxYULgGHBJz/KuZt3AMkm2AhcALzXLu4AvAR+tqv/dU37XMvvUJnB47jkOfO5R7w2QxtQwIXAE2JtkKsl5wEFgtq/MLJ2OX4BrgUeqqpK8A/hvwM1V9T+6havqOHAyyZXNqKCPAg+u7qNoHPXeJOZdwtL4WfaO4ao6leQm4GFgC3BXVT2Z5FZgvqpmgTuBe5MsAi/TCQqAm4BLgVuS3NKs+0BVvQD8LPAF4C3Abzc/2oS6N4kNY9D0Es5FJK2foaaNqKqHgIf61t3S8/o14MMDtvtl4JeX2Oc88O6VHKw2v0HTSxgC0vrxjmGNHaeXkDaOISBJLWYISFKLGQKS1GI+T0Ar1v8kMUmTyxDQivU/SWyQblB0R/pIGk+GgM7JSp4kZm1BGl+GgNbNSm4SA07XLqw5SBvHENBY6K0tWHOQNo4hoLFw3b7d3hksjYBDRCWpxQwBSWoxm4O06TgTqTQ8awIa2jAPiFk4fnLkD5HpvY9h4fjJMwJB0pkMAQ1tubH/+y/feXp456jvD3AmUmk4NgdpRXrH/vf/he0IH2nyWBOQpBYzBCSpxQwBSWoxQ0CSWswQkKQWc3SQluWzAaTNy5qAluWzAaTNy5qAhrLSZwNImgzWBDQSC8dPjnRqCUkd1gS04XyAjDQ+DAFtuLWaXsLZQqXVMwQ0sXo7rLtNS+MWAgaVxp19AlrSMFNHj9q4zxbqtNYad4aAljTM0FA7eJc37kGldrM5SG/qzYaG2sErTT5DQOfM5wdIk8/mIElqsaFCIMnVSZ5Kspjk5gHvn5/k/ub9uSR7mvUXJvlKkr9McnvfNr/f7PPx5ud71+QTadMbh+cYS5vFss1BSbYAdwBXAUeBI0lmq2qhp9gNwCtVdWmSg8BtwAHgNeCTwLubn34fqar5VX4GtUhv34NzGUmrN0yfwBXAYlU9DZDkPmA/0BsC+4FfbF4/ANyeJFX1beCrSS5du0PWehvnWUPth5DW1jDNQTuB53uWjzbrBpapqlPAq8CFQ+z715umoE8myRDltQGcNVRqj1GODvpIVR1L8nbgN4F/CtzTXyjJjcCNALt3+xfgRnHWUKkdhqkJHAMu6Vne1awbWCbJVuAC4KU322lVHWt+/wVwmE6z06Byh6pqpqpmLr744iEOV5I0rGFC4AiwN8lUkvOAg8BsX5lZ4Prm9bXAI1VVS+0wydYkFzWvvwv4IPCNlR68JGl1lm0OqqpTSW4CHga2AHdV1ZNJbgXmq2oWuBO4N8ki8DKdoAAgybPANuC8JNcAHwC+BTzcBMAW4MvAr63lB5MkLW+oPoGqegh4qG/dLT2vXwM+vMS2e5bY7XuHO0RJ0npx2gidNs5DQyWtD6eN0GkODZXax5qAzuDQUKldrAlIUosZApLUYoaAxpozhkrryz4Bje2oIGcMldafIaCxHRXkjKHS+jMEBDgqSGorQ0AaQrfJrGv/5TutpWhTsGNYGkK3yQw6ndW9gSBNMkNAGlK3yWycOs+l1bI5qMXGdVSQpI1jTaDFxnVUkKSNY02g5RwVJLWbNQFJajFDQJJazBCQpBazT6CFHBUkqcuaQAs5KkhSlzWBlnJU0GBOD6G2sSYg9XB6CLWNISD1cXoItYkhIEktZp9AizgqSFI/awIt4qggSf2sCbSMo4Ik9bImIEktZk1Ardd7b4D9JWobawItcHjuOQ587tHT4991pt57A+wvUdtYE2gBO4SXZ1+J2soQaAm/5CQNYnOQJLWYNQG1Tv8kcXYGq82GqgkkuTrJU0kWk9w84P3zk9zfvD+XZE+z/sIkX0nyl0lu79vmvUmeaLb5TJKsySfSaXYID9bbEQx2Bqvdlq0JJNkC3AFcBRwFjiSZraqFnmI3AK9U1aVJDgK3AQeA14BPAu9ufnp9FvgZYA54CLga+O3VfRz1skN4afaRSB3D1ASuABar6umqeh24D9jfV2Y/cHfz+gHg/UlSVd+uqq/SCYPTkuwAtlXV16qqgHuAa1bxObSE7pedc+JLGmSYENgJPN+zfLRZN7BMVZ0CXgUuXGafR5fZJwBJbkwyn2T+xIkTQxyuJGlYY98xXFWHgEMAMzMzNeLDmRiH555j7pmX2Te1fdSHMha8K1gabJiawDHgkp7lXc26gWWSbAUuAF5aZp+7ltmnVqH7hWdfQId3BUuDDVMTOALsTTJF54v6IHBdX5lZ4HrgUeBa4JGmrX+gqjqe5GSSK+l0DH8U+I/ncPx6E/umttsX0MPOYOlsy4ZAVZ1KchPwMLAFuKuqnkxyKzBfVbPAncC9SRaBl+kEBQBJngW2AecluQb4QDOy6GeBLwBvoTMqyJFBWhe9NQBJZxqqT6CqHqIzjLN33S09r18DPrzEtnuWWD/P2cNGpTXV2+xjE5B0trHvGNbK+AjJM123b7dNYtKbcO6gTcYbxCSthDWBTcgOUEnDsiYgSS1mTWCTsC9A0rmwJrBJ2Bcg6VxYE9hE7AuQtFLWBLRpLBw/6fMTpBUyBDaB7mRxbbb/8p2n+0JsEpOGZ3PQJuBkcd4UJp0rawKbhJPFSToXhoAktZghIEktZghIUosZApLUYo4OmkC9z8tt84ggSatnCEyg3uflStJqGAITykniJK0F+wQkqcUMAUlqMUNAklrMEJCkFrNjeEL0Dgude+Zl9k1tH/ERSdoMDIEJ4bBQSevB5qAJMr1jm0NDJa0pQ2DCLRw/aQ1B0jmzOWiC9U4Z4fQRks6FITDBfJqWpNWyOUiSWswQkKQWMwQkqcUMAUlqMTuGJ0x3FJCjgSStBUNgwjgiSNJaGqo5KMnVSZ5Kspjk5gHvn5/k/ub9uSR7et77RLP+qST/sGf9s0meSPJ4kvk1+TSSpBVZtiaQZAtwB3AVcBQ4kmS2qhZ6it0AvFJVlyY5CNwGHEgyDRwELgPeBXw5yfdV1RvNdn+/ql5cw88jSVqBYWoCVwCLVfV0Vb0O3Afs7yuzH7i7ef0A8P4kadbfV1XfqapngMVmf1qBw3PPMffMy6M+DEmb0DAhsBN4vmf5aLNuYJmqOgW8Cly4zLYF/PckjyW5cal/PMmNSeaTzJ84cWKIw918ulNI2xksaa2NsmP4h6vqWJLvBX43yZ9U1R/0F6qqQ8AhgJmZmdrogxyl7jMEFo6fZN/UdjuEJa25YWoCx4BLepZ3NesGlkmyFbgAeOnNtq2q7u8XgC9hM9FZugEwvWObtQBJ62KYEDgC7E0yleQ8Oh29s31lZoHrm9fXAo9UVTXrDzajh6aAvcAfJXlbkrcDJHkb8AHgG6v/OJtHtx9gesc27v/Y+6wFSFoXyzYHVdWpJDcBDwNbgLuq6skktwLzVTUL3Ancm2QReJlOUNCU+yKwAJwCPl5VbyR5J/ClTt8xW4HDVfU76/D5Jpb9AJI2wlB9AlX1EPBQ37pbel6/Bnx4iW1/BfiVvnVPAz+w0oNtA/sBJG0k5w4aM/YDSNpIThsxhrr9AJK03qwJSFKLGQJj4vDccxz43KPeGSxpQxkCY+LBx48ZAJI2nCEgSS1mx/AIdYeDOgpI0qgYAiNkE5CkUbM5SJJazBCQpBYzBCSpxQyBMbNvarsdxZI2jB3DY2Tf1Hani5C0oawJjIGF4ydZOH5y1IchqYWsCYyB6R3bAJ8dIGnjGQIj0n1ymE1AkkbJENhg3buEuzeJ+de/pFEyBDbQ4bnn+IUvPQH89SggnxwmaZQMgQ3SGwD/5ife45e/pLHg6KANYABIGleGwAZ48PFjgAEgafzYHLRO+qeJ3je13QCQNHYMgTXWP/pHksaZzUFrzACQNEkMAUlqMUNgnTkvkKRxZp/AOnNeIEnjzBBYQ935gMA7giVNBkNgDfSPCPJ+AEmTwhA4R90vfsC//iVNLENghXr/6n/7d29lesc2v/wlTSxDYBndL/2/ev0N3nreljPuAZjesc1nAUiaaIbAMrz5S9Jm5n0CktRiQ9UEklwNfBrYAny+qv5t3/vnA/cA7wVeAg5U1bPNe58AbgDeAH6uqh4eZp8bqbeTt1/3Rq99U9tPNwl1fzv2X9KkWzYEkmwB7gCuAo4CR5LMVtVCT7EbgFeq6tIkB4HbgANJpoGDwGXAu4AvJ/m+Zpvl9rmmfurzc3x18cWBX+a9o3v6Te/YZqevpE1rmJrAFcBiVT0NkOQ+YD/Q+4W9H/jF5vUDwO1J0qy/r6q+AzyTZLHZH0Psc8380n95kq8uvggwsH3f0T1aysLxkxz43KMsHD95+u7v3vXDbH8u20m9pt+1jU996LJ12fcwIbATeL5n+Siwb6kyVXUqyavAhc36r/Vt221DWW6fACS5EbgRYPfu1X9Jv2fnBWc16fjlr0F6m/u6NcL+9cs51+2kjTL2o4Oq6hBwCGBmZqbOZR+f+tBl65ai2ryu27d74B8IS60/1/1JozTM6KBjwCU9y7uadQPLJNkKXECng3ipbYfZpyRpnQ0TAkeAvUmmkpxHp6N3tq/MLHB98/pa4JGqqmb9wSTnJ5kC9gJ/NOQ+JUnrbNnmoKaN/ybgYTrDOe+qqieT3ArMV9UscCdwb9Px+zKdL3Wacl+k0+F7Cvh4Vb0BMGifa//xJElvJp0/2CfDzMxMzc/Pj/owJGmiJHmsqmYGvecdw5LUYoaAJLWYISBJLWYISFKLTVTHcJITwLdGfRyrdBHw4qgPYsx4Ts7mOTmb5+Rsw56Tv1VVFw96Y6JCYDNIMr9UL31beU7O5jk5m+fkbGtxTmwOkqQWMwQkqcUMgY13aNQHMIY8J2fznJzNc3K2VZ8T+wQkqcWsCUhSixkCktRihsA6SnJJkq8kWUjyZJKfb9ZvT/K7Sf5X8/t7Rn2sGy3JliRfT/Jfm+WpJHNJFpPc30wx3hpJ3pHkgSR/kuSbSd7X9uskyb9o/t98I8lvJPnutl0nSe5K8kKSb/SsG3hdpOMzzbn5n0l+aJh/wxBYX6eAf1lV08CVwMeTTAM3A79XVXuB32uW2+bngW/2LN8G/GpVXQq8AtwwkqManU8Dv1NV3w/8AJ1z09rrJMlO4OeAmap6N50p5w/SvuvkC8DVfeuWui5+jM4zW/bSeSTvZ4f5BwyBdVRVx6vqj5vXf0HnP/ZOYD9wd1PsbuCakRzgiCTZBfwj4PPNcoAfBR5oirTqnCS5APh7dJ7LQVW9XlV/TsuvEzrPO3lL87TCtwLHadl1UlV/QOcZLb2Wui72A/dUx9eAdyTZsdy/YQhskCR7gB8E5oB3VtXx5q0/Bd45quMakf8A/Gvg/zXLFwJ/XlWnmuWjdMKyLaaAE8CvN01kn0/yNlp8nVTVMeDfAc/R+fJ/FXiMdl8nXUtdFzuB53vKDXV+DIENkORvAr8J/POqOtn7XvMYztaM003yQeCFqnps1McyRrYCPwR8tqp+EPg2fU0/LbxOvofOX7ZTwLuAt3F2s0jrrcV1YQissyTfRScA/lNV/Vaz+s+61bTm9wujOr4R+LvAjyd5FriPTvX+03Sqrt3Hne4Cjo3m8EbiKHC0quaa5QfohEKbr5N/ADxTVSeq6v8Cv0Xn2mnzddK11HVxDLikp9xQ58cQWEdNW/edwDer6t/3vDULXN+8vh54cKOPbVSq6hNVtauq9tDp6Hukqj4CfAW4tinWtnPyp8DzSf5Os+r9dJ7L3drrhE4z0JVJ3tr8P+qek9ZeJz2Wui5mgY82o4SuBF7taTZakncMr6MkPwz8IfAEf93+/Qt0+gW+COymMzX2P6mq/s6fTS/JjwD/qqo+mORv06kZbAe+DvxUVX1nhIe3oZJcTqej/DzgaeCn6fyR1trrJMkvAQfojLL7OvDP6LRxt+Y6SfIbwI/QmTL6z4BPAf+ZAddFE5a302k2+yvgp6tq2YeyGwKS1GI2B0lSixkCktRihoAktZghIEktZghIUosZApLUYoaAJLXY/wfYOr5kES7VVwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + } + } + ], + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": null, + "source": [], + "outputs": [], + "metadata": {} + } + ], + "metadata": { + "kernelspec": { + "name": "python3", + "display_name": "Python 3.9.4 64-bit ('edjPySDM': conda)" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.4" + }, + "interpreter": { + "hash": "a75ec040163759486af1b34912b938f3bf396ad0d146e1bd1595d18e46184c90" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file