PySDM_examples.Grabowski_and_Pawlowska_2023.simulation

  1import numpy as np
  2from PySDM_examples.utils import BasicSimulation
  3
  4from PySDM import Builder
  5from PySDM.backends import CPU
  6from PySDM.backends.impl_numba.test_helpers import scipy_ode_condensation_solver
  7from PySDM.dynamics import AmbientThermodynamics, Condensation
  8from PySDM.environments import Parcel
  9from PySDM.initialisation import equilibrate_wet_radii
 10from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
 11from PySDM.physics import si
 12
 13
 14class Simulation(BasicSimulation):
 15    def __init__(
 16        self,
 17        settings,
 18        products=None,
 19        scipy_solver=False,
 20        sampling_class=ConstantMultiplicity,
 21    ):
 22        builder = Builder(
 23            n_sd=settings.n_sd,
 24            backend=CPU(
 25                formulae=settings.formulae, override_jit_flags={"parallel": False}
 26            ),
 27            environment=Parcel(
 28                dt=settings.timestep,
 29                p0=settings.initial_pressure,
 30                initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio,
 31                T0=settings.initial_temperature,
 32                w=settings.vertical_velocity,
 33                mass_of_dry_air=44 * si.kg,
 34            ),
 35        )
 36        builder.add_dynamic(AmbientThermodynamics())
 37        builder.add_dynamic(
 38            Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x)
 39        )
 40        for attribute in (
 41            "critical saturation",
 42            "equilibrium saturation",
 43            "critical volume",
 44        ):
 45            builder.request_attribute(attribute)
 46
 47        env = builder.particulator.environment
 48        volume = env.mass_of_dry_air / settings.initial_air_density
 49        attributes = {
 50            k: np.empty(0)
 51            for k in ("dry volume", "kappa times dry volume", "multiplicity")
 52        }
 53
 54        assert len(settings.aerosol_modes_by_kappa.keys()) == 1
 55        kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0]
 56        spectrum = settings.aerosol_modes_by_kappa[kappa]
 57
 58        r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd)
 59        v_dry = settings.formulae.trivia.volume(radius=r_dry)
 60        attributes["multiplicity"] = np.append(
 61            attributes["multiplicity"], n_per_volume * volume
 62        )
 63        attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
 64        attributes["kappa times dry volume"] = np.append(
 65            attributes["kappa times dry volume"], v_dry * kappa
 66        )
 67        r_wet = equilibrate_wet_radii(
 68            r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]),
 69            environment=env,
 70            kappa_times_dry_volume=attributes["kappa times dry volume"],
 71        )
 72        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
 73
 74        super().__init__(
 75            particulator=builder.build(attributes=attributes, products=products)
 76        )
 77        if scipy_solver:
 78            scipy_ode_condensation_solver.patch_particulator(self.particulator)
 79
 80        self.output_attributes = {
 81            "volume": tuple([] for _ in range(self.particulator.n_sd)),
 82            "dry volume": tuple([] for _ in range(self.particulator.n_sd)),
 83            "critical saturation": tuple([] for _ in range(self.particulator.n_sd)),
 84            "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)),
 85            "critical volume": tuple([] for _ in range(self.particulator.n_sd)),
 86            "multiplicity": tuple([] for _ in range(self.particulator.n_sd)),
 87        }
 88        self.settings = settings
 89
 90        self.__sanity_checks(attributes, volume)
 91
 92    def __sanity_checks(self, attributes, volume):
 93        for attribute in attributes.values():
 94            assert attribute.shape[0] == self.particulator.n_sd
 95        np.testing.assert_approx_equal(
 96            sum(attributes["multiplicity"]) / volume,
 97            sum(
 98                mode.norm_factor
 99                for mode in self.settings.aerosol_modes_by_kappa.values()
100            ),
101            significant=4,
102        )
103
104    def _save(self, output):
105        for key, attr in self.output_attributes.items():
106            attr_data = self.particulator.attributes[key].to_ndarray()
107            for drop_id in range(self.particulator.n_sd):
108                attr[drop_id].append(attr_data[drop_id])
109        super()._save(output)
110
111    def run(self):
112        output_products = super()._run(
113            self.settings.nt, self.settings.steps_per_output_interval
114        )
115        return {"products": output_products, "attributes": self.output_attributes}
 15class Simulation(BasicSimulation):
 16    def __init__(
 17        self,
 18        settings,
 19        products=None,
 20        scipy_solver=False,
 21        sampling_class=ConstantMultiplicity,
 22    ):
 23        builder = Builder(
 24            n_sd=settings.n_sd,
 25            backend=CPU(
 26                formulae=settings.formulae, override_jit_flags={"parallel": False}
 27            ),
 28            environment=Parcel(
 29                dt=settings.timestep,
 30                p0=settings.initial_pressure,
 31                initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio,
 32                T0=settings.initial_temperature,
 33                w=settings.vertical_velocity,
 34                mass_of_dry_air=44 * si.kg,
 35            ),
 36        )
 37        builder.add_dynamic(AmbientThermodynamics())
 38        builder.add_dynamic(
 39            Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x)
 40        )
 41        for attribute in (
 42            "critical saturation",
 43            "equilibrium saturation",
 44            "critical volume",
 45        ):
 46            builder.request_attribute(attribute)
 47
 48        env = builder.particulator.environment
 49        volume = env.mass_of_dry_air / settings.initial_air_density
 50        attributes = {
 51            k: np.empty(0)
 52            for k in ("dry volume", "kappa times dry volume", "multiplicity")
 53        }
 54
 55        assert len(settings.aerosol_modes_by_kappa.keys()) == 1
 56        kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0]
 57        spectrum = settings.aerosol_modes_by_kappa[kappa]
 58
 59        r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd)
 60        v_dry = settings.formulae.trivia.volume(radius=r_dry)
 61        attributes["multiplicity"] = np.append(
 62            attributes["multiplicity"], n_per_volume * volume
 63        )
 64        attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
 65        attributes["kappa times dry volume"] = np.append(
 66            attributes["kappa times dry volume"], v_dry * kappa
 67        )
 68        r_wet = equilibrate_wet_radii(
 69            r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]),
 70            environment=env,
 71            kappa_times_dry_volume=attributes["kappa times dry volume"],
 72        )
 73        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
 74
 75        super().__init__(
 76            particulator=builder.build(attributes=attributes, products=products)
 77        )
 78        if scipy_solver:
 79            scipy_ode_condensation_solver.patch_particulator(self.particulator)
 80
 81        self.output_attributes = {
 82            "volume": tuple([] for _ in range(self.particulator.n_sd)),
 83            "dry volume": tuple([] for _ in range(self.particulator.n_sd)),
 84            "critical saturation": tuple([] for _ in range(self.particulator.n_sd)),
 85            "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)),
 86            "critical volume": tuple([] for _ in range(self.particulator.n_sd)),
 87            "multiplicity": tuple([] for _ in range(self.particulator.n_sd)),
 88        }
 89        self.settings = settings
 90
 91        self.__sanity_checks(attributes, volume)
 92
 93    def __sanity_checks(self, attributes, volume):
 94        for attribute in attributes.values():
 95            assert attribute.shape[0] == self.particulator.n_sd
 96        np.testing.assert_approx_equal(
 97            sum(attributes["multiplicity"]) / volume,
 98            sum(
 99                mode.norm_factor
100                for mode in self.settings.aerosol_modes_by_kappa.values()
101            ),
102            significant=4,
103        )
104
105    def _save(self, output):
106        for key, attr in self.output_attributes.items():
107            attr_data = self.particulator.attributes[key].to_ndarray()
108            for drop_id in range(self.particulator.n_sd):
109                attr[drop_id].append(attr_data[drop_id])
110        super()._save(output)
111
112    def run(self):
113        output_products = super()._run(
114            self.settings.nt, self.settings.steps_per_output_interval
115        )
116        return {"products": output_products, "attributes": self.output_attributes}
Simulation( settings, products=None, scipy_solver=False, sampling_class=<class 'PySDM.initialisation.sampling.spectral_sampling.ConstantMultiplicity'>)
16    def __init__(
17        self,
18        settings,
19        products=None,
20        scipy_solver=False,
21        sampling_class=ConstantMultiplicity,
22    ):
23        builder = Builder(
24            n_sd=settings.n_sd,
25            backend=CPU(
26                formulae=settings.formulae, override_jit_flags={"parallel": False}
27            ),
28            environment=Parcel(
29                dt=settings.timestep,
30                p0=settings.initial_pressure,
31                initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio,
32                T0=settings.initial_temperature,
33                w=settings.vertical_velocity,
34                mass_of_dry_air=44 * si.kg,
35            ),
36        )
37        builder.add_dynamic(AmbientThermodynamics())
38        builder.add_dynamic(
39            Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x)
40        )
41        for attribute in (
42            "critical saturation",
43            "equilibrium saturation",
44            "critical volume",
45        ):
46            builder.request_attribute(attribute)
47
48        env = builder.particulator.environment
49        volume = env.mass_of_dry_air / settings.initial_air_density
50        attributes = {
51            k: np.empty(0)
52            for k in ("dry volume", "kappa times dry volume", "multiplicity")
53        }
54
55        assert len(settings.aerosol_modes_by_kappa.keys()) == 1
56        kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0]
57        spectrum = settings.aerosol_modes_by_kappa[kappa]
58
59        r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd)
60        v_dry = settings.formulae.trivia.volume(radius=r_dry)
61        attributes["multiplicity"] = np.append(
62            attributes["multiplicity"], n_per_volume * volume
63        )
64        attributes["dry volume"] = np.append(attributes["dry volume"], v_dry)
65        attributes["kappa times dry volume"] = np.append(
66            attributes["kappa times dry volume"], v_dry * kappa
67        )
68        r_wet = equilibrate_wet_radii(
69            r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]),
70            environment=env,
71            kappa_times_dry_volume=attributes["kappa times dry volume"],
72        )
73        attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet)
74
75        super().__init__(
76            particulator=builder.build(attributes=attributes, products=products)
77        )
78        if scipy_solver:
79            scipy_ode_condensation_solver.patch_particulator(self.particulator)
80
81        self.output_attributes = {
82            "volume": tuple([] for _ in range(self.particulator.n_sd)),
83            "dry volume": tuple([] for _ in range(self.particulator.n_sd)),
84            "critical saturation": tuple([] for _ in range(self.particulator.n_sd)),
85            "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)),
86            "critical volume": tuple([] for _ in range(self.particulator.n_sd)),
87            "multiplicity": tuple([] for _ in range(self.particulator.n_sd)),
88        }
89        self.settings = settings
90
91        self.__sanity_checks(attributes, volume)
output_attributes
settings
def run(self):
112    def run(self):
113        output_products = super()._run(
114            self.settings.nt, self.settings.steps_per_output_interval
115        )
116        return {"products": output_products, "attributes": self.output_attributes}