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            attributes["signed water mass"] = attributes.pop("water mass")
158
159            if self.settings.freezing_inp_spec is None:
160                immersed_surface_area = formulae.trivia.sphere_surface(
161                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
162                )
163            else:
164                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
165                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
166                )
167
168            if self.settings.freezing_singular:
169                attributes["freezing temperature"] = (
170                    formulae.freezing_temperature_spectrum.invcdf(
171                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
172                        immersed_surface_area,
173                    )
174                )
175            else:
176                attributes["immersed surface area"] = immersed_surface_area
177
178            if self.settings.freezing_inp_frac != 1:
179                assert self.settings.n_sd % 2 == 0
180                assert 0 < self.settings.freezing_inp_frac < 1
181                freezing_attribute = {
182                    True: "freezing temperature",
183                    False: "immersed surface area",
184                }[self.settings.freezing_singular]
185                for name, array in attributes.items():
186                    if array.shape[-1] != self.settings.n_sd // 2:
187                        raise AssertionError(f"attribute >>{name}<< has wrong size")
188                    array = array.copy()
189                    full_shape = list(array.shape)
190                    orig = slice(None, full_shape[-1])
191                    copy = slice(orig.stop, None)
192                    full_shape[-1] *= 2
193                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
194                    if name == freezing_attribute:
195                        attributes[name][orig] = array
196                        attributes[name][copy] = 0
197                    elif name == "multiplicity":
198                        attributes[name][orig] = array * self.settings.freezing_inp_frac
199                        attributes[name][copy] = array * (
200                            1 - self.settings.freezing_inp_frac
201                        )
202                    elif len(array.shape) > 1:  # particle positions
203                        # TODO #599: seed
204                        for dim, _ in enumerate(array.shape):
205                            # only to make particles not shadow each other in visualisations
206                            attributes[name][dim, orig] = array[dim, :]
207                            attributes[name][dim, copy] = np.random.permutation(
208                                array[dim, :]
209                            )
210                    else:
211                        attributes[name][orig] = array
212                        attributes[name][copy] = array
213
214                non_zero_per_gridbox = np.count_nonzero(
215                    attributes[freezing_attribute]
216                ) / np.prod(self.settings.grid)
217                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
218
219        self.particulator = builder.build(attributes, tuple(products))
220
221        if self.SpinUp is not None:
222            self.SpinUp(self.particulator, self.settings.n_spin_up)
223        if self.storage is not None:
224            self.storage.init(self.settings)
225
226    def run(self, controller=DummyController(), vtk_exporter=None):
227        with controller:
228            for step in self.settings.output_steps:
229                if controller.panic:
230                    break
231
232                self.particulator.run(step - self.particulator.n_steps)
233
234                self.store(step)
235
236                if vtk_exporter is not None:
237                    vtk_exporter.export_attributes(self.particulator)
238                    vtk_exporter.export_products(self.particulator)
239
240                controller.set_percent(step / self.settings.output_steps[-1])
241
242    def store(self, step):
243        for name, product in self.particulator.products.items():
244            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            attributes["signed water mass"] = attributes.pop("water mass")
159
160            if self.settings.freezing_inp_spec is None:
161                immersed_surface_area = formulae.trivia.sphere_surface(
162                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
163                )
164            else:
165                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
166                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
167                )
168
169            if self.settings.freezing_singular:
170                attributes["freezing temperature"] = (
171                    formulae.freezing_temperature_spectrum.invcdf(
172                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
173                        immersed_surface_area,
174                    )
175                )
176            else:
177                attributes["immersed surface area"] = immersed_surface_area
178
179            if self.settings.freezing_inp_frac != 1:
180                assert self.settings.n_sd % 2 == 0
181                assert 0 < self.settings.freezing_inp_frac < 1
182                freezing_attribute = {
183                    True: "freezing temperature",
184                    False: "immersed surface area",
185                }[self.settings.freezing_singular]
186                for name, array in attributes.items():
187                    if array.shape[-1] != self.settings.n_sd // 2:
188                        raise AssertionError(f"attribute >>{name}<< has wrong size")
189                    array = array.copy()
190                    full_shape = list(array.shape)
191                    orig = slice(None, full_shape[-1])
192                    copy = slice(orig.stop, None)
193                    full_shape[-1] *= 2
194                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
195                    if name == freezing_attribute:
196                        attributes[name][orig] = array
197                        attributes[name][copy] = 0
198                    elif name == "multiplicity":
199                        attributes[name][orig] = array * self.settings.freezing_inp_frac
200                        attributes[name][copy] = array * (
201                            1 - self.settings.freezing_inp_frac
202                        )
203                    elif len(array.shape) > 1:  # particle positions
204                        # TODO #599: seed
205                        for dim, _ in enumerate(array.shape):
206                            # only to make particles not shadow each other in visualisations
207                            attributes[name][dim, orig] = array[dim, :]
208                            attributes[name][dim, copy] = np.random.permutation(
209                                array[dim, :]
210                            )
211                    else:
212                        attributes[name][orig] = array
213                        attributes[name][copy] = array
214
215                non_zero_per_gridbox = np.count_nonzero(
216                    attributes[freezing_attribute]
217                ) / np.prod(self.settings.grid)
218                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
219
220        self.particulator = builder.build(attributes, tuple(products))
221
222        if self.SpinUp is not None:
223            self.SpinUp(self.particulator, self.settings.n_spin_up)
224        if self.storage is not None:
225            self.storage.init(self.settings)
226
227    def run(self, controller=DummyController(), vtk_exporter=None):
228        with controller:
229            for step in self.settings.output_steps:
230                if controller.panic:
231                    break
232
233                self.particulator.run(step - self.particulator.n_steps)
234
235                self.store(step)
236
237                if vtk_exporter is not None:
238                    vtk_exporter.export_attributes(self.particulator)
239                    vtk_exporter.export_products(self.particulator)
240
241                controller.set_percent(step / self.settings.output_steps[-1])
242
243    def store(self, step):
244        for name, product in self.particulator.products.items():
245            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            attributes["signed water mass"] = attributes.pop("water mass")
159
160            if self.settings.freezing_inp_spec is None:
161                immersed_surface_area = formulae.trivia.sphere_surface(
162                    diameter=2 * formulae.trivia.radius(volume=attributes["dry volume"])
163                )
164            else:
165                immersed_surface_area = self.settings.freezing_inp_spec.percentiles(
166                    np.random.random(attributes["dry volume"].size),  # TODO #599: seed
167                )
168
169            if self.settings.freezing_singular:
170                attributes["freezing temperature"] = (
171                    formulae.freezing_temperature_spectrum.invcdf(
172                        np.random.random(immersed_surface_area.size),  # TODO #599: seed
173                        immersed_surface_area,
174                    )
175                )
176            else:
177                attributes["immersed surface area"] = immersed_surface_area
178
179            if self.settings.freezing_inp_frac != 1:
180                assert self.settings.n_sd % 2 == 0
181                assert 0 < self.settings.freezing_inp_frac < 1
182                freezing_attribute = {
183                    True: "freezing temperature",
184                    False: "immersed surface area",
185                }[self.settings.freezing_singular]
186                for name, array in attributes.items():
187                    if array.shape[-1] != self.settings.n_sd // 2:
188                        raise AssertionError(f"attribute >>{name}<< has wrong size")
189                    array = array.copy()
190                    full_shape = list(array.shape)
191                    orig = slice(None, full_shape[-1])
192                    copy = slice(orig.stop, None)
193                    full_shape[-1] *= 2
194                    attributes[name] = np.empty(full_shape, dtype=array.dtype)
195                    if name == freezing_attribute:
196                        attributes[name][orig] = array
197                        attributes[name][copy] = 0
198                    elif name == "multiplicity":
199                        attributes[name][orig] = array * self.settings.freezing_inp_frac
200                        attributes[name][copy] = array * (
201                            1 - self.settings.freezing_inp_frac
202                        )
203                    elif len(array.shape) > 1:  # particle positions
204                        # TODO #599: seed
205                        for dim, _ in enumerate(array.shape):
206                            # only to make particles not shadow each other in visualisations
207                            attributes[name][dim, orig] = array[dim, :]
208                            attributes[name][dim, copy] = np.random.permutation(
209                                array[dim, :]
210                            )
211                    else:
212                        attributes[name][orig] = array
213                        attributes[name][copy] = array
214
215                non_zero_per_gridbox = np.count_nonzero(
216                    attributes[freezing_attribute]
217                ) / np.prod(self.settings.grid)
218                assert non_zero_per_gridbox == self.settings.n_sd_per_gridbox / 2
219
220        self.particulator = builder.build(attributes, tuple(products))
221
222        if self.SpinUp is not None:
223            self.SpinUp(self.particulator, self.settings.n_spin_up)
224        if self.storage is not None:
225            self.storage.init(self.settings)
def run( self, controller=<PySDM_examples.utils.dummy_controller.DummyController object>, vtk_exporter=None):
227    def run(self, controller=DummyController(), vtk_exporter=None):
228        with controller:
229            for step in self.settings.output_steps:
230                if controller.panic:
231                    break
232
233                self.particulator.run(step - self.particulator.n_steps)
234
235                self.store(step)
236
237                if vtk_exporter is not None:
238                    vtk_exporter.export_attributes(self.particulator)
239                    vtk_exporter.export_products(self.particulator)
240
241                controller.set_percent(step / self.settings.output_steps[-1])
def store(self, step):
243    def store(self, step):
244        for name, product in self.particulator.products.items():
245            self.storage.save(product.get(), step, name)