PySDM_examples.Yang_et_al_2018.simulation

  1import numpy as np
  2
  3import PySDM.products as PySDM_products
  4from PySDM import Builder, Formulae
  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        )
 34
 35        condensation = Condensation(
 36            adaptive=settings.adaptive,
 37            rtol_x=settings.rtol_x,
 38            rtol_thd=settings.rtol_thd,
 39            dt_cond_range=settings.dt_cond_range,
 40        )
 41
 42        builder = Builder(
 43            backend=backend(
 44                formulae=self.formulae, override_jit_flags={"parallel": False}
 45            ),
 46            n_sd=settings.n_sd,
 47            environment=env,
 48            dynamics=(
 49                AmbientThermodynamics(),
 50                condensation,
 51            ),
 52        )
 53
 54        environment = builder.particulator.environment
 55
 56        products = [
 57            PySDM_products.ParticleSizeSpectrumPerVolume(
 58                name="Particles Wet Size Spectrum",
 59                radius_bins_edges=settings.r_bins_edges,
 60            ),
 61            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
 62            PySDM_products.CondensationTimestepMax(name="dt_cond_max"),
 63            PySDM_products.RipeningRate(),
 64            PySDM_products.MeanRadius(
 65                name="r_mean_gt_1_um", radius_range=(1 * si.um, np.inf)
 66            ),
 67            PySDM_products.ActivatedMeanRadius(
 68                name="r_act", count_activated=True, count_unactivated=False
 69            ),
 70            PySDM_products.Time(name="t"),
 71        ]
 72
 73        attributes = environment.init_attributes(
 74            n_in_dv=settings.n, kappa=settings.kappa, r_dry=settings.r_dry
 75        )
 76
 77        self.particulator = builder.build(attributes, products)
 78
 79        self.n_steps = settings.n_steps
 80
 81    def save(self, output):
 82        _sp = self.particulator
 83        cell_id = 0
 84        output["r_bins_values"].append(
 85            (_sp.products["Particles Wet Size Spectrum"].get()).T
 86        )
 87        volume = _sp.attributes["volume"].to_ndarray()
 88        output["r"].append((self.formulae.trivia.radius(volume=volume)).T)
 89        output["S"].append(_sp.environment["RH"][cell_id] - 1)
 90        output["t"].append(_sp.products["t"].get())
 91        for key in ("water_vapour_mixing_ratio", "T", "z"):
 92            output[key].append(_sp.environment[key][cell_id])
 93        for key in (
 94            "dt_cond_max",
 95            "dt_cond_min",
 96            "ripening rate",
 97            "r_mean_gt_1_um",
 98            "r_act",
 99        ):
100            output[key].append(_sp.products[key].get()[cell_id].copy())
101
102    def run(self):
103        output = {
104            key: []
105            for key in (
106                "r",
107                "S",
108                "z",
109                "t",
110                "water_vapour_mixing_ratio",
111                "T",
112                "r_bins_values",
113                "dt_cond_max",
114                "dt_cond_min",
115                "ripening rate",
116                "r_mean_gt_1_um",
117                "r_act",
118            )
119        }
120
121        self.save(output)
122        for _ in range(self.n_steps):
123            self.particulator.run(self.n_substeps)
124            self.save(output)
125        for k, v in output.items():
126            output[k] = np.asarray(v)
127        return output