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