PySDM_examples.Szumowski_et_al_1998.simulation

  1import numpy as np
  2from PySDM_examples.Szumowski_et_al_1998.make_default_product_collection import (
  3    make_default_product_collection,
  4)
  5from PySDM_examples.Szumowski_et_al_1998.mpdata_2d import MPDATA_2D
  6from PySDM_examples.utils import DummyController
  7
  8from PySDM.backends import CPU
  9from PySDM.builder import Builder
 10from PySDM.dynamics import (
 11    AmbientThermodynamics,
 12    Coalescence,
 13    Collision,
 14    Condensation,
 15    Displacement,
 16    EulerianAdvection,
 17    Freezing,
 18)
 19from PySDM.environments import Kinematic2D
 20from PySDM.initialisation.sampling import spatial_sampling
 21
 22
 23class Simulation:
 24    def __init__(self, settings, storage, SpinUp, backend_class=CPU):
 25        self.settings = settings
 26        self.storage = storage
 27        self.particulator = None
 28        self.backend_class = backend_class
 29        self.SpinUp = SpinUp
 30
 31    @property
 32    def products(self):
 33        return self.particulator.products
 34
 35    def reinit(self, products=None):
 36        formulae = self.settings.formulae
 37        backend = self.backend_class(formulae=formulae)
 38        environment = Kinematic2D(
 39            dt=self.settings.dt,
 40            grid=self.settings.grid,
 41            size=self.settings.size,
 42            rhod_of=self.settings.rhod_of_zZ,
 43            mixed_phase=self.settings.processes["freezing"],
 44        )
 45        builder = Builder(
 46            n_sd=self.settings.n_sd, backend=backend, environment=environment
 47        )
 48
 49        if products is not None:
 50            products = list(products)
 51        else:
 52            products = make_default_product_collection(self.settings)
 53
 54        if self.settings.processes["fluid advection"]:
 55            builder.add_dynamic(AmbientThermodynamics())
 56        if self.settings.processes["condensation"]:
 57            kwargs = {}
 58            if not self.settings.condensation_adaptive:
 59                kwargs["substeps"] = (self.settings.condensation_substeps,)
 60            condensation = Condensation(
 61                rtol_x=self.settings.condensation_rtol_x,
 62                rtol_thd=self.settings.condensation_rtol_thd,
 63                adaptive=self.settings.condensation_adaptive,
 64                dt_cond_range=self.settings.condensation_dt_cond_range,
 65                schedule=self.settings.condensation_schedule,
 66                **kwargs,
 67            )
 68            builder.add_dynamic(condensation)
 69        if self.settings.processes["fluid advection"]:
 70            initial_profiles = {
 71                "th": self.settings.initial_dry_potential_temperature_profile,
 72                "water_vapour_mixing_ratio": self.settings.initial_vapour_mixing_ratio_profile,
 73            }
 74            advectees = dict(
 75                (
 76                    key,
 77                    np.repeat(
 78                        profile.reshape(1, -1),
 79                        builder.particulator.environment.mesh.grid[0],
 80                        axis=0,
 81                    ),
 82                )
 83                for key, profile in initial_profiles.items()
 84            )
 85            solver = MPDATA_2D(
 86                advectees=advectees,
 87                stream_function=self.settings.stream_function,
 88                rhod_of_zZ=self.settings.rhod_of_zZ,
 89                dt=self.settings.dt,
 90                grid=self.settings.grid,
 91                size=self.settings.size,
 92                n_iters=self.settings.mpdata_iters,
 93                infinite_gauge=self.settings.mpdata_iga,
 94                nonoscillatory=self.settings.mpdata_fct,
 95                third_order_terms=self.settings.mpdata_tot,
 96            )
 97            builder.add_dynamic(EulerianAdvection(solver))
 98        if self.settings.processes["particle advection"]:
 99            builder.add_dynamic(
100                Displacement(
101                    enable_sedimentation=self.settings.processes["sedimentation"],
102                    adaptive=self.settings.displacement_adaptive,
103                    rtol=self.settings.displacement_rtol,
104                )
105            )
106        if (
107            self.settings.processes["coalescence"]
108            and self.settings.processes["breakup"]
109        ):
110            builder.add_dynamic(
111                Collision(
112                    collision_kernel=self.settings.kernel,
113                    enable_breakup=self.settings.processes["breakup"],
114                    coalescence_efficiency=self.settings.coalescence_efficiency,
115                    breakup_efficiency=self.settings.breakup_efficiency,
116                    fragmentation_function=self.settings.breakup_fragmentation,
117                    adaptive=self.settings.coalescence_adaptive,
118                    dt_coal_range=self.settings.coalescence_dt_coal_range,
119                    substeps=self.settings.coalescence_substeps,
120                    optimized_random=self.settings.coalescence_optimized_random,
121                )
122            )
123        elif (
124            self.settings.processes["coalescence"]
125            and not self.settings.processes["breakup"]
126        ):
127            builder.add_dynamic(
128                Coalescence(
129                    collision_kernel=self.settings.kernel,
130                    adaptive=self.settings.coalescence_adaptive,
131                    dt_coal_range=self.settings.coalescence_dt_coal_range,
132                    substeps=self.settings.coalescence_substeps,
133                    optimized_random=self.settings.coalescence_optimized_random,
134                )
135            )
136        assert not (
137            self.settings.processes["breakup"]
138            and not self.settings.processes["coalescence"]
139        )
140        if self.settings.processes["freezing"]:
141            builder.add_dynamic(
142                Freezing(
143                    singular=self.settings.freezing_singular,
144                    thaw=self.settings.freezing_thaw,
145                )
146            )
147
148        attributes = builder.particulator.environment.init_attributes(
149            spatial_discretisation=spatial_sampling.Pseudorandom(),
150            dry_radius_spectrum=self.settings.spectrum_per_mass_of_dry_air,
151            kappa=self.settings.kappa,
152            n_sd=self.settings.n_sd
153            // (2 if self.settings.freezing_inp_frac != 1 else 1),
154        )
155
156        if self.settings.processes["freezing"]:
157            if self.settings.freezing_inp_spec is None:
158                immersed_surface_area = formulae.trivia.sphere_surface(
159                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
160                )
161            else:
162                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
163                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
164                )
165
166            if self.settings.freezing_singular:
167                attributes["freezing temperature"] = (
168                    formulae.freezing_temperature_spectrum.invcdf(
169                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
170                        immersed_surface_area,
171                    )
172                )
173            else:
174                attributes["immersed surface area"] = immersed_surface_area
175
176            if self.settings.freezing_inp_frac != 1:
177                assert self.settings.n_sd % 2 == 0
178                assert 0 < self.settings.freezing_inp_frac < 1
179                freezing_attribute = {
180                    True: "freezing temperature",
181                    False: "immersed surface area",
182                }[self.settings.freezing_singular]
183                for name, array in attributes.items():
184                    if array.shape[-1] != self.settings.n_sd // 2:
185                        raise AssertionError(f"attribute >>{name}<< has wrong size")
186                    array = array.copy()
187                    full_shape = list(array.shape)
188                    orig = slice(None, full_shape[-1])
189                    copy = slice(orig.stop, None)
190                    full_shape[-1] *= 2
191                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
192                    if name == freezing_attribute:
193                        attributes[name][orig] = array
194                        attributes[name][copy] = 0
195                    elif name == "multiplicity":
196                        attributes[name][orig] = array * self.settings.freezing_inp_frac
197                        attributes[name][copy] = array * (
198                            1 - self.settings.freezing_inp_frac
199                        )
200                    elif len(array.shape) > 1:  # particle positions
201                        # TODO #599: seed
202                        for dim, _ in enumerate(array.shape):
203                            # only to make particles not shadow each other in visualisations
204                            attributes[name][dim, orig] = array[dim, :]
205                            attributes[name][dim, copy] = np.random.permutation(
206                                array[dim, :]
207                            )
208                    else:
209                        attributes[name][orig] = array
210                        attributes[name][copy] = array
211
212                non_zero_per_gridbox = np.count_nonzero(
213                    attributes[freezing_attribute]
214                ) / np.prod(self.settings.grid)
215                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
216
217        self.particulator = builder.build(attributes, tuple(products))
218
219        if self.SpinUp is not None:
220            self.SpinUp(self.particulator, self.settings.n_spin_up)
221        if self.storage is not None:
222            self.storage.init(self.settings)
223
224    def run(self, controller=DummyController(), vtk_exporter=None):
225        with controller:
226            for step in self.settings.output_steps:
227                if controller.panic:
228                    break
229
230                self.particulator.run(step - self.particulator.n_steps)
231
232                self.store(step)
233
234                if vtk_exporter is not None:
235                    vtk_exporter.export_attributes(self.particulator)
236                    vtk_exporter.export_products(self.particulator)
237
238                controller.set_percent(step / self.settings.output_steps[-1])
239
240    def store(self, step):
241        for name, product in self.particulator.products.items():
242            self.storage.save(product.get(), step, name)
class Simulation:
 24class Simulation:
 25    def __init__(self, settings, storage, SpinUp, backend_class=CPU):
 26        self.settings = settings
 27        self.storage = storage
 28        self.particulator = None
 29        self.backend_class = backend_class
 30        self.SpinUp = SpinUp
 31
 32    @property
 33    def products(self):
 34        return self.particulator.products
 35
 36    def reinit(self, products=None):
 37        formulae = self.settings.formulae
 38        backend = self.backend_class(formulae=formulae)
 39        environment = Kinematic2D(
 40            dt=self.settings.dt,
 41            grid=self.settings.grid,
 42            size=self.settings.size,
 43            rhod_of=self.settings.rhod_of_zZ,
 44            mixed_phase=self.settings.processes["freezing"],
 45        )
 46        builder = Builder(
 47            n_sd=self.settings.n_sd, backend=backend, environment=environment
 48        )
 49
 50        if products is not None:
 51            products = list(products)
 52        else:
 53            products = make_default_product_collection(self.settings)
 54
 55        if self.settings.processes["fluid advection"]:
 56            builder.add_dynamic(AmbientThermodynamics())
 57        if self.settings.processes["condensation"]:
 58            kwargs = {}
 59            if not self.settings.condensation_adaptive:
 60                kwargs["substeps"] = (self.settings.condensation_substeps,)
 61            condensation = Condensation(
 62                rtol_x=self.settings.condensation_rtol_x,
 63                rtol_thd=self.settings.condensation_rtol_thd,
 64                adaptive=self.settings.condensation_adaptive,
 65                dt_cond_range=self.settings.condensation_dt_cond_range,
 66                schedule=self.settings.condensation_schedule,
 67                **kwargs,
 68            )
 69            builder.add_dynamic(condensation)
 70        if self.settings.processes["fluid advection"]:
 71            initial_profiles = {
 72                "th": self.settings.initial_dry_potential_temperature_profile,
 73                "water_vapour_mixing_ratio": self.settings.initial_vapour_mixing_ratio_profile,
 74            }
 75            advectees = dict(
 76                (
 77                    key,
 78                    np.repeat(
 79                        profile.reshape(1, -1),
 80                        builder.particulator.environment.mesh.grid[0],
 81                        axis=0,
 82                    ),
 83                )
 84                for key, profile in initial_profiles.items()
 85            )
 86            solver = MPDATA_2D(
 87                advectees=advectees,
 88                stream_function=self.settings.stream_function,
 89                rhod_of_zZ=self.settings.rhod_of_zZ,
 90                dt=self.settings.dt,
 91                grid=self.settings.grid,
 92                size=self.settings.size,
 93                n_iters=self.settings.mpdata_iters,
 94                infinite_gauge=self.settings.mpdata_iga,
 95                nonoscillatory=self.settings.mpdata_fct,
 96                third_order_terms=self.settings.mpdata_tot,
 97            )
 98            builder.add_dynamic(EulerianAdvection(solver))
 99        if self.settings.processes["particle advection"]:
100            builder.add_dynamic(
101                Displacement(
102                    enable_sedimentation=self.settings.processes["sedimentation"],
103                    adaptive=self.settings.displacement_adaptive,
104                    rtol=self.settings.displacement_rtol,
105                )
106            )
107        if (
108            self.settings.processes["coalescence"]
109            and self.settings.processes["breakup"]
110        ):
111            builder.add_dynamic(
112                Collision(
113                    collision_kernel=self.settings.kernel,
114                    enable_breakup=self.settings.processes["breakup"],
115                    coalescence_efficiency=self.settings.coalescence_efficiency,
116                    breakup_efficiency=self.settings.breakup_efficiency,
117                    fragmentation_function=self.settings.breakup_fragmentation,
118                    adaptive=self.settings.coalescence_adaptive,
119                    dt_coal_range=self.settings.coalescence_dt_coal_range,
120                    substeps=self.settings.coalescence_substeps,
121                    optimized_random=self.settings.coalescence_optimized_random,
122                )
123            )
124        elif (
125            self.settings.processes["coalescence"]
126            and not self.settings.processes["breakup"]
127        ):
128            builder.add_dynamic(
129                Coalescence(
130                    collision_kernel=self.settings.kernel,
131                    adaptive=self.settings.coalescence_adaptive,
132                    dt_coal_range=self.settings.coalescence_dt_coal_range,
133                    substeps=self.settings.coalescence_substeps,
134                    optimized_random=self.settings.coalescence_optimized_random,
135                )
136            )
137        assert not (
138            self.settings.processes["breakup"]
139            and not self.settings.processes["coalescence"]
140        )
141        if self.settings.processes["freezing"]:
142            builder.add_dynamic(
143                Freezing(
144                    singular=self.settings.freezing_singular,
145                    thaw=self.settings.freezing_thaw,
146                )
147            )
148
149        attributes = builder.particulator.environment.init_attributes(
150            spatial_discretisation=spatial_sampling.Pseudorandom(),
151            dry_radius_spectrum=self.settings.spectrum_per_mass_of_dry_air,
152            kappa=self.settings.kappa,
153            n_sd=self.settings.n_sd
154            // (2 if self.settings.freezing_inp_frac != 1 else 1),
155        )
156
157        if self.settings.processes["freezing"]:
158            if self.settings.freezing_inp_spec is None:
159                immersed_surface_area = formulae.trivia.sphere_surface(
160                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
161                )
162            else:
163                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
164                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
165                )
166
167            if self.settings.freezing_singular:
168                attributes["freezing temperature"] = (
169                    formulae.freezing_temperature_spectrum.invcdf(
170                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
171                        immersed_surface_area,
172                    )
173                )
174            else:
175                attributes["immersed surface area"] = immersed_surface_area
176
177            if self.settings.freezing_inp_frac != 1:
178                assert self.settings.n_sd % 2 == 0
179                assert 0 < self.settings.freezing_inp_frac < 1
180                freezing_attribute = {
181                    True: "freezing temperature",
182                    False: "immersed surface area",
183                }[self.settings.freezing_singular]
184                for name, array in attributes.items():
185                    if array.shape[-1] != self.settings.n_sd // 2:
186                        raise AssertionError(f"attribute >>{name}<< has wrong size")
187                    array = array.copy()
188                    full_shape = list(array.shape)
189                    orig = slice(None, full_shape[-1])
190                    copy = slice(orig.stop, None)
191                    full_shape[-1] *= 2
192                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
193                    if name == freezing_attribute:
194                        attributes[name][orig] = array
195                        attributes[name][copy] = 0
196                    elif name == "multiplicity":
197                        attributes[name][orig] = array * self.settings.freezing_inp_frac
198                        attributes[name][copy] = array * (
199                            1 - self.settings.freezing_inp_frac
200                        )
201                    elif len(array.shape) > 1:  # particle positions
202                        # TODO #599: seed
203                        for dim, _ in enumerate(array.shape):
204                            # only to make particles not shadow each other in visualisations
205                            attributes[name][dim, orig] = array[dim, :]
206                            attributes[name][dim, copy] = np.random.permutation(
207                                array[dim, :]
208                            )
209                    else:
210                        attributes[name][orig] = array
211                        attributes[name][copy] = array
212
213                non_zero_per_gridbox = np.count_nonzero(
214                    attributes[freezing_attribute]
215                ) / np.prod(self.settings.grid)
216                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
217
218        self.particulator = builder.build(attributes, tuple(products))
219
220        if self.SpinUp is not None:
221            self.SpinUp(self.particulator, self.settings.n_spin_up)
222        if self.storage is not None:
223            self.storage.init(self.settings)
224
225    def run(self, controller=DummyController(), vtk_exporter=None):
226        with controller:
227            for step in self.settings.output_steps:
228                if controller.panic:
229                    break
230
231                self.particulator.run(step - self.particulator.n_steps)
232
233                self.store(step)
234
235                if vtk_exporter is not None:
236                    vtk_exporter.export_attributes(self.particulator)
237                    vtk_exporter.export_products(self.particulator)
238
239                controller.set_percent(step / self.settings.output_steps[-1])
240
241    def store(self, step):
242        for name, product in self.particulator.products.items():
243            self.storage.save(product.get(), step, name)
Simulation( settings, storage, SpinUp, backend_class=<class 'PySDM.backends.numba.Numba'>)
25    def __init__(self, settings, storage, SpinUp, backend_class=CPU):
26        self.settings = settings
27        self.storage = storage
28        self.particulator = None
29        self.backend_class = backend_class
30        self.SpinUp = SpinUp
settings
storage
particulator
backend_class
SpinUp
products
32    @property
33    def products(self):
34        return self.particulator.products
def reinit(self, products=None):
 36    def reinit(self, products=None):
 37        formulae = self.settings.formulae
 38        backend = self.backend_class(formulae=formulae)
 39        environment = Kinematic2D(
 40            dt=self.settings.dt,
 41            grid=self.settings.grid,
 42            size=self.settings.size,
 43            rhod_of=self.settings.rhod_of_zZ,
 44            mixed_phase=self.settings.processes["freezing"],
 45        )
 46        builder = Builder(
 47            n_sd=self.settings.n_sd, backend=backend, environment=environment
 48        )
 49
 50        if products is not None:
 51            products = list(products)
 52        else:
 53            products = make_default_product_collection(self.settings)
 54
 55        if self.settings.processes["fluid advection"]:
 56            builder.add_dynamic(AmbientThermodynamics())
 57        if self.settings.processes["condensation"]:
 58            kwargs = {}
 59            if not self.settings.condensation_adaptive:
 60                kwargs["substeps"] = (self.settings.condensation_substeps,)
 61            condensation = Condensation(
 62                rtol_x=self.settings.condensation_rtol_x,
 63                rtol_thd=self.settings.condensation_rtol_thd,
 64                adaptive=self.settings.condensation_adaptive,
 65                dt_cond_range=self.settings.condensation_dt_cond_range,
 66                schedule=self.settings.condensation_schedule,
 67                **kwargs,
 68            )
 69            builder.add_dynamic(condensation)
 70        if self.settings.processes["fluid advection"]:
 71            initial_profiles = {
 72                "th": self.settings.initial_dry_potential_temperature_profile,
 73                "water_vapour_mixing_ratio": self.settings.initial_vapour_mixing_ratio_profile,
 74            }
 75            advectees = dict(
 76                (
 77                    key,
 78                    np.repeat(
 79                        profile.reshape(1, -1),
 80                        builder.particulator.environment.mesh.grid[0],
 81                        axis=0,
 82                    ),
 83                )
 84                for key, profile in initial_profiles.items()
 85            )
 86            solver = MPDATA_2D(
 87                advectees=advectees,
 88                stream_function=self.settings.stream_function,
 89                rhod_of_zZ=self.settings.rhod_of_zZ,
 90                dt=self.settings.dt,
 91                grid=self.settings.grid,
 92                size=self.settings.size,
 93                n_iters=self.settings.mpdata_iters,
 94                infinite_gauge=self.settings.mpdata_iga,
 95                nonoscillatory=self.settings.mpdata_fct,
 96                third_order_terms=self.settings.mpdata_tot,
 97            )
 98            builder.add_dynamic(EulerianAdvection(solver))
 99        if self.settings.processes["particle advection"]:
100            builder.add_dynamic(
101                Displacement(
102                    enable_sedimentation=self.settings.processes["sedimentation"],
103                    adaptive=self.settings.displacement_adaptive,
104                    rtol=self.settings.displacement_rtol,
105                )
106            )
107        if (
108            self.settings.processes["coalescence"]
109            and self.settings.processes["breakup"]
110        ):
111            builder.add_dynamic(
112                Collision(
113                    collision_kernel=self.settings.kernel,
114                    enable_breakup=self.settings.processes["breakup"],
115                    coalescence_efficiency=self.settings.coalescence_efficiency,
116                    breakup_efficiency=self.settings.breakup_efficiency,
117                    fragmentation_function=self.settings.breakup_fragmentation,
118                    adaptive=self.settings.coalescence_adaptive,
119                    dt_coal_range=self.settings.coalescence_dt_coal_range,
120                    substeps=self.settings.coalescence_substeps,
121                    optimized_random=self.settings.coalescence_optimized_random,
122                )
123            )
124        elif (
125            self.settings.processes["coalescence"]
126            and not self.settings.processes["breakup"]
127        ):
128            builder.add_dynamic(
129                Coalescence(
130                    collision_kernel=self.settings.kernel,
131                    adaptive=self.settings.coalescence_adaptive,
132                    dt_coal_range=self.settings.coalescence_dt_coal_range,
133                    substeps=self.settings.coalescence_substeps,
134                    optimized_random=self.settings.coalescence_optimized_random,
135                )
136            )
137        assert not (
138            self.settings.processes["breakup"]
139            and not self.settings.processes["coalescence"]
140        )
141        if self.settings.processes["freezing"]:
142            builder.add_dynamic(
143                Freezing(
144                    singular=self.settings.freezing_singular,
145                    thaw=self.settings.freezing_thaw,
146                )
147            )
148
149        attributes = builder.particulator.environment.init_attributes(
150            spatial_discretisation=spatial_sampling.Pseudorandom(),
151            dry_radius_spectrum=self.settings.spectrum_per_mass_of_dry_air,
152            kappa=self.settings.kappa,
153            n_sd=self.settings.n_sd
154            // (2 if self.settings.freezing_inp_frac != 1 else 1),
155        )
156
157        if self.settings.processes["freezing"]:
158            if self.settings.freezing_inp_spec is None:
159                immersed_surface_area = formulae.trivia.sphere_surface(
160                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
161                )
162            else:
163                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
164                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
165                )
166
167            if self.settings.freezing_singular:
168                attributes["freezing temperature"] = (
169                    formulae.freezing_temperature_spectrum.invcdf(
170                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
171                        immersed_surface_area,
172                    )
173                )
174            else:
175                attributes["immersed surface area"] = immersed_surface_area
176
177            if self.settings.freezing_inp_frac != 1:
178                assert self.settings.n_sd % 2 == 0
179                assert 0 < self.settings.freezing_inp_frac < 1
180                freezing_attribute = {
181                    True: "freezing temperature",
182                    False: "immersed surface area",
183                }[self.settings.freezing_singular]
184                for name, array in attributes.items():
185                    if array.shape[-1] != self.settings.n_sd // 2:
186                        raise AssertionError(f"attribute >>{name}<< has wrong size")
187                    array = array.copy()
188                    full_shape = list(array.shape)
189                    orig = slice(None, full_shape[-1])
190                    copy = slice(orig.stop, None)
191                    full_shape[-1] *= 2
192                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
193                    if name == freezing_attribute:
194                        attributes[name][orig] = array
195                        attributes[name][copy] = 0
196                    elif name == "multiplicity":
197                        attributes[name][orig] = array * self.settings.freezing_inp_frac
198                        attributes[name][copy] = array * (
199                            1 - self.settings.freezing_inp_frac
200                        )
201                    elif len(array.shape) > 1:  # particle positions
202                        # TODO #599: seed
203                        for dim, _ in enumerate(array.shape):
204                            # only to make particles not shadow each other in visualisations
205                            attributes[name][dim, orig] = array[dim, :]
206                            attributes[name][dim, copy] = np.random.permutation(
207                                array[dim, :]
208                            )
209                    else:
210                        attributes[name][orig] = array
211                        attributes[name][copy] = array
212
213                non_zero_per_gridbox = np.count_nonzero(
214                    attributes[freezing_attribute]
215                ) / np.prod(self.settings.grid)
216                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
217
218        self.particulator = builder.build(attributes, tuple(products))
219
220        if self.SpinUp is not None:
221            self.SpinUp(self.particulator, self.settings.n_spin_up)
222        if self.storage is not None:
223            self.storage.init(self.settings)
def run( self, controller=<PySDM_examples.utils.dummy_controller.DummyController object>, vtk_exporter=None):
225    def run(self, controller=DummyController(), vtk_exporter=None):
226        with controller:
227            for step in self.settings.output_steps:
228                if controller.panic:
229                    break
230
231                self.particulator.run(step - self.particulator.n_steps)
232
233                self.store(step)
234
235                if vtk_exporter is not None:
236                    vtk_exporter.export_attributes(self.particulator)
237                    vtk_exporter.export_products(self.particulator)
238
239                controller.set_percent(step / self.settings.output_steps[-1])
def store(self, step):
241    def store(self, step):
242        for name, product in self.particulator.products.items():
243            self.storage.save(product.get(), step, name)