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        ]
 73
 74        attributes = environment.init_attributes(
 75            n_in_dv=settings.n, kappa=settings.kappa, r_dry=settings.r_dry
 76        )
 77
 78        self.particulator = builder.build(attributes, products)
 79
 80        self.n_steps = settings.n_steps
 81
 82    def save(self, output):
 83        _sp = self.particulator
 84        cell_id = 0
 85        output["r_bins_values"].append(
 86            _sp.products["Particles Wet Size Spectrum"].get()
 87        )
 88        volume = _sp.attributes["volume"].to_ndarray()
 89        output["r"].append(self.formulae.trivia.radius(volume=volume))
 90        output["S"].append(_sp.environment["RH"][cell_id] - 1)
 91        for key in ("water_vapour_mixing_ratio", "T", "z", "t"):
 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        return output