PySDM_examples.Arabas_and_Shima_2017.simulation

  1import numpy as np
  2
  3import PySDM.products as PySDM_products
  4from PySDM import Particulator
  5from PySDM.backends import Numba
  6from PySDM.dynamics import AmbientThermodynamics, Condensation
  7from PySDM.environments import Parcel
  8from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii
  9from PySDM.physics import constants as const
 10
 11
 12class Simulation:
 13    def __init__(self, settings, backend=Numba):
 14        t_half = settings.z_half / settings.w_avg
 15
 16        dt_output = (2 * t_half) / settings.n_output
 17        self.n_substeps = 1
 18        while dt_output / self.n_substeps >= settings.dt_max:  # TODO #334 dt_max
 19            self.n_substeps += 1
 20
 21        attributes = {}
 22        r_dry = np.array([settings.r_dry])
 23        attributes["dry volume"] = settings.formulae.trivia.volume(radius=r_dry)
 24        attributes["kappa times dry volume"] = attributes["dry volume"] * settings.kappa
 25        attributes["multiplicity"] = np.array([settings.n_in_dv], dtype=np.int64)
 26        environment = Parcel(
 27            dt=dt_output / self.n_substeps,
 28            mass_of_dry_air=settings.mass_of_dry_air,
 29            p0=settings.p0,
 30            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 31            T0=settings.T0,
 32            w=settings.w,
 33            backend=backend(
 34                formulae=settings.formulae,
 35                **(
 36                    {"override_jit_flags": {"parallel": False}}
 37                    if backend is Numba
 38                    else {}
 39                ),
 40            ),
 41        )
 42        r_wet = equilibrate_wet_radii(
 43            r_dry=r_dry,
 44            environment=environment,
 45            kappa_times_dry_volume=attributes["kappa times dry volume"],
 46        )
 47        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
 48        products = [
 49            PySDM_products.MeanRadius(name="radius_m1", unit="um"),
 50            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
 51            PySDM_products.ParcelDisplacement(name="z"),
 52            PySDM_products.AmbientRelativeHumidity(name="RH"),
 53            PySDM_products.PeakSaturation(name="S_max"),
 54            PySDM_products.Time(name="t"),
 55            PySDM_products.ActivatingRate(unit="s^-1 mg^-1", name="activating_rate"),
 56            PySDM_products.DeactivatingRate(
 57                unit="s^-1 mg^-1", name="deactivating_rate"
 58            ),
 59            PySDM_products.RipeningRate(unit="s^-1 mg^-1", name="ripening_rate"),
 60        ]
 61
 62        self.particulator = Particulator(
 63            n_sd=1,
 64            dynamics=[
 65                AmbientThermodynamics(),
 66                Condensation(
 67                    rtol_x=settings.rtol_x,
 68                    rtol_thd=settings.rtol_thd,
 69                    dt_cond_range=settings.dt_cond_range,
 70                ),
 71            ],
 72            attributes=attributes,
 73            products=products,
 74            environment=environment,
 75        )
 76
 77        self.n_output = settings.n_output
 78
 79    def save(self, output):
 80        cell_id = 0
 81        output["r"].append(
 82            self.particulator.products["radius_m1"].get(unit=const.si.m)[cell_id]
 83        )
 84        output["dt_cond_min"].append(
 85            self.particulator.products["dt_cond_min"].get()[cell_id]
 86        )
 87        output["z"].append(self.particulator.products["z"].get()[cell_id])
 88        output["RH"].append(self.particulator.products["RH"].get()[cell_id])
 89        output["t"].append(self.particulator.products["t"].get())
 90
 91        for event in ("activating", "deactivating", "ripening"):
 92            output[event + "_rate"].append(
 93                self.particulator.products[event + "_rate"].get()[cell_id]
 94            )
 95
 96    def run(self):
 97        output = {
 98            "r": [],
 99            "RH": [],
100            "z": [],
101            "t": [],
102            "dt_cond_min": [],
103            "activating_rate": [],
104            "deactivating_rate": [],
105            "ripening_rate": [],
106        }
107
108        self.save(output)
109        for _ in range(self.n_output):
110            self.particulator.run(self.n_substeps)
111            self.save(output)
112
113        return output