PySDM_examples.featured.Yang_et_al_2018.simulation

  1import numpy as np
  2
  3import PySDM.products as PySDM_products
  4from PySDM import Formulae, Particulator
  5from PySDM.backends import CPU
  6from PySDM.dynamics import AmbientThermodynamics, Condensation
  7from PySDM.environments import Parcel
  8from PySDM.physics import si
  9
 10
 11class Simulation:
 12    def __init__(self, settings, backend=CPU):
 13        dt_output = (
 14            settings.total_time / settings.n_steps
 15        )  # TODO #334 overwritten in notebook
 16        self.n_substeps = 1  # TODO #334 use condensation substeps
 17        while dt_output / self.n_substeps >= settings.dt_max:
 18            self.n_substeps += 1
 19        self.formulae = Formulae(
 20            diffusion_coordinate=settings.coord,
 21            saturation_vapour_pressure="AugustRocheMagnus",
 22        )
 23        self.bins_edges = self.formulae.trivia.volume(settings.r_bins_edges)
 24
 25        env = Parcel(
 26            dt=dt_output / self.n_substeps,
 27            mass_of_dry_air=settings.mass_of_dry_air,
 28            p0=settings.p0,
 29            initial_relative_humidity=settings.RH0,
 30            T0=settings.T0,
 31            w=settings.w,
 32            z0=settings.z0,
 33            backend=backend(formulae=self.formulae),
 34        )
 35
 36        condensation = Condensation(
 37            adaptive=settings.adaptive,
 38            rtol_x=settings.rtol_x,
 39            rtol_thd=settings.rtol_thd,
 40            dt_cond_range=settings.dt_cond_range,
 41        )
 42
 43        products = [
 44            PySDM_products.ParticleSizeSpectrumPerVolume(
 45                name="Particles Wet Size Spectrum",
 46                radius_bins_edges=settings.r_bins_edges,
 47            ),
 48            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
 49            PySDM_products.CondensationTimestepMax(name="dt_cond_max"),
 50            PySDM_products.RipeningRate(),
 51            PySDM_products.MeanRadius(
 52                name="r_mean_gt_1_um", radius_range=(1 * si.um, np.inf)
 53            ),
 54            PySDM_products.ActivatedMeanRadius(
 55                name="r_act", count_activated=True, count_unactivated=False
 56            ),
 57            PySDM_products.Time(name="t"),
 58        ]
 59
 60        attributes = env.init_attributes(
 61            n_in_dv=settings.n, kappa=settings.kappa, r_dry=settings.r_dry
 62        )
 63
 64        self.particulator = Particulator(
 65            attributes=attributes,
 66            products=products,
 67            n_sd=settings.n_sd,
 68            environment=env,
 69            dynamics=(
 70                AmbientThermodynamics(),
 71                condensation,
 72            ),
 73        )
 74
 75        self.n_steps = settings.n_steps
 76
 77    def save(self, output):
 78        _sp = self.particulator
 79        cell_id = 0
 80        output["r_bins_values"].append(
 81            (_sp.products["Particles Wet Size Spectrum"].get()).T
 82        )
 83        volume = _sp.attributes["volume"].to_ndarray()
 84        output["r"].append((self.formulae.trivia.radius(volume=volume)).T)
 85        output["S"].append(_sp.environment["RH"][cell_id] - 1)
 86        output["t"].append(_sp.products["t"].get())
 87        for key in ("water_vapour_mixing_ratio", "T", "z"):
 88            output[key].append(_sp.environment[key][cell_id])
 89        for key in (
 90            "dt_cond_max",
 91            "dt_cond_min",
 92            "ripening rate",
 93            "r_mean_gt_1_um",
 94            "r_act",
 95        ):
 96            output[key].append(_sp.products[key].get()[cell_id].copy())
 97
 98    def run(self):
 99        output = {
100            key: []
101            for key in (
102                "r",
103                "S",
104                "z",
105                "t",
106                "water_vapour_mixing_ratio",
107                "T",
108                "r_bins_values",
109                "dt_cond_max",
110                "dt_cond_min",
111                "ripening rate",
112                "r_mean_gt_1_um",
113                "r_act",
114            )
115        }
116
117        self.save(output)
118        for _ in range(self.n_steps):
119            self.particulator.advance(self.n_substeps)
120            self.save(output)
121        for k, v in output.items():
122            output[k] = np.asarray(v)
123        return output