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