PySDM_examples.Arabas_and_Shima_2017.simulation

  1import numpy as np
  2
  3import PySDM.products as PySDM_products
  4from PySDM.backends import CPU
  5from PySDM.builder import Builder
  6from PySDM.dynamics import AmbientThermodynamics, Condensation
  7from PySDM.environments import Parcel
  8from PySDM.initialisation import equilibrate_wet_radii
  9from PySDM.physics import constants as const
 10
 11
 12class Simulation:
 13    def __init__(self, settings, backend=CPU):
 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        builder = Builder(
 22            backend=backend(
 23                formulae=settings.formulae,
 24                **(
 25                    {"override_jit_flags": {"parallel": False}}
 26                    if backend == CPU
 27                    else {}
 28                )
 29            ),
 30            n_sd=1,
 31            environment=Parcel(
 32                dt=dt_output / self.n_substeps,
 33                mass_of_dry_air=settings.mass_of_dry_air,
 34                p0=settings.p0,
 35                initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 36                T0=settings.T0,
 37                w=settings.w,
 38            ),
 39        )
 40
 41        builder.add_dynamic(AmbientThermodynamics())
 42        builder.add_dynamic(
 43            Condensation(
 44                rtol_x=settings.rtol_x,
 45                rtol_thd=settings.rtol_thd,
 46                dt_cond_range=settings.dt_cond_range,
 47            )
 48        )
 49        attributes = {}
 50        r_dry = np.array([settings.r_dry])
 51        attributes["dry volume"] = settings.formulae.trivia.volume(radius=r_dry)
 52        attributes["kappa times dry volume"] = attributes["dry volume"] * settings.kappa
 53        attributes["multiplicity"] = np.array([settings.n_in_dv], dtype=np.int64)
 54        environment = builder.particulator.environment
 55        r_wet = equilibrate_wet_radii(
 56            r_dry=r_dry,
 57            environment=environment,
 58            kappa_times_dry_volume=attributes["kappa times dry volume"],
 59        )
 60        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
 61        products = [
 62            PySDM_products.MeanRadius(name="radius_m1", unit="um"),
 63            PySDM_products.CondensationTimestepMin(name="dt_cond_min"),
 64            PySDM_products.ParcelDisplacement(name="z"),
 65            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
 66            PySDM_products.Time(name="t"),
 67            PySDM_products.ActivatingRate(unit="s^-1 mg^-1", name="activating_rate"),
 68            PySDM_products.DeactivatingRate(
 69                unit="s^-1 mg^-1", name="deactivating_rate"
 70            ),
 71            PySDM_products.RipeningRate(unit="s^-1 mg^-1", name="ripening_rate"),
 72            PySDM_products.PeakSupersaturation(unit="%", name="S_max"),
 73        ]
 74
 75        self.particulator = builder.build(attributes, products)
 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["S"].append(self.particulator.products["RH"].get()[cell_id] / 100 - 1)
 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            "S": [],
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