PySDM_examples.sandbox.seeding.simulation

  1import numpy as np
  2
  3from PySDM_examples.sandbox.seeding.settings import Settings
  4
  5from PySDM import Particulator
  6from PySDM.backends import CPU
  7from PySDM.environments import Parcel
  8from PySDM.dynamics import Condensation, AmbientThermodynamics, Coalescence, Seeding
  9from PySDM.dynamics.collisions.collision_kernels import Geometric
 10from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity
 11from PySDM import products
 12from PySDM.physics import si
 13
 14
 15class Simulation:
 16    def __init__(self, settings: Settings):
 17        environment = Parcel(
 18            dt=settings.timestep,
 19            mass_of_dry_air=settings.mass_of_dry_air,
 20            w=settings.updraft,
 21            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 22            p0=settings.initial_total_pressure,
 23            T0=settings.initial_temperature,
 24            backend=CPU(formulae=settings.formulae),
 25        )
 26        r_dry, n_in_dv = ConstantMultiplicity(
 27            settings.initial_aerosol_dry_radii
 28        ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend)
 29        attributes = environment.init_attributes(
 30            n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry
 31        )
 32        self.particulator = Particulator(
 33            environment=environment,
 34            n_sd=settings.n_sd_seeding + settings.n_sd_initial,
 35            dynamics=[
 36                AmbientThermodynamics(),
 37                Condensation(),
 38            ]
 39            + (
 40                []
 41                if not settings.enable_collisions
 42                else [
 43                    Coalescence(collision_kernel=Geometric()),
 44                ]
 45            )
 46            + [
 47                Seeding(
 48                    **{
 49                        k: getattr(settings, k)
 50                        for k in (
 51                            "super_droplet_injection_rate",
 52                            "seeded_particle_multiplicity",
 53                            "seeded_particle_extensive_attributes",
 54                        )
 55                    }
 56                ),
 57            ],
 58            attributes={
 59                k: np.pad(
 60                    array=v,
 61                    pad_width=(0, settings.n_sd_seeding),
 62                    mode="constant",
 63                    constant_values=np.nan if k == "multiplicity" else 0,
 64                )
 65                for k, v in attributes.items()
 66            },
 67            products=(
 68                products.SuperDropletCountPerGridbox(name="sd_count"),
 69                products.Time(),
 70                products.WaterMixingRatio(
 71                    radius_range=(settings.rain_water_radius_threshold, np.inf),
 72                    name="rain water mixing ratio",
 73                ),
 74                products.EffectiveRadius(
 75                    name="r_eff",
 76                    unit="um",
 77                    radius_range=(0.5 * si.um, 25 * si.um),
 78                ),
 79                products.ParticleConcentration(
 80                    name="n_drop",
 81                    unit="cm^-3",
 82                    radius_range=(0.5 * si.um, 25 * si.um),
 83                ),
 84            ),
 85        )
 86        self.n_steps = int(settings.t_max // settings.timestep)
 87
 88    def run(self):
 89        output = {
 90            "attributes": {"water mass": []},
 91            "products": {key: [] for key in self.particulator.products},
 92        }
 93        for step in range(self.n_steps + 1):
 94            if step != 0:
 95                self.particulator.advance(n_steps=1)
 96            for key, attr in output["attributes"].items():
 97                data = self.particulator.attributes[key].to_ndarray(raw=True)
 98                data[data == 0] = np.nan
 99                attr.append(data)
100            for key, prod in output["products"].items():
101                value = self.particulator.products[key].get()
102                if not isinstance(value, float):
103                    (value,) = value
104                prod.append(float(value))
105        for out in ("attributes", "products"):
106            for key, val in output[out].items():
107                output[out][key] = np.array(val)
108        return output
class Simulation:
 16class Simulation:
 17    def __init__(self, settings: Settings):
 18        environment = Parcel(
 19            dt=settings.timestep,
 20            mass_of_dry_air=settings.mass_of_dry_air,
 21            w=settings.updraft,
 22            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
 23            p0=settings.initial_total_pressure,
 24            T0=settings.initial_temperature,
 25            backend=CPU(formulae=settings.formulae),
 26        )
 27        r_dry, n_in_dv = ConstantMultiplicity(
 28            settings.initial_aerosol_dry_radii
 29        ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend)
 30        attributes = environment.init_attributes(
 31            n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry
 32        )
 33        self.particulator = Particulator(
 34            environment=environment,
 35            n_sd=settings.n_sd_seeding + settings.n_sd_initial,
 36            dynamics=[
 37                AmbientThermodynamics(),
 38                Condensation(),
 39            ]
 40            + (
 41                []
 42                if not settings.enable_collisions
 43                else [
 44                    Coalescence(collision_kernel=Geometric()),
 45                ]
 46            )
 47            + [
 48                Seeding(
 49                    **{
 50                        k: getattr(settings, k)
 51                        for k in (
 52                            "super_droplet_injection_rate",
 53                            "seeded_particle_multiplicity",
 54                            "seeded_particle_extensive_attributes",
 55                        )
 56                    }
 57                ),
 58            ],
 59            attributes={
 60                k: np.pad(
 61                    array=v,
 62                    pad_width=(0, settings.n_sd_seeding),
 63                    mode="constant",
 64                    constant_values=np.nan if k == "multiplicity" else 0,
 65                )
 66                for k, v in attributes.items()
 67            },
 68            products=(
 69                products.SuperDropletCountPerGridbox(name="sd_count"),
 70                products.Time(),
 71                products.WaterMixingRatio(
 72                    radius_range=(settings.rain_water_radius_threshold, np.inf),
 73                    name="rain water mixing ratio",
 74                ),
 75                products.EffectiveRadius(
 76                    name="r_eff",
 77                    unit="um",
 78                    radius_range=(0.5 * si.um, 25 * si.um),
 79                ),
 80                products.ParticleConcentration(
 81                    name="n_drop",
 82                    unit="cm^-3",
 83                    radius_range=(0.5 * si.um, 25 * si.um),
 84                ),
 85            ),
 86        )
 87        self.n_steps = int(settings.t_max // settings.timestep)
 88
 89    def run(self):
 90        output = {
 91            "attributes": {"water mass": []},
 92            "products": {key: [] for key in self.particulator.products},
 93        }
 94        for step in range(self.n_steps + 1):
 95            if step != 0:
 96                self.particulator.advance(n_steps=1)
 97            for key, attr in output["attributes"].items():
 98                data = self.particulator.attributes[key].to_ndarray(raw=True)
 99                data[data == 0] = np.nan
100                attr.append(data)
101            for key, prod in output["products"].items():
102                value = self.particulator.products[key].get()
103                if not isinstance(value, float):
104                    (value,) = value
105                prod.append(float(value))
106        for out in ("attributes", "products"):
107            for key, val in output[out].items():
108                output[out][key] = np.array(val)
109        return output
Simulation(settings: PySDM_examples.sandbox.seeding.settings.Settings)
17    def __init__(self, settings: Settings):
18        environment = Parcel(
19            dt=settings.timestep,
20            mass_of_dry_air=settings.mass_of_dry_air,
21            w=settings.updraft,
22            initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio,
23            p0=settings.initial_total_pressure,
24            T0=settings.initial_temperature,
25            backend=CPU(formulae=settings.formulae),
26        )
27        r_dry, n_in_dv = ConstantMultiplicity(
28            settings.initial_aerosol_dry_radii
29        ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend)
30        attributes = environment.init_attributes(
31            n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry
32        )
33        self.particulator = Particulator(
34            environment=environment,
35            n_sd=settings.n_sd_seeding + settings.n_sd_initial,
36            dynamics=[
37                AmbientThermodynamics(),
38                Condensation(),
39            ]
40            + (
41                []
42                if not settings.enable_collisions
43                else [
44                    Coalescence(collision_kernel=Geometric()),
45                ]
46            )
47            + [
48                Seeding(
49                    **{
50                        k: getattr(settings, k)
51                        for k in (
52                            "super_droplet_injection_rate",
53                            "seeded_particle_multiplicity",
54                            "seeded_particle_extensive_attributes",
55                        )
56                    }
57                ),
58            ],
59            attributes={
60                k: np.pad(
61                    array=v,
62                    pad_width=(0, settings.n_sd_seeding),
63                    mode="constant",
64                    constant_values=np.nan if k == "multiplicity" else 0,
65                )
66                for k, v in attributes.items()
67            },
68            products=(
69                products.SuperDropletCountPerGridbox(name="sd_count"),
70                products.Time(),
71                products.WaterMixingRatio(
72                    radius_range=(settings.rain_water_radius_threshold, np.inf),
73                    name="rain water mixing ratio",
74                ),
75                products.EffectiveRadius(
76                    name="r_eff",
77                    unit="um",
78                    radius_range=(0.5 * si.um, 25 * si.um),
79                ),
80                products.ParticleConcentration(
81                    name="n_drop",
82                    unit="cm^-3",
83                    radius_range=(0.5 * si.um, 25 * si.um),
84                ),
85            ),
86        )
87        self.n_steps = int(settings.t_max // settings.timestep)
particulator
n_steps
def run(self):
 89    def run(self):
 90        output = {
 91            "attributes": {"water mass": []},
 92            "products": {key: [] for key in self.particulator.products},
 93        }
 94        for step in range(self.n_steps + 1):
 95            if step != 0:
 96                self.particulator.advance(n_steps=1)
 97            for key, attr in output["attributes"].items():
 98                data = self.particulator.attributes[key].to_ndarray(raw=True)
 99                data[data == 0] = np.nan
100                attr.append(data)
101            for key, prod in output["products"].items():
102                value = self.particulator.products[key].get()
103                if not isinstance(value, float):
104                    (value,) = value
105                prod.append(float(value))
106        for out in ("attributes", "products"):
107            for key, val in output[out].items():
108                output[out][key] = np.array(val)
109        return output