PySDM_examples.featured.Arabas_and_Shima_2017.simulation

View Source

  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(formulae=settings.formulae),
 34        )
 35        r_wet = equilibrate_wet_radii(
 36            r_dry=r_dry,
 37            environment=environment,
 38            kappa_times_dry_volume=attributes["kappa times dry volume"],
 39        )
 40        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
 41        products = [
 42            PySDM_products.MeanRadius(name="radius_m1", unit="um"),
 43            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
 44            PySDM_products.ParcelDisplacement(name="z"),
 45            PySDM_products.AmbientRelativeHumidity(name="RH"),
 46            PySDM_products.PeakSaturation(name="S_max"),
 47            PySDM_products.Time(name="t"),
 48            PySDM_products.ActivatingRate(unit="s^-1 mg^-1", name="activating_rate"),
 49            PySDM_products.DeactivatingRate(
 50                unit="s^-1 mg^-1", name="deactivating_rate"
 51            ),
 52            PySDM_products.RipeningRate(unit="s^-1 mg^-1", name="ripening_rate"),
 53        ]
 54
 55        self.particulator = Particulator(
 56            n_sd=1,
 57            dynamics=[
 58                AmbientThermodynamics(),
 59                Condensation(
 60                    rtol_x=settings.rtol_x,
 61                    rtol_thd=settings.rtol_thd,
 62                    dt_cond_range=settings.dt_cond_range,
 63                ),
 64            ],
 65            attributes=attributes,
 66            products=products,
 67            environment=environment,
 68        )
 69
 70        self.n_output = settings.n_output
 71
 72    def save(self, output):
 73        cell_id = 0
 74        output["r"].append(
 75            self.particulator.products["radius_m1"].get(unit=const.si.m)[cell_id]
 76        )
 77        output["dt_cond_min"].append(
 78            self.particulator.products["dt_cond_min"].get()[cell_id]
 79        )
 80        output["z"].append(self.particulator.products["z"].get()[cell_id])
 81        output["RH"].append(self.particulator.products["RH"].get()[cell_id])
 82        output["t"].append(self.particulator.products["t"].get())
 83
 84        for event in ("activating", "deactivating", "ripening"):
 85            output[event + "_rate"].append(
 86                self.particulator.products[event + "_rate"].get()[cell_id]
 87            )
 88
 89    def run(self):
 90        output = {
 91            "r": [],
 92            "RH": [],
 93            "z": [],
 94            "t": [],
 95            "dt_cond_min": [],
 96            "activating_rate": [],
 97            "deactivating_rate": [],
 98            "ripening_rate": [],
 99        }
100
101        self.save(output)
102        for _ in range(self.n_output):
103            self.particulator.advance(self.n_substeps)
104            self.save(output)
105
106        return output

class Simulation: View Source

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

Simulation(settings, backend=<class 'PySDM.backends.Numba'>) View Source

14    def __init__(self, settings, backend=Numba):
15        t_half = settings.z_half / settings.w_avg
16
17        dt_output = (2 * t_half) / settings.n_output
18        self.n_substeps = 1
19        while dt_output / self.n_substeps >= settings.dt_max:  # TODO #334 dt_max
20            self.n_substeps += 1
21
22        attributes = {}
23        r_dry = np.array([settings.r_dry])
24        attributes["dry volume"] = settings.formulae.trivia.volume(radius=r_dry)
25        attributes["kappa times dry volume"] = attributes["dry volume"] * settings.kappa
26        attributes["multiplicity"] = np.array([settings.n_in_dv], dtype=np.int64)
27        environment = Parcel(
28            dt=dt_output / self.n_substeps,
29            mass_of_dry_air=settings.mass_of_dry_air,
30            p0=settings.p0,
31            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
32            T0=settings.T0,
33            w=settings.w,
34            backend=backend(formulae=settings.formulae),
35        )
36        r_wet = equilibrate_wet_radii(
37            r_dry=r_dry,
38            environment=environment,
39            kappa_times_dry_volume=attributes["kappa times dry volume"],
40        )
41        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
42        products = [
43            PySDM_products.MeanRadius(name="radius_m1", unit="um"),
44            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
45            PySDM_products.ParcelDisplacement(name="z"),
46            PySDM_products.AmbientRelativeHumidity(name="RH"),
47            PySDM_products.PeakSaturation(name="S_max"),
48            PySDM_products.Time(name="t"),
49            PySDM_products.ActivatingRate(unit="s^-1 mg^-1", name="activating_rate"),
50            PySDM_products.DeactivatingRate(
51                unit="s^-1 mg^-1", name="deactivating_rate"
52            ),
53            PySDM_products.RipeningRate(unit="s^-1 mg^-1", name="ripening_rate"),
54        ]
55
56        self.particulator = Particulator(
57            n_sd=1,
58            dynamics=[
59                AmbientThermodynamics(),
60                Condensation(
61                    rtol_x=settings.rtol_x,
62                    rtol_thd=settings.rtol_thd,
63                    dt_cond_range=settings.dt_cond_range,
64                ),
65            ],
66            attributes=attributes,
67            products=products,
68            environment=environment,
69        )
70
71        self.n_output = settings.n_output

n_substeps

particulator

n_output

def save(self, output): View Source

73    def save(self, output):
74        cell_id = 0
75        output["r"].append(
76            self.particulator.products["radius_m1"].get(unit=const.si.m)[cell_id]
77        )
78        output["dt_cond_min"].append(
79            self.particulator.products["dt_cond_min"].get()[cell_id]
80        )
81        output["z"].append(self.particulator.products["z"].get()[cell_id])
82        output["RH"].append(self.particulator.products["RH"].get()[cell_id])
83        output["t"].append(self.particulator.products["t"].get())
84
85        for event in ("activating", "deactivating", "ripening"):
86            output[event + "_rate"].append(
87                self.particulator.products[event + "_rate"].get()[cell_id]
88            )

def run(self): View Source

 90    def run(self):
 91        output = {
 92            "r": [],
 93            "RH": [],
 94            "z": [],
 95            "t": [],
 96            "dt_cond_min": [],
 97            "activating_rate": [],
 98            "deactivating_rate": [],
 99            "ripening_rate": [],
100        }
101
102        self.save(output)
103        for _ in range(self.n_output):
104            self.particulator.advance(self.n_substeps)
105            self.save(output)
106
107        return output