PySDM_examples.sandbox.Spichtinger_et_al_2023.simulation

  1import numpy as np
  2
  3from PySDM_examples.utils import BasicSimulation
  4
  5import PySDM.products as PySDM_products
  6from PySDM.backends import Numba
  7from PySDM import Particulator
  8from PySDM.dynamics import (
  9    AmbientThermodynamics,
 10    Condensation,
 11    Freezing,
 12    VapourDepositionOnIce,
 13)
 14from PySDM.environments import Parcel
 15from PySDM.initialisation import discretise_multiplicities
 16from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii
 17
 18
 19class Simulation(BasicSimulation):
 20    def __init__(self, settings, backend=Numba):
 21
 22        dt = settings.dt
 23
 24        formulae = settings.formulae
 25
 26        env = Parcel(
 27            mixed_phase=True,
 28            dt=dt,
 29            mass_of_dry_air=settings.mass_of_dry_air,
 30            p0=settings.initial_pressure,
 31            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 32            T0=settings.initial_temperature,
 33            w=settings.w_updraft,
 34            backend=backend(formulae=settings.formulae),
 35        )
 36
 37        self.n_sd = settings.n_sd
 38        self.multiplicities = discretise_multiplicities(
 39            settings.specific_concentration * env.mass_of_dry_air
 40        )
 41        self.r_dry = settings.r_dry
 42        v_dry = settings.formulae.trivia.volume(radius=self.r_dry)
 43        kappa = settings.kappa
 44
 45        self.r_wet = equilibrate_wet_radii(
 46            r_dry=self.r_dry,
 47            environment=env,
 48            kappa_times_dry_volume=kappa * v_dry,
 49        )
 50
 51        attributes = {
 52            "multiplicity": self.multiplicities,
 53            "dry volume": v_dry,
 54            "kappa times dry volume": kappa * v_dry,
 55            "signed water mass": formulae.particle_shape_and_density.radius_to_mass(
 56                self.r_wet
 57            ),
 58        }
 59
 60        products = [
 61            PySDM_products.Time(name="t"),
 62            PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"),
 63            PySDM_products.ParticleConcentration(
 64                name="n_i", unit="1/m**3", radius_range=(-np.inf, 0)
 65            ),
 66        ]
 67
 68        self.n_output = settings.n_output
 69        self.n_substeps = int(settings.t_duration / dt / self.n_output)
 70        super().__init__(
 71            Particulator(
 72                attributes=attributes,
 73                products=products,
 74                n_sd=settings.n_sd,
 75                environment=env,
 76                dynamics=(
 77                    AmbientThermodynamics(),
 78                    Condensation(),
 79                    VapourDepositionOnIce(),
 80                    Freezing(
 81                        homogeneous_freezing="time-dependent", immersion_freezing=None
 82                    ),
 83                ),
 84            )
 85        )
 86
 87    def save(self, output):
 88        cell_id = 0
 89        output["t"].append(self.particulator.products["t"].get())
 90        output["ni"].append(self.particulator.products["n_i"].get()[cell_id])
 91        output["RHi"].append(self.particulator.products["RH_ice"].get()[cell_id])
 92
 93    def run(self):
 94        output = {
 95            "t": [],
 96            "ni": [],
 97            "RHi": [],
 98        }
 99
100        self.save(output)
101
102        RHi_old = self.particulator.products["RH_ice"].get()[0].copy()
103        for _ in range(self.n_output):
104
105            self.particulator.advance(self.n_substeps)
106
107            self.save(output)
108
109            RHi = self.particulator.products["RH_ice"].get()[0].copy()
110            dRHi = (RHi_old - RHi) / RHi_old
111            if dRHi > 0.0 and RHi < 130.0:
112                print("break")
113                break
114            RHi_old = RHi
115
116        return output["ni"][-1]