PySDM_examples.Shipway_and_Hill_2012.simulation

View Source

  1from collections import namedtuple
  2
  3import numpy as np
  4from PySDM_examples.Shipway_and_Hill_2012.mpdata_1d import MPDATA_1D
  5
  6import PySDM.products as PySDM_products
  7from PySDM import Builder
  8from PySDM.backends import CPU
  9from PySDM.dynamics import (
 10    AmbientThermodynamics,
 11    Coalescence,
 12    Condensation,
 13    Displacement,
 14    EulerianAdvection,
 15)
 16from PySDM.environments.kinematic_1d import Kinematic1D
 17from PySDM.impl.mesh import Mesh
 18from PySDM.initialisation.sampling import spatial_sampling, spectral_sampling
 19
 20
 21class Simulation:
 22    def __init__(self, settings, backend=CPU):
 23        self.nt = settings.nt
 24        self.z0 = -settings.particle_reservoir_depth
 25        self.save_spec_and_attr_times = settings.save_spec_and_attr_times
 26        self.number_of_bins = settings.number_of_bins
 27
 28        self.particulator = None
 29        self.output_attributes = None
 30        self.output_products = None
 31
 32        self.mesh = Mesh(
 33            grid=(settings.nz,),
 34            size=(settings.z_max + settings.particle_reservoir_depth,),
 35        )
 36
 37        env = Kinematic1D(
 38            dt=settings.dt,
 39            mesh=self.mesh,
 40            thd_of_z=settings.thd,
 41            rhod_of_z=settings.rhod,
 42            z0=-settings.particle_reservoir_depth,
 43        )
 44
 45        def zZ_to_z_above_reservoir(zZ):
 46            z_above_reservoir = zZ * (settings.nz * settings.dz) + self.z0
 47            return z_above_reservoir
 48
 49        mpdata = MPDATA_1D(
 50            nz=settings.nz,
 51            dt=settings.dt,
 52            mpdata_settings=settings.mpdata_settings,
 53            advector_of_t=lambda t: settings.rho_times_w(t) * settings.dt / settings.dz,
 54            advectee_of_zZ_at_t0=lambda zZ: settings.water_vapour_mixing_ratio(
 55                zZ_to_z_above_reservoir(zZ)
 56            ),
 57            g_factor_of_zZ=lambda zZ: settings.rhod(zZ_to_z_above_reservoir(zZ)),
 58        )
 59
 60        _extra_nz = settings.particle_reservoir_depth // settings.dz
 61        _z_vec = settings.dz * np.linspace(
 62            -_extra_nz, settings.nz - _extra_nz, settings.nz + 1
 63        )
 64        self.g_factor_vec = settings.rhod(_z_vec)
 65
 66        self.builder = Builder(
 67            n_sd=settings.n_sd,
 68            backend=backend(formulae=settings.formulae),
 69            environment=env,
 70        )
 71        self.builder.add_dynamic(AmbientThermodynamics())
 72
 73        if settings.enable_condensation:
 74            self.builder.add_dynamic(
 75                Condensation(
 76                    adaptive=settings.condensation_adaptive,
 77                    rtol_thd=settings.condensation_rtol_thd,
 78                    rtol_x=settings.condensation_rtol_x,
 79                    update_thd=settings.condensation_update_thd,
 80                )
 81            )
 82        self.builder.add_dynamic(EulerianAdvection(mpdata))
 83
 84        self.products = []
 85        if settings.precip:
 86            self.add_collision_dynamic(self.builder, settings, self.products)
 87
 88        displacement = Displacement(
 89            enable_sedimentation=settings.precip,
 90            precipitation_counting_level_index=int(
 91                settings.particle_reservoir_depth / settings.dz
 92            ),
 93        )
 94        self.builder.add_dynamic(displacement)
 95        self.attributes = self.builder.particulator.environment.init_attributes(
 96            spatial_discretisation=spatial_sampling.Pseudorandom(),
 97            spectral_discretisation=spectral_sampling.ConstantMultiplicity(
 98                spectrum=settings.wet_radius_spectrum_per_mass_of_dry_air
 99            ),
100            kappa=settings.kappa,
101            collisions_only=not settings.enable_condensation,
102            z_part=settings.z_part,
103        )
104        self.products += [
105            PySDM_products.WaterMixingRatio(
106                name="cloud water mixing ratio",
107                unit="g/kg",
108                radius_range=settings.cloud_water_radius_range,
109            ),
110            PySDM_products.WaterMixingRatio(
111                name="rain water mixing ratio",
112                unit="g/kg",
113                radius_range=settings.rain_water_radius_range,
114            ),
115            PySDM_products.AmbientDryAirDensity(name="rhod"),
116            PySDM_products.AmbientDryAirPotentialTemperature(name="thd"),
117            PySDM_products.ParticleSizeSpectrumPerVolume(
118                name="wet spectrum", radius_bins_edges=settings.r_bins_edges
119            ),
120            PySDM_products.ParticleConcentration(
121                name="nc", radius_range=settings.cloud_water_radius_range
122            ),
123            PySDM_products.ParticleConcentration(
124                name="nr", radius_range=settings.rain_water_radius_range
125            ),
126            PySDM_products.ParticleConcentration(
127                name="na", radius_range=(0, settings.cloud_water_radius_range[0])
128            ),
129            PySDM_products.MeanRadius(),
130            PySDM_products.EffectiveRadius(
131                radius_range=settings.cloud_water_radius_range
132            ),
133            PySDM_products.SuperDropletCountPerGridbox(),
134            PySDM_products.AveragedTerminalVelocity(
135                name="rain averaged terminal velocity",
136                radius_range=settings.rain_water_radius_range,
137            ),
138            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
139            PySDM_products.AmbientPressure(name="p"),
140            PySDM_products.AmbientTemperature(name="T"),
141            PySDM_products.AmbientWaterVapourMixingRatio(
142                name="water_vapour_mixing_ratio"
143            ),
144        ]
145        if settings.enable_condensation:
146            self.products.extend(
147                [
148                    PySDM_products.RipeningRate(name="ripening"),
149                    PySDM_products.ActivatingRate(name="activating"),
150                    PySDM_products.DeactivatingRate(name="deactivating"),
151                    PySDM_products.PeakSupersaturation(unit="%"),
152                    PySDM_products.ParticleSizeSpectrumPerVolume(
153                        name="dry spectrum",
154                        radius_bins_edges=settings.r_bins_edges_dry,
155                        dry=True,
156                    ),
157                ]
158            )
159        if settings.precip:
160            self.products.extend(
161                [
162                    PySDM_products.CollisionRatePerGridbox(
163                        name="collision_rate",
164                    ),
165                    PySDM_products.CollisionRateDeficitPerGridbox(
166                        name="collision_deficit",
167                    ),
168                    PySDM_products.CoalescenceRatePerGridbox(
169                        name="coalescence_rate",
170                    ),
171                    PySDM_products.SurfacePrecipitation(),
172                ]
173            )
174        self.particulator = self.builder.build(
175            attributes=self.attributes, products=tuple(self.products)
176        )
177
178        self.output_attributes = {
179            "cell origin": [],
180            "position in cell": [],
181            "radius": [],
182            "multiplicity": [],
183        }
184        self.output_products = {}
185        for k, v in self.particulator.products.items():
186            if len(v.shape) == 0:
187                self.output_products[k] = np.zeros(self.nt + 1)
188            elif len(v.shape) == 1:
189                self.output_products[k] = np.zeros((self.mesh.grid[-1], self.nt + 1))
190            elif len(v.shape) == 2:
191                number_of_time_sections = len(self.save_spec_and_attr_times)
192                self.output_products[k] = np.zeros(
193                    (self.mesh.grid[-1], self.number_of_bins, number_of_time_sections)
194                )
195
196    @staticmethod
197    def add_collision_dynamic(builder, settings, _):
198        builder.add_dynamic(
199            Coalescence(
200                collision_kernel=settings.collision_kernel,
201                adaptive=settings.coalescence_adaptive,
202            )
203        )
204
205    def save_scalar(self, step):
206        for k, v in self.particulator.products.items():
207            if len(v.shape) > 1:
208                continue
209            if len(v.shape) == 1:
210                self.output_products[k][:, step] = v.get()
211            else:
212                self.output_products[k][step] = v.get()
213
214    def save_spectrum(self, index):
215        for k, v in self.particulator.products.items():
216            if len(v.shape) == 2:
217                self.output_products[k][:, :, index] = v.get()
218
219    def save_attributes(self):
220        for k, v in self.output_attributes.items():
221            v.append(self.particulator.attributes[k].to_ndarray())
222
223    def save(self, step):
224        self.save_scalar(step)
225        time = step * self.particulator.dt
226        if len(self.save_spec_and_attr_times) > 0 and (
227            np.min(
228                np.abs(
229                    np.ones_like(self.save_spec_and_attr_times) * time
230                    - np.array(self.save_spec_and_attr_times)
231                )
232            )
233            < 0.1
234        ):
235            save_index = np.argmin(
236                np.abs(
237                    np.ones_like(self.save_spec_and_attr_times) * time
238                    - np.array(self.save_spec_and_attr_times)
239                )
240            )
241            self.save_spectrum(save_index)
242            self.save_attributes()
243
244    def run(self):
245        mesh = self.particulator.mesh
246
247        assert "t" not in self.output_products and "z" not in self.output_products
248        self.output_products["t"] = np.linspace(
249            0, self.nt * self.particulator.dt, self.nt + 1, endpoint=True
250        )
251        self.output_products["z"] = np.linspace(
252            self.z0 + mesh.dz / 2,
253            self.z0 + (mesh.grid[-1] - 1 / 2) * mesh.dz,
254            mesh.grid[-1],
255            endpoint=True,
256        )
257
258        self.save(0)
259        for step in range(self.nt):
260            mpdata = self.particulator.dynamics["EulerianAdvection"].solvers
261            mpdata.update_advector_field()
262            if "Displacement" in self.particulator.dynamics:
263                self.particulator.dynamics["Displacement"].upload_courant_field(
264                    (mpdata.advector / self.g_factor_vec,)
265                )
266            self.particulator.run(steps=1)
267            self.save(step + 1)
268
269        Outputs = namedtuple("Outputs", "products attributes")
270        output_results = Outputs(self.output_products, self.output_attributes)
271        return output_results

