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