PySDM_examples.Grabowski_and_Pawlowska_2023.settings

  1import numpy as np
  2from pystrict import strict
  3
  4from PySDM import Formulae
  5from PySDM.dynamics import condensation
  6from PySDM.initialisation.spectra import Lognormal, Sum
  7from PySDM.physics import si
  8
  9condensation_tolerance = condensation.DEFAULTS.rtol_thd / 100
 10
 11
 12@strict
 13class Settings:
 14    def __init__(
 15        self,
 16        *,
 17        aerosol: str,
 18        vertical_velocity: float,
 19        dt: float,
 20        n_sd: int,
 21        initial_temperature: float = 283 * si.K,
 22        initial_pressure: float = 900 * si.mbar,
 23        initial_relative_humidity: float = 0.97,
 24        displacement: float = 1000 * si.m,
 25        mass_accommodation_coefficient: float = 0.3,
 26        rtol_thd: float = condensation_tolerance,
 27        rtol_x: float = condensation_tolerance,
 28    ):
 29        self.formulae = Formulae(constants={"MAC": mass_accommodation_coefficient})
 30        self.n_sd = n_sd
 31        self.aerosol_modes_by_kappa = {
 32            "pristine": {
 33                1.28: Sum(
 34                    (
 35                        Lognormal(
 36                            norm_factor=125 / si.cm**3, m_mode=11 * si.nm, s_geom=1.2
 37                        ),
 38                        Lognormal(
 39                            norm_factor=65 / si.cm**3, m_mode=60 * si.nm, s_geom=1.7
 40                        ),
 41                    )
 42                )
 43            },
 44            "polluted": {
 45                1.28: Sum(
 46                    (
 47                        Lognormal(
 48                            norm_factor=160 / si.cm**3, m_mode=29 * si.nm, s_geom=1.36
 49                        ),
 50                        Lognormal(
 51                            norm_factor=380 / si.cm**3, m_mode=71 * si.nm, s_geom=1.57
 52                        ),
 53                    )
 54                )
 55            },
 56        }[aerosol]
 57
 58        self.vertical_velocity = vertical_velocity
 59        self.initial_pressure = initial_pressure
 60        self.initial_temperature = initial_temperature
 61        self.initial_relative_humidity = initial_relative_humidity
 62        self.t_max = displacement / vertical_velocity
 63        self.timestep = dt
 64        self.output_interval = self.timestep
 65        self.rtol_thd = rtol_thd
 66        self.rtol_x = rtol_x
 67
 68        self.initial_vapour_mixing_ratio = (
 69            self.formulae.trivia.water_vapour_mixing_ratio(
 70                self.initial_pressure,
 71                self.initial_relative_humidity,
 72                self.formulae.saturation_vapour_pressure.pvs_water(
 73                    self.initial_temperature
 74                ),
 75            )
 76        )
 77
 78    @property
 79    def initial_air_density(self):
 80        const = self.formulae.constants
 81        dry_air_density = (
 82            self.formulae.trivia.p_d(
 83                self.initial_pressure, self.initial_vapour_mixing_ratio
 84            )
 85            / self.initial_temperature
 86            / const.Rd
 87        )
 88        return dry_air_density * (1 + self.initial_vapour_mixing_ratio)
 89
 90    @property
 91    def nt(self) -> int:
 92        nt = self.t_max / self.timestep
 93        nt_int = round(nt)
 94        np.testing.assert_almost_equal(nt, nt_int)
 95        return nt_int
 96
 97    @property
 98    def steps_per_output_interval(self) -> int:
 99        return int(self.output_interval / self.timestep)
100
101    @property
102    def output_steps(self) -> np.ndarray:
103        return np.arange(0, self.nt + 1, self.steps_per_output_interval)