class Simulation: View Source

 22class Simulation:
 23    def __init__(self, settings, backend=CPU):
 24        self.nt = settings.nt
 25        self.z0 = -settings.particle_reservoir_depth
 26        self.save_spec_and_attr_times = settings.save_spec_and_attr_times
 27        self.number_of_bins = settings.number_of_bins
 28
 29        self.particulator = None
 30        self.output_attributes = None
 31        self.output_products = None
 32
 33        self.mesh = Mesh(
 34            grid=(settings.nz,),
 35            size=(settings.z_max + settings.particle_reservoir_depth,),
 36        )
 37
 38        env = Kinematic1D(
 39            dt=settings.dt,
 40            mesh=self.mesh,
 41            thd_of_z=settings.thd,
 42            rhod_of_z=settings.rhod,
 43            z0=-settings.particle_reservoir_depth,
 44        )
 45
 46        def zZ_to_z_above_reservoir(zZ):
 47            z_above_reservoir = zZ * (settings.nz * settings.dz) + self.z0
 48            return z_above_reservoir
 49
 50        mpdata = MPDATA_1D(
 51            nz=settings.nz,
 52            dt=settings.dt,
 53            mpdata_settings=settings.mpdata_settings,
 54            advector_of_t=lambda t: settings.rho_times_w(t) * settings.dt / settings.dz,
 55            advectee_of_zZ_at_t0=lambda zZ: settings.water_vapour_mixing_ratio(
 56                zZ_to_z_above_reservoir(zZ)
 57            ),
 58            g_factor_of_zZ=lambda zZ: settings.rhod(zZ_to_z_above_reservoir(zZ)),
 59        )
 60
 61        _extra_nz = settings.particle_reservoir_depth // settings.dz
 62        _z_vec = settings.dz * np.linspace(
 63            -_extra_nz, settings.nz - _extra_nz, settings.nz + 1
 64        )
 65        self.g_factor_vec = settings.rhod(_z_vec)
 66
 67        self.builder = Builder(
 68            n_sd=settings.n_sd,
 69            backend=backend(formulae=settings.formulae),
 70            environment=env,
 71        )
 72        self.builder.add_dynamic(AmbientThermodynamics())
 73
 74        if settings.enable_condensation:
 75            self.builder.add_dynamic(
 76                Condensation(
 77                    adaptive=settings.condensation_adaptive,
 78                    rtol_thd=settings.condensation_rtol_thd,
 79                    rtol_x=settings.condensation_rtol_x,
 80                    update_thd=settings.condensation_update_thd,
 81                )
 82            )
 83        self.builder.add_dynamic(EulerianAdvection(mpdata))
 84
 85        self.products = []
 86        if settings.precip:
 87            self.add_collision_dynamic(self.builder, settings, self.products)
 88
 89        displacement = Displacement(
 90            enable_sedimentation=settings.precip,
 91            precipitation_counting_level_index=int(
 92                settings.particle_reservoir_depth / settings.dz
 93            ),
 94        )
 95        self.builder.add_dynamic(displacement)
 96        self.attributes = self.builder.particulator.environment.init_attributes(
 97            spatial_discretisation=spatial_sampling.Pseudorandom(),
 98            spectral_discretisation=spectral_sampling.ConstantMultiplicity(
 99                spectrum=settings.wet_radius_spectrum_per_mass_of_dry_air
100            ),
101            kappa=settings.kappa,
102            collisions_only=not settings.enable_condensation,
103            z_part=settings.z_part,
104        )
105        self.products += [
106            PySDM_products.WaterMixingRatio(
107                name="cloud water mixing ratio",
108                unit="g/kg",
109                radius_range=settings.cloud_water_radius_range,
110            ),
111            PySDM_products.WaterMixingRatio(
112                name="rain water mixing ratio",
113                unit="g/kg",
114                radius_range=settings.rain_water_radius_range,
115            ),
116            PySDM_products.AmbientDryAirDensity(name="rhod"),
117            PySDM_products.AmbientDryAirPotentialTemperature(name="thd"),
118            PySDM_products.ParticleSizeSpectrumPerVolume(
119                name="wet spectrum", radius_bins_edges=settings.r_bins_edges
120            ),
121            PySDM_products.ParticleConcentration(
122                name="nc", radius_range=settings.cloud_water_radius_range
123            ),
124            PySDM_products.ParticleConcentration(
125                name="nr", radius_range=settings.rain_water_radius_range
126            ),
127            PySDM_products.ParticleConcentration(
128                name="na", radius_range=(0, settings.cloud_water_radius_range[0])
129            ),
130            PySDM_products.MeanRadius(),
131            PySDM_products.EffectiveRadius(
132                radius_range=settings.cloud_water_radius_range
133            ),
134            PySDM_products.SuperDropletCountPerGridbox(),
135            PySDM_products.AveragedTerminalVelocity(
136                name="rain averaged terminal velocity",
137                radius_range=settings.rain_water_radius_range,
138            ),
139            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
140            PySDM_products.AmbientPressure(name="p"),
141            PySDM_products.AmbientTemperature(name="T"),
142            PySDM_products.AmbientWaterVapourMixingRatio(
143                name="water_vapour_mixing_ratio"
144            ),
145        ]
146        if settings.enable_condensation:
147            self.products.extend(
148                [
149                    PySDM_products.RipeningRate(name="ripening"),
150                    PySDM_products.ActivatingRate(name="activating"),
151                    PySDM_products.DeactivatingRate(name="deactivating"),
152                    PySDM_products.PeakSupersaturation(unit="%"),
153                    PySDM_products.ParticleSizeSpectrumPerVolume(
154                        name="dry spectrum",
155                        radius_bins_edges=settings.r_bins_edges_dry,
156                        dry=True,
157                    ),
158                ]
159            )
160        if settings.precip:
161            self.products.extend(
162                [
163                    PySDM_products.CollisionRatePerGridbox(
164                        name="collision_rate",
165                    ),
166                    PySDM_products.CollisionRateDeficitPerGridbox(
167                        name="collision_deficit",
168                    ),
169                    PySDM_products.CoalescenceRatePerGridbox(
170                        name="coalescence_rate",
171                    ),
172                    PySDM_products.SurfacePrecipitation(),
173                ]
174            )
175        self.particulator = self.builder.build(
176            attributes=self.attributes, products=tuple(self.products)
177        )
178
179        self.output_attributes = {
180            "cell origin": [],
181            "position in cell": [],
182            "radius": [],
183            "multiplicity": [],
184        }
185        self.output_products = {}
186        for k, v in self.particulator.products.items():
187            if len(v.shape) == 0:
188                self.output_products[k] = np.zeros(self.nt + 1)
189            elif len(v.shape) == 1:
190                self.output_products[k] = np.zeros((self.mesh.grid[-1], self.nt + 1))
191            elif len(v.shape) == 2:
192                number_of_time_sections = len(self.save_spec_and_attr_times)
193                self.output_products[k] = np.zeros(
194                    (self.mesh.grid[-1], self.number_of_bins, number_of_time_sections)
195                )
196
197    @staticmethod
198    def add_collision_dynamic(builder, settings, _):
199        builder.add_dynamic(
200            Coalescence(
201                collision_kernel=settings.collision_kernel,
202                adaptive=settings.coalescence_adaptive,
203            )
204        )
205
206    def save_scalar(self, step):
207        for k, v in self.particulator.products.items():
208            if len(v.shape) > 1:
209                continue
210            if len(v.shape) == 1:
211                self.output_products[k][:, step] = v.get()
212            else:
213                self.output_products[k][step] = v.get()
214
215    def save_spectrum(self, index):
216        for k, v in self.particulator.products.items():
217            if len(v.shape) == 2:
218                self.output_products[k][:, :, index] = v.get()
219
220    def save_attributes(self):
221        for k, v in self.output_attributes.items():
222            v.append(self.particulator.attributes[k].to_ndarray())
223
224    def save(self, step):
225        self.save_scalar(step)
226        time = step * self.particulator.dt
227        if len(self.save_spec_and_attr_times) > 0 and (
228            np.min(
229                np.abs(
230                    np.ones_like(self.save_spec_and_attr_times) * time
231                    - np.array(self.save_spec_and_attr_times)
232                )
233            )
234            < 0.1
235        ):
236            save_index = np.argmin(
237                np.abs(
238                    np.ones_like(self.save_spec_and_attr_times) * time
239                    - np.array(self.save_spec_and_attr_times)
240                )
241            )
242            self.save_spectrum(save_index)
243            self.save_attributes()
244
245    def run(self):
246        mesh = self.particulator.mesh
247
248        assert "t" not in self.output_products and "z" not in self.output_products
249        self.output_products["t"] = np.linspace(
250            0, self.nt * self.particulator.dt, self.nt + 1, endpoint=True
251        )
252        self.output_products["z"] = np.linspace(
253            self.z0 + mesh.dz / 2,
254            self.z0 + (mesh.grid[-1] - 1 / 2) * mesh.dz,
255            mesh.grid[-1],
256            endpoint=True,
257        )
258
259        self.save(0)
260        for step in range(self.nt):
261            mpdata = self.particulator.dynamics["EulerianAdvection"].solvers
262            mpdata.update_advector_field()
263            if "Displacement" in self.particulator.dynamics:
264                self.particulator.dynamics["Displacement"].upload_courant_field(
265                    (mpdata.advector / self.g_factor_vec,)
266                )
267            self.particulator.run(steps=1)
268            self.save(step + 1)
269
270        Outputs = namedtuple("Outputs", "products attributes")
271        output_results = Outputs(self.output_products, self.output_attributes)
272        return output_results

