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.builder import Builder
  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    builder = Builder(
 31        n_sd=settings.n_sd,
 32        backend=backend_class(settings.formulae),
 33        environment=Box(dv=settings.dv, dt=settings.dt),
 34        dynamics=(
 35            Collision(
 36                collision_kernel=settings.kernel,
 37                coalescence_efficiency=settings.coal_eff,
 38                breakup_efficiency=settings.break_eff,
 39                fragmentation_function=settings.fragmentation,
 40                adaptive=settings.adaptive,
 41                warn_overflows=settings.warn_overflows,
 42            ),
 43        ),
 44    )
 45    builder.particulator.environment["rhod"] = 1.0
 46    attributes = {}
 47    if sample_in_radius:
 48        diams, attributes["multiplicity"] = Logarithmic(
 49            settings.spectrum
 50        ).sample_deterministic(settings.n_sd)
 51        radii = diams / 2
 52        attributes["volume"] = Formulae().trivia.volume(radius=radii)
 53    else:
 54        attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(
 55            settings.spectrum
 56        ).sample_deterministic(settings.n_sd)
 57    products = (
 58        ParticleVolumeVersusRadiusLogarithmSpectrum(
 59            radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr"
 60        ),
 61        NumberSizeSpectrum(radius_bins_edges=settings.radius_bins_edges, name="N(v)"),
 62        CollisionRatePerGridbox(name="cr"),
 63        CollisionRateDeficitPerGridbox(name="crd"),
 64        CoalescenceRatePerGridbox(name="cor"),
 65        BreakupRatePerGridbox(name="br"),
 66    )
 67    core = builder.build(attributes, products)
 68
 69    if steps is None:
 70        steps = settings.output_steps
 71    y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
 72    if return_nv:
 73        y2 = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
 74    else:
 75        y2 = None
 76
 77    rates = np.zeros((len(steps), 4))
 78    # run
 79    for i, step in enumerate(steps):
 80        core.run(step - core.n_steps)
 81        y[i] = core.products["dv/dlnr"].get()[0]
 82        if return_nv:
 83            (y2[i],) = core.products["N(v)"].get()
 84        (rates[i, 0],) = core.products["cr"].get()
 85        (rates[i, 1],) = core.products["crd"].get()
 86        (rates[i, 2],) = core.products["cor"].get()
 87        (rates[i, 3],) = core.products["br"].get()
 88
 89    x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0]
 90
 91    return namedtuple("_", ("x", "y", "y2", "rates"))(x=x, y=y, y2=y2, rates=rates)
 92
 93
 94def run_box_NObreakup(settings, steps=None, backend_class=CPU):
 95    builder = Builder(
 96        n_sd=settings.n_sd,
 97        backend=backend_class(settings.formulae),
 98        environment=Box(dv=settings.dv, dt=settings.dt),
 99        dynamics=(
100            Coalescence(
101                collision_kernel=settings.kernel,
102                coalescence_efficiency=settings.coal_eff,
103                adaptive=settings.adaptive,
104            ),
105        ),
106    )
107    builder.particulator.environment["rhod"] = 1.0
108    attributes = {}
109    attributes["volume"], attributes["multiplicity"] = ConstantMultiplicity(
110        settings.spectrum
111    ).sample_deterministic(settings.n_sd)
112    products = (
113        ParticleVolumeVersusRadiusLogarithmSpectrum(
114            radius_bins_edges=settings.radius_bins_edges, name="dv/dlnr"
115        ),
116        CollisionRatePerGridbox(name="cr"),
117        CollisionRateDeficitPerGridbox(name="crd"),
118        CoalescenceRatePerGridbox(name="cor"),
119    )
120    core = builder.build(attributes, products)
121
122    if steps is None:
123        steps = settings.output_steps
124    y = np.ndarray((len(steps), len(settings.radius_bins_edges) - 1))
125    rates = np.zeros((len(steps), 4))
126    # run
127    for i, step in enumerate(steps):
128        core.run(step - core.n_steps)
129        (y[i],) = core.products["dv/dlnr"].get()
130        (rates[i, 0],) = core.products["cr"].get()
131        (rates[i, 1],) = core.products["crd"].get()
132        (rates[i, 2],) = core.products["cor"].get()
133
134    x = (settings.radius_bins_edges[:-1] / si.micrometres,)[0]
135
136    return (x, y, rates)