PySDM_examples.Morrison_and_Grabowski_2007.common

  1import numba
  2import numpy
  3import numpy as np  # pylint: disable=reimported
  4import PyMPDATA
  5import scipy
  6
  7import PySDM
  8from PySDM import Formulae
  9from PySDM.dynamics import collisions, condensation, displacement
 10from PySDM.dynamics.collisions.breakup_efficiencies import ConstEb
 11from PySDM.dynamics.collisions.breakup_fragmentations import AlwaysN
 12from PySDM.dynamics.collisions.coalescence_efficiencies import ConstEc
 13from PySDM.dynamics.collisions.collision_kernels import Geometric
 14from PySDM.initialisation import spectra
 15from PySDM.physics import si
 16
 17
 18class Common:
 19    def __init__(self, formulae: Formulae):
 20        self.formulae = formulae
 21        const = formulae.constants
 22
 23        self.condensation_rtol_x = condensation.DEFAULTS.rtol_x
 24        self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd
 25        self.condensation_adaptive = True
 26        self.condensation_substeps = -1
 27        self.condensation_dt_cond_range = condensation.DEFAULTS.cond_range
 28        self.condensation_schedule = condensation.DEFAULTS.schedule
 29
 30        self.coalescence_adaptive = True
 31        self.coalescence_dt_coal_range = collisions.collision.DEFAULTS.dt_coal_range
 32        self.coalescence_optimized_random = True
 33        self.coalescence_substeps = 1
 34        self.kernel = Geometric(collection_efficiency=1)
 35        self.coalescence_efficiency = ConstEc(Ec=1.0)
 36        self.breakup_efficiency = ConstEb(Eb=1.0)
 37        self.breakup_fragmentation = AlwaysN(n=2)
 38
 39        self.freezing_singular = True
 40        self.freezing_thaw = False
 41        self.freezing_inp_spec = None
 42
 43        self.displacement_adaptive = displacement.DEFAULTS.adaptive
 44        self.displacement_rtol = displacement.DEFAULTS.rtol
 45        self.freezing_inp_frac = 1
 46
 47        self.n_sd_per_gridbox = 20
 48
 49        self.aerosol_radius_threshold = 0.5 * si.micrometre
 50        self.drizzle_radius_threshold = 25 * si.micrometre
 51
 52        self.r_bins_edges = np.logspace(
 53            np.log10(0.001 * si.micrometre),
 54            np.log10(100 * si.micrometre),
 55            64,
 56            endpoint=True,
 57        )
 58        self.T_bins_edges = np.linspace(const.T0 - 40, const.T0 - 20, 64, endpoint=True)
 59
 60        # TODO #599
 61        n_bins_per_phase = 25
 62        solid_phase_radii = (
 63            np.linspace(-n_bins_per_phase, -1, n_bins_per_phase + 1) * si.um
 64        )
 65        liquid_phase_radii = (
 66            np.linspace(0, n_bins_per_phase, n_bins_per_phase + 1) * si.um
 67        )
 68        self.terminal_velocity_radius_bin_edges = np.concatenate(
 69            [solid_phase_radii, liquid_phase_radii]
 70        )
 71
 72        self.output_interval = 1 * si.minute
 73        self.spin_up_time = 0
 74
 75        self.mode_1 = spectra.Lognormal(
 76            norm_factor=60 / si.centimetre**3 / const.rho_STP,
 77            m_mode=0.04 * si.micrometre,
 78            s_geom=1.4,
 79        )
 80        self.mode_2 = spectra.Lognormal(
 81            norm_factor=40 / si.centimetre**3 / const.rho_STP,
 82            m_mode=0.15 * si.micrometre,
 83            s_geom=1.6,
 84        )
 85        self.spectrum_per_mass_of_dry_air = spectra.Sum((self.mode_1, self.mode_2))
 86        self.kappa = 1  # TODO #441!
 87
 88        self.processes = {
 89            "particle advection": True,
 90            "fluid advection": True,
 91            "coalescence": True,
 92            "condensation": True,
 93            "sedimentation": True,
 94            "breakup": False,
 95            "freezing": False,
 96        }
 97
 98        self.mpdata_iters = 2
 99        self.mpdata_iga = True
100        self.mpdata_fct = True
101        self.mpdata_tot = True
102
103        key_packages = [PySDM, PyMPDATA, numba, numpy, scipy]
104        try:
105            import ThrustRTC  # pylint: disable=import-outside-toplevel
106
107            key_packages.append(ThrustRTC)
108        except:  # pylint: disable=bare-except
109            pass
110        self.versions = {}
111        for pkg in key_packages:
112            try:
113                self.versions[pkg.__name__] = pkg.__version__
114            except AttributeError:
115                pass
116        self.versions = str(self.versions)
117
118        self.dt = None
119        self.simulation_time = None
120        self.grid = None
121        self.p0 = None
122        self.initial_water_vapour_mixing_ratio = None
123        self.th_std0 = None
124        self.size = None
125
126    @property
127    def n_steps(self) -> int:
128        return int(self.simulation_time / self.dt)  # TODO #413
129
130    @property
131    def steps_per_output_interval(self) -> int:
132        return int(self.output_interval / self.dt)
133
134    @property
135    def n_spin_up(self) -> int:
136        return int(self.spin_up_time / self.dt)
137
138    @property
139    def output_steps(self) -> np.ndarray:
140        return np.arange(0, self.n_steps + 1, self.steps_per_output_interval)
141
142    @property
143    def n_sd(self):
144        return self.grid[0] * self.grid[1] * self.n_sd_per_gridbox
145
146    @property
147    def initial_vapour_mixing_ratio_profile(self):
148        return np.full(self.grid[-1], self.initial_water_vapour_mixing_ratio)
149
150    @property
151    def initial_dry_potential_temperature_profile(self):
152        return np.full(
153            self.grid[-1],
154            self.formulae.state_variable_triplet.th_dry(
155                self.th_std0, self.initial_water_vapour_mixing_ratio
156            ),
157        )