Simulation(settings, backend=<class 'PySDM.backends.numba.Numba'>) View Source

 23    def __init__(self, settings, backend=CPU):
 24        self.nt = settings.nt
 25        self.z0 = -settings.particle_reservoir_depth
 26        self.save_spec_and_attr_times = settings.save_spec_and_attr_times
 27        self.number_of_bins = settings.number_of_bins
 28
 29        self.particulator = None
 30        self.output_attributes = None
 31        self.output_products = None
 32
 33        self.mesh = Mesh(
 34            grid=(settings.nz,),
 35            size=(settings.z_max + settings.particle_reservoir_depth,),
 36        )
 37
 38        env = Kinematic1D(
 39            dt=settings.dt,
 40            mesh=self.mesh,
 41            thd_of_z=settings.thd,
 42            rhod_of_z=settings.rhod,
 43            z0=-settings.particle_reservoir_depth,
 44        )
 45
 46        def zZ_to_z_above_reservoir(zZ):
 47            z_above_reservoir = zZ * (settings.nz * settings.dz) + self.z0
 48            return z_above_reservoir
 49
 50        mpdata = MPDATA_1D(
 51            nz=settings.nz,
 52            dt=settings.dt,
 53            mpdata_settings=settings.mpdata_settings,
 54            advector_of_t=lambda t: settings.rho_times_w(t) * settings.dt / settings.dz,
 55            advectee_of_zZ_at_t0=lambda zZ: settings.water_vapour_mixing_ratio(
 56                zZ_to_z_above_reservoir(zZ)
 57            ),
 58            g_factor_of_zZ=lambda zZ: settings.rhod(zZ_to_z_above_reservoir(zZ)),
 59        )
 60
 61        _extra_nz = settings.particle_reservoir_depth // settings.dz
 62        _z_vec = settings.dz * np.linspace(
 63            -_extra_nz, settings.nz - _extra_nz, settings.nz + 1
 64        )
 65        self.g_factor_vec = settings.rhod(_z_vec)
 66
 67        self.builder = Builder(
 68            n_sd=settings.n_sd,
 69            backend=backend(formulae=settings.formulae),
 70            environment=env,
 71        )
 72        self.builder.add_dynamic(AmbientThermodynamics())
 73
 74        if settings.enable_condensation:
 75            self.builder.add_dynamic(
 76                Condensation(
 77                    adaptive=settings.condensation_adaptive,
 78                    rtol_thd=settings.condensation_rtol_thd,
 79                    rtol_x=settings.condensation_rtol_x,
 80                    update_thd=settings.condensation_update_thd,
 81                )
 82            )
 83        self.builder.add_dynamic(EulerianAdvection(mpdata))
 84
 85        self.products = []
 86        if settings.precip:
 87            self.add_collision_dynamic(self.builder, settings, self.products)
 88
 89        displacement = Displacement(
 90            enable_sedimentation=settings.precip,
 91            precipitation_counting_level_index=int(
 92                settings.particle_reservoir_depth / settings.dz
 93            ),
 94        )
 95        self.builder.add_dynamic(displacement)
 96        self.attributes = self.builder.particulator.environment.init_attributes(
 97            spatial_discretisation=spatial_sampling.Pseudorandom(),
 98            spectral_discretisation=spectral_sampling.ConstantMultiplicity(
 99                spectrum=settings.wet_radius_spectrum_per_mass_of_dry_air
100            ),
101            kappa=settings.kappa,
102            collisions_only=not settings.enable_condensation,
103            z_part=settings.z_part,
104        )
105        self.products += [
106            PySDM_products.WaterMixingRatio(
107                name="cloud water mixing ratio",
108                unit="g/kg",
109                radius_range=settings.cloud_water_radius_range,
110            ),
111            PySDM_products.WaterMixingRatio(
112                name="rain water mixing ratio",
113                unit="g/kg",
114                radius_range=settings.rain_water_radius_range,
115            ),
116            PySDM_products.AmbientDryAirDensity(name="rhod"),
117            PySDM_products.AmbientDryAirPotentialTemperature(name="thd"),
118            PySDM_products.ParticleSizeSpectrumPerVolume(
119                name="wet spectrum", radius_bins_edges=settings.r_bins_edges
120            ),
121            PySDM_products.ParticleConcentration(
122                name="nc", radius_range=settings.cloud_water_radius_range
123            ),
124            PySDM_products.ParticleConcentration(
125                name="nr", radius_range=settings.rain_water_radius_range
126            ),
127            PySDM_products.ParticleConcentration(
128                name="na", radius_range=(0, settings.cloud_water_radius_range[0])
129            ),
130            PySDM_products.MeanRadius(),
131            PySDM_products.EffectiveRadius(
132                radius_range=settings.cloud_water_radius_range
133            ),
134            PySDM_products.SuperDropletCountPerGridbox(),
135            PySDM_products.AveragedTerminalVelocity(
136                name="rain averaged terminal velocity",
137                radius_range=settings.rain_water_radius_range,
138            ),
139            PySDM_products.AmbientRelativeHumidity(name="RH", unit="%"),
140            PySDM_products.AmbientPressure(name="p"),
141            PySDM_products.AmbientTemperature(name="T"),
142            PySDM_products.AmbientWaterVapourMixingRatio(
143                name="water_vapour_mixing_ratio"
144            ),
145        ]
146        if settings.enable_condensation:
147            self.products.extend(
148                [
149                    PySDM_products.RipeningRate(name="ripening"),
150                    PySDM_products.ActivatingRate(name="activating"),
151                    PySDM_products.DeactivatingRate(name="deactivating"),
152                    PySDM_products.PeakSupersaturation(unit="%"),
153                    PySDM_products.ParticleSizeSpectrumPerVolume(
154                        name="dry spectrum",
155                        radius_bins_edges=settings.r_bins_edges_dry,
156                        dry=True,
157                    ),
158                ]
159            )
160        if settings.precip:
161            self.products.extend(
162                [
163                    PySDM_products.CollisionRatePerGridbox(
164                        name="collision_rate",
165                    ),
166                    PySDM_products.CollisionRateDeficitPerGridbox(
167                        name="collision_deficit",
168                    ),
169                    PySDM_products.CoalescenceRatePerGridbox(
170                        name="coalescence_rate",
171                    ),
172                    PySDM_products.SurfacePrecipitation(),
173                ]
174            )
175        self.particulator = self.builder.build(
176            attributes=self.attributes, products=tuple(self.products)
177        )
178
179        self.output_attributes = {
180            "cell origin": [],
181            "position in cell": [],
182            "radius": [],
183            "multiplicity": [],
184        }
185        self.output_products = {}
186        for k, v in self.particulator.products.items():
187            if len(v.shape) == 0:
188                self.output_products[k] = np.zeros(self.nt + 1)
189            elif len(v.shape) == 1:
190                self.output_products[k] = np.zeros((self.mesh.grid[-1], self.nt + 1))
191            elif len(v.shape) == 2:
192                number_of_time_sections = len(self.save_spec_and_attr_times)
193                self.output_products[k] = np.zeros(
194                    (self.mesh.grid[-1], self.number_of_bins, number_of_time_sections)
195                )

save_spec_and_attr_times

number_of_bins

particulator

output_attributes

output_products

mesh

g_factor_vec

builder

products

attributes

@staticmethod

def add_collision_dynamic(builder, settings, _): View Source

197    @staticmethod
198    def add_collision_dynamic(builder, settings, _):
199        builder.add_dynamic(
200            Coalescence(
201                collision_kernel=settings.collision_kernel,
202                adaptive=settings.coalescence_adaptive,
203            )
204        )

def save_scalar(self, step): View Source

206    def save_scalar(self, step):
207        for k, v in self.particulator.products.items():
208            if len(v.shape) > 1:
209                continue
210            if len(v.shape) == 1:
211                self.output_products[k][:, step] = v.get()
212            else:
213                self.output_products[k][step] = v.get()

def save_spectrum(self, index): View Source

215    def save_spectrum(self, index):
216        for k, v in self.particulator.products.items():
217            if len(v.shape) == 2:
218                self.output_products[k][:, :, index] = v.get()

def save_attributes(self): View Source

220    def save_attributes(self):
221        for k, v in self.output_attributes.items():
222            v.append(self.particulator.attributes[k].to_ndarray())

def save(self, step): View Source

224    def save(self, step):
225        self.save_scalar(step)
226        time = step * self.particulator.dt
227        if len(self.save_spec_and_attr_times) > 0 and (
228            np.min(
229                np.abs(
230                    np.ones_like(self.save_spec_and_attr_times) * time
231                    - np.array(self.save_spec_and_attr_times)
232                )
233            )
234            < 0.1
235        ):
236            save_index = np.argmin(
237                np.abs(
238                    np.ones_like(self.save_spec_and_attr_times) * time
239                    - np.array(self.save_spec_and_attr_times)
240                )
241            )
242            self.save_spectrum(save_index)
243            self.save_attributes()

def run(self): View Source

245    def run(self):
246        mesh = self.particulator.mesh
247
248        assert "t" not in self.output_products and "z" not in self.output_products
249        self.output_products["t"] = np.linspace(
250            0, self.nt * self.particulator.dt, self.nt + 1, endpoint=True
251        )
252        self.output_products["z"] = np.linspace(
253            self.z0 + mesh.dz / 2,
254            self.z0 + (mesh.grid[-1] - 1 / 2) * mesh.dz,
255            mesh.grid[-1],
256            endpoint=True,
257        )
258
259        self.save(0)
260        for step in range(self.nt):
261            mpdata = self.particulator.dynamics["EulerianAdvection"].solvers
262            mpdata.update_advector_field()
263            if "Displacement" in self.particulator.dynamics:
264                self.particulator.dynamics["Displacement"].upload_courant_field(
265                    (mpdata.advector / self.g_factor_vec,)
266                )
267            self.particulator.run(steps=1)
268            self.save(step + 1)
269
270        Outputs = namedtuple("Outputs", "products attributes")
271        output_results = Outputs(self.output_products, self.output_attributes)
272        return output_results