PySDM_examples.deJong_Mackay_et_al_2023.simulation_0D

  1from collections import namedtuple
  2
  3import numpy as np
  4
  5from PySDM.backends import CPU
  6from PySDM.particulator import Particulator
  7from PySDM.dynamics import Coalescence, Collision
  8from PySDM.environments import Box
  9from PySDM.formulae import Formulae
 10from PySDM.initialisation.sampling.spectral_sampling import (
 11    ConstantMultiplicity,
 12    Logarithmic,
 13)
 14from PySDM.physics import si
 15from PySDM.products.collision.collision_rates import (
 16    BreakupRatePerGridbox,
 17    CoalescenceRatePerGridbox,
 18    CollisionRateDeficitPerGridbox,
 19    CollisionRatePerGridbox,
 20)
 21from PySDM.products.size_spectral import (
 22    NumberSizeSpectrum,
 23    ParticleVolumeVersusRadiusLogarithmSpectrum,
 24)
 25
 26
 27def run_box_breakup(
 28    settings, steps=None, backend_class=CPU, sample_in_radius=False, return_nv=False
 29):
 30    environment = Box(
 31        dv=settings.dv, dt=settings.dt, backend=backend_class(settings.formulae)
 32    )
 33    environment["rhod"] = 1.0
 34    attributes = {}
 35    if sample_in_radius:
 36        diams, attributes["multiplicity"] = Logarithmic(
 37            settings.spectrum
 38        ).sample_deterministic(settings.n_sd)
 39        radii = diams / 2
 40        attributes["volume"] = Formulae().trivia.volume(radius=radii)
 41    else:
 42        attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(
 43            settings.spectrum
 44        ).sample_deterministic(settings.n_sd)
 45    products = (
 46        ParticleVolumeVersusRadiusLogarithmSpectrum(
 47            radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr"
 48        ),
 49        NumberSizeSpectrum(radius_bins_edges=settings.radius_bins_edges, name="N(v)"),
 50        CollisionRatePerGridbox(name="cr"),
 51        CollisionRateDeficitPerGridbox(name="crd"),
 52        CoalescenceRatePerGridbox(name="cor"),
 53        BreakupRatePerGridbox(name="br"),
 54    )
 55    core = Particulator(
 56        n_sd=settings.n_sd,
 57        environment=environment,
 58        dynamics=(
 59            Collision(
 60                collision_kernel=settings.kernel,
 61                coalescence_efficiency=settings.coal_eff,
 62                breakup_efficiency=settings.break_eff,
 63                fragmentation_function=settings.fragmentation,
 64                adaptive=settings.adaptive,
 65                warn_overflows=settings.warn_overflows,
 66            ),
 67        ),
 68        attributes=attributes,
 69        products=products,
 70    )
 71
 72    if steps is None:
 73        steps = settings.output_steps
 74    y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
 75    if return_nv:
 76        y2 = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
 77    else:
 78        y2 = None
 79
 80    rates = np.zeros((len(steps), 4))
 81    # run
 82    for i, step in enumerate(steps):
 83        core.run(step - core.n_steps)
 84        y[i] = core.products["dv/dlnr"].get()[0]
 85        if return_nv:
 86            (y2[i],) = core.products["N(v)"].get()
 87        (rates[i, 0],) = core.products["cr"].get()
 88        (rates[i, 1],) = core.products["crd"].get()
 89        (rates[i, 2],) = core.products["cor"].get()
 90        (rates[i, 3],) = core.products["br"].get()
 91
 92    x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0]
 93
 94    return namedtuple("_", ("x", "y", "y2", "rates"))(x=x, y=y, y2=y2, rates=rates)
 95
 96
 97def run_box_NObreakup(settings, steps=None, backend_class=CPU):
 98    environment = Box(
 99        dv=settings.dv, dt=settings.dt, backend=backend_class(settings.formulae)
100    )
101    environment["rhod"] = 1.0
102    attributes = {}
103    attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(
104        settings.spectrum
105    ).sample_deterministic(settings.n_sd)
106    products = (
107        ParticleVolumeVersusRadiusLogarithmSpectrum(
108            radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr"
109        ),
110        CollisionRatePerGridbox(name="cr"),
111        CollisionRateDeficitPerGridbox(name="crd"),
112        CoalescenceRatePerGridbox(name="cor"),
113    )
114    core = Particulator(
115        n_sd=settings.n_sd,
116        environment=environment,
117        dynamics=(
118            Coalescence(
119                collision_kernel=settings.kernel,
120                coalescence_efficiency=settings.coal_eff,
121                adaptive=settings.adaptive,
122            ),
123        ),
124        attributes=attributes,
125        products=products,
126    )
127
128    if steps is None:
129        steps = settings.output_steps
130    y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
131    rates = np.zeros((len(steps), 4))
132    # run
133    for i, step in enumerate(steps):
134        core.run(step - core.n_steps)
135        (y[i],) = core.products["dv/dlnr"].get()
136        (rates[i, 0],) = core.products["cr"].get()
137        (rates[i, 1],) = core.products["crd"].get()
138        (rates[i, 2],) = core.products["cor"].get()
139
140    x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0]
141
142    return (x, y, rates)