PySDM.particulator
The very class exposing PySDM.particulator.Particulator.run() method for launching simulations
1""" 2The very class exposing `PySDM.particulator.Particulator.run()` method for launching simulations 3""" 4 5import numpy as np 6 7from PySDM.backends.impl_common.backend_methods import BackendMethods 8from PySDM.backends.impl_common.freezing_attributes import ( 9 SingularAttributes, 10 TimeDependentAttributes, 11) 12from PySDM.backends.impl_common.index import make_Index 13from PySDM.backends.impl_common.indexed_storage import make_IndexedStorage 14from PySDM.backends.impl_common.pair_indicator import make_PairIndicator 15from PySDM.backends.impl_common.pairwise_storage import make_PairwiseStorage 16from PySDM.impl.particle_attributes import ParticleAttributes 17 18 19class Particulator: # pylint: disable=too-many-public-methods,too-many-instance-attributes 20 def __init__(self, n_sd, backend): 21 assert isinstance(backend, BackendMethods) 22 self.__n_sd = n_sd 23 24 self.backend = backend 25 self.formulae = backend.formulae 26 self.environment = None 27 self.attributes: (ParticleAttributes, None) = None 28 self.dynamics = {} 29 self.products = {} 30 self.observers = [] 31 32 self.n_steps = 0 33 34 self.sorting_scheme = "default" 35 self.condensation_solver = None 36 37 self.Index = make_Index(backend) # pylint: disable=invalid-name 38 self.PairIndicator = make_PairIndicator(backend) # pylint: disable=invalid-name 39 self.PairwiseStorage = make_PairwiseStorage( # pylint: disable=invalid-name 40 backend 41 ) 42 self.IndexedStorage = make_IndexedStorage( # pylint: disable=invalid-name 43 backend 44 ) 45 46 self.timers = {} 47 self.null = self.Storage.empty(0, dtype=float) 48 49 def run(self, steps): 50 for _ in range(steps): 51 for key, dynamic in self.dynamics.items(): 52 with self.timers[key]: 53 dynamic() 54 self.n_steps += 1 55 self._notify_observers() 56 57 def _notify_observers(self): 58 reversed_order_so_that_environment_is_last = reversed(self.observers) 59 for observer in reversed_order_so_that_environment_is_last: 60 observer.notify() 61 62 @property 63 def Storage(self): 64 return self.backend.Storage 65 66 @property 67 def Random(self): 68 return self.backend.Random 69 70 @property 71 def n_sd(self) -> int: 72 return self.__n_sd 73 74 @property 75 def dt(self) -> float: 76 if self.environment is not None: 77 return self.environment.dt 78 return None 79 80 @property 81 def mesh(self): 82 if self.environment is not None: 83 return self.environment.mesh 84 return None 85 86 def normalize(self, prob, norm_factor): 87 self.backend.normalize( 88 prob=prob, 89 cell_id=self.attributes["cell id"], 90 cell_idx=self.attributes.cell_idx, 91 cell_start=self.attributes.cell_start, 92 norm_factor=norm_factor, 93 timestep=self.dt, 94 dv=self.mesh.dv, 95 ) 96 97 def update_TpRH(self): 98 self.backend.temperature_pressure_rh( 99 # input 100 rhod=self.environment.get_predicted("rhod"), 101 thd=self.environment.get_predicted("thd"), 102 water_vapour_mixing_ratio=self.environment.get_predicted( 103 "water_vapour_mixing_ratio" 104 ), 105 # output 106 T=self.environment.get_predicted("T"), 107 p=self.environment.get_predicted("p"), 108 RH=self.environment.get_predicted("RH"), 109 ) 110 111 def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): 112 """Updates droplet volumes by simulating condensation driven by prior changes 113 in environment thermodynamic state, updates the environment state. 114 In the case of parcel environment, condensation is driven solely by changes in 115 the dry-air density (theta and water_vapour_mixing_ratio should not be changed 116 by other dynamics). 117 In the case of prescribed-flow/kinematic environments, the dry-air density is 118 constant in time throughout the simulation. 119 This function should only change environment's predicted `thd` and 120 `water_vapour_mixing_ratio` (and not `rhod`). 121 """ 122 self.backend.condensation( 123 solver=self.condensation_solver, 124 n_cell=self.mesh.n_cell, 125 cell_start_arg=self.attributes.cell_start, 126 water_mass=self.attributes["signed water mass"], 127 multiplicity=self.attributes["multiplicity"], 128 vdry=self.attributes["dry volume"], 129 idx=self.attributes._ParticleAttributes__idx, 130 rhod=self.environment["rhod"], 131 thd=self.environment["thd"], 132 water_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 133 dv=self.environment.dv, 134 prhod=self.environment.get_predicted("rhod"), 135 pthd=self.environment.get_predicted("thd"), 136 predicted_water_vapour_mixing_ratio=self.environment.get_predicted( 137 "water_vapour_mixing_ratio" 138 ), 139 kappa=self.attributes["kappa"], 140 f_org=self.attributes["dry volume organic fraction"], 141 rtol_x=rtol_x, 142 rtol_thd=rtol_thd, 143 v_cr=self.attributes["critical volume"], 144 timestep=self.dt, 145 counters=counters, 146 cell_order=cell_order, 147 RH_max=RH_max, 148 success=success, 149 cell_id=self.attributes["cell id"], 150 reynolds_number=self.attributes["Reynolds number"], 151 air_density=self.environment["air density"], 152 air_dynamic_viscosity=self.environment["air dynamic viscosity"], 153 ) 154 self.attributes.mark_updated("signed water mass") 155 156 def collision_coalescence_breakup( 157 self, 158 *, 159 enable_breakup, 160 gamma, 161 rand, 162 Ec, 163 Eb, 164 fragment_mass, 165 coalescence_rate, 166 breakup_rate, 167 breakup_rate_deficit, 168 is_first_in_pair, 169 warn_overflows, 170 max_multiplicity, 171 ): 172 # pylint: disable=too-many-locals 173 idx = self.attributes._ParticleAttributes__idx 174 healthy = self.attributes._ParticleAttributes__healthy_memory 175 cell_id = self.attributes["cell id"] 176 multiplicity = self.attributes["multiplicity"] 177 attributes = self.attributes.get_extensive_attribute_storage() 178 if enable_breakup: 179 self.backend.collision_coalescence_breakup( 180 multiplicity=multiplicity, 181 idx=idx, 182 attributes=attributes, 183 gamma=gamma, 184 rand=rand, 185 Ec=Ec, 186 Eb=Eb, 187 fragment_mass=fragment_mass, 188 healthy=healthy, 189 cell_id=cell_id, 190 coalescence_rate=coalescence_rate, 191 breakup_rate=breakup_rate, 192 breakup_rate_deficit=breakup_rate_deficit, 193 is_first_in_pair=is_first_in_pair, 194 warn_overflows=warn_overflows, 195 particle_mass=self.attributes["water mass"], 196 max_multiplicity=max_multiplicity, 197 ) 198 else: 199 self.backend.collision_coalescence( 200 multiplicity=multiplicity, 201 idx=idx, 202 attributes=attributes, 203 gamma=gamma, 204 healthy=healthy, 205 cell_id=cell_id, 206 coalescence_rate=coalescence_rate, 207 is_first_in_pair=is_first_in_pair, 208 ) 209 self.attributes.sanitize() 210 self.attributes.mark_updated("multiplicity") 211 for key in self.attributes.get_extensive_attribute_keys(): 212 self.attributes.mark_updated(key) 213 214 def oxidation( 215 self, 216 *, 217 kinetic_consts, 218 timestep, 219 equilibrium_consts, 220 dissociation_factors, 221 do_chemistry_flag, 222 ): 223 self.backend.oxidation( 224 n_sd=self.n_sd, 225 cell_ids=self.attributes["cell id"], 226 do_chemistry_flag=do_chemistry_flag, 227 k0=kinetic_consts["k0"], 228 k1=kinetic_consts["k1"], 229 k2=kinetic_consts["k2"], 230 k3=kinetic_consts["k3"], 231 K_SO2=equilibrium_consts["K_SO2"], 232 K_HSO3=equilibrium_consts["K_HSO3"], 233 dissociation_factor_SO2=dissociation_factors["SO2"], 234 timestep=timestep, 235 # input 236 droplet_volume=self.attributes["volume"], 237 pH=self.attributes["pH"], 238 # output 239 moles_O3=self.attributes["moles_O3"], 240 moles_H2O2=self.attributes["moles_H2O2"], 241 moles_S_IV=self.attributes["moles_S_IV"], 242 moles_S_VI=self.attributes["moles_S_VI"], 243 ) 244 for attr in ("moles_S_IV", "moles_S_VI", "moles_H2O2", "moles_O3"): 245 self.attributes.mark_updated(attr) 246 247 def dissolution( 248 self, 249 *, 250 gaseous_compounds, 251 system_type, 252 dissociation_factors, 253 timestep, 254 environment_mixing_ratios, 255 do_chemistry_flag, 256 ): 257 self.backend.dissolution( 258 n_cell=self.mesh.n_cell, 259 n_threads=1, 260 cell_order=np.arange(self.mesh.n_cell), 261 cell_start_arg=self.attributes.cell_start, 262 idx=self.attributes._ParticleAttributes__idx, 263 do_chemistry_flag=do_chemistry_flag, 264 mole_amounts={ 265 key: self.attributes["moles_" + key] for key in gaseous_compounds.keys() 266 }, 267 env_mixing_ratio=environment_mixing_ratios, 268 # note: assuming condensation was called 269 env_p=self.environment.get_predicted("p"), 270 env_T=self.environment.get_predicted("T"), 271 env_rho_d=self.environment.get_predicted("rhod"), 272 timestep=timestep, 273 dv=self.mesh.dv, 274 droplet_volume=self.attributes["volume"], 275 multiplicity=self.attributes["multiplicity"], 276 system_type=system_type, 277 dissociation_factors=dissociation_factors, 278 ) 279 for key in gaseous_compounds.keys(): 280 self.attributes.mark_updated(f"moles_{key}") 281 282 def chem_recalculate_cell_data(self, *, equilibrium_consts, kinetic_consts): 283 self.backend.chem_recalculate_cell_data( 284 equilibrium_consts=equilibrium_consts, 285 kinetic_consts=kinetic_consts, 286 temperature=self.environment.get_predicted("T"), 287 ) 288 289 def chem_recalculate_drop_data(self, *, dissociation_factors, equilibrium_consts): 290 self.backend.chem_recalculate_drop_data( 291 dissociation_factors=dissociation_factors, 292 equilibrium_consts=equilibrium_consts, 293 pH=self.attributes["pH"], 294 cell_id=self.attributes["cell id"], 295 ) 296 297 def recalculate_cell_id(self): 298 if not self.attributes.has_attribute("cell origin"): 299 return 300 self.backend.cell_id( 301 self.attributes["cell id"], 302 self.attributes["cell origin"], 303 self.backend.Storage.from_ndarray(self.environment.mesh.strides), 304 ) 305 self.attributes._ParticleAttributes__sorted = False 306 307 def sort_within_pair_by_attr(self, is_first_in_pair, attr_name): 308 self.backend.sort_within_pair_by_attr( 309 self.attributes._ParticleAttributes__idx, 310 is_first_in_pair, 311 self.attributes[attr_name], 312 ) 313 314 def moments( 315 self, 316 *, 317 moment_0, 318 moments, 319 specs: dict, 320 attr_name="signed water mass", 321 attr_range=(-np.inf, np.inf), 322 weighting_attribute="water mass", 323 weighting_rank=0, 324 skip_division_by_m0=False, 325 ): 326 """ 327 Writes to `moment_0` and `moment` the zero-th and the k-th statistical moments 328 of particle attributes computed filtering by value of the attribute `attr_name` 329 to fall within `attr_range`. The moment ranks are defined by `specs`. 330 331 Parameters: 332 specs: e.g., `specs={'volume': (1,2,3), 'kappa': (1)}` computes three moments 333 of volume and one moment of kappa 334 skip_division_by_m0: if set to `True`, the values written to `moments` are 335 multiplied by the 0-th moment (e.g., total volume instead of mean volume) 336 """ 337 if len(specs) == 0: 338 raise ValueError("empty specs passed") 339 attr_data, ranks = [], [] 340 for attr in specs: 341 for rank in specs[attr]: 342 attr_data.append(self.attributes[attr]) 343 ranks.append(rank) 344 assert len(set(attr_data)) <= 1 345 if len(attr_data) == 0: 346 attr_data = self.backend.Storage.empty((0,), dtype=float) 347 else: 348 attr_data = attr_data[0] 349 350 ranks = self.backend.Storage.from_ndarray(np.array(ranks, dtype=float)) 351 352 self.backend.moments( 353 moment_0=moment_0, 354 moments=moments, 355 multiplicity=self.attributes["multiplicity"], 356 attr_data=attr_data, 357 cell_id=self.attributes["cell id"], 358 idx=self.attributes._ParticleAttributes__idx, 359 length=self.attributes.super_droplet_count, 360 ranks=ranks, 361 min_x=attr_range[0], 362 max_x=attr_range[1], 363 x_attr=self.attributes[attr_name], 364 weighting_attribute=self.attributes[weighting_attribute], 365 weighting_rank=weighting_rank, 366 skip_division_by_m0=skip_division_by_m0, 367 ) 368 369 def spectrum_moments( 370 self, 371 *, 372 moment_0, 373 moments, 374 attr, 375 rank, 376 attr_bins, 377 attr_name="water mass", 378 weighting_attribute="water mass", 379 weighting_rank=0, 380 ): 381 attr_data = self.attributes[attr] 382 self.backend.spectrum_moments( 383 moment_0=moment_0, 384 moments=moments, 385 multiplicity=self.attributes["multiplicity"], 386 attr_data=attr_data, 387 cell_id=self.attributes["cell id"], 388 idx=self.attributes._ParticleAttributes__idx, 389 length=self.attributes.super_droplet_count, 390 rank=rank, 391 x_bins=attr_bins, 392 x_attr=self.attributes[attr_name], 393 weighting_attribute=self.attributes[weighting_attribute], 394 weighting_rank=weighting_rank, 395 ) 396 397 def adaptive_sdm_end(self, dt_left): 398 return self.backend.adaptive_sdm_end(dt_left, self.attributes.cell_start) 399 400 def remove_precipitated( 401 self, *, displacement, precipitation_counting_level_index 402 ) -> float: 403 rainfall_mass = self.backend.flag_precipitated( 404 cell_origin=self.attributes["cell origin"], 405 position_in_cell=self.attributes["position in cell"], 406 water_mass=self.attributes["water mass"], 407 multiplicity=self.attributes["multiplicity"], 408 idx=self.attributes._ParticleAttributes__idx, 409 length=self.attributes.super_droplet_count, 410 healthy=self.attributes._ParticleAttributes__healthy_memory, 411 precipitation_counting_level_index=precipitation_counting_level_index, 412 displacement=displacement, 413 ) 414 self.attributes.sanitize() 415 return rainfall_mass 416 417 def flag_out_of_column(self): 418 self.backend.flag_out_of_column( 419 cell_origin=self.attributes["cell origin"], 420 position_in_cell=self.attributes["position in cell"], 421 idx=self.attributes._ParticleAttributes__idx, 422 length=self.attributes.super_droplet_count, 423 healthy=self.attributes._ParticleAttributes__healthy_memory, 424 domain_top_level_index=self.mesh.grid[-1], 425 ) 426 self.attributes.sanitize() 427 428 def calculate_displacement( 429 self, *, displacement, courant, cell_origin, position_in_cell, n_substeps 430 ): 431 for dim in range(len(self.environment.mesh.grid)): 432 self.backend.calculate_displacement( 433 dim=dim, 434 displacement=displacement, 435 courant=courant[dim], 436 cell_origin=cell_origin, 437 position_in_cell=position_in_cell, 438 n_substeps=n_substeps, 439 ) 440 441 def isotopic_fractionation(self, heavy_isotopes: tuple): 442 self.backend.isotopic_fractionation() 443 for isotope in heavy_isotopes: 444 self.attributes.mark_updated(f"moles_{isotope}") 445 446 def seeding( 447 self, 448 *, 449 seeded_particle_index, 450 seeded_particle_multiplicity, 451 seeded_particle_extensive_attributes, 452 number_of_super_particles_to_inject, 453 ): 454 n_null = self.n_sd - self.attributes.super_droplet_count 455 if n_null == 0: 456 raise ValueError( 457 "No available seeds to inject. Please provide particles with nan filled attributes." 458 ) 459 460 if number_of_super_particles_to_inject > n_null: 461 raise ValueError( 462 "Trying to inject more super particles than space available." 463 ) 464 465 if number_of_super_particles_to_inject > len(seeded_particle_multiplicity): 466 raise ValueError( 467 "Trying to inject multiple super particles with the same attributes. \ 468 Instead increase multiplicity of injected particles." 469 ) 470 471 self.backend.seeding( 472 idx=self.attributes._ParticleAttributes__idx, 473 multiplicity=self.attributes["multiplicity"], 474 extensive_attributes=self.attributes.get_extensive_attribute_storage(), 475 seeded_particle_index=seeded_particle_index, 476 seeded_particle_multiplicity=seeded_particle_multiplicity, 477 seeded_particle_extensive_attributes=seeded_particle_extensive_attributes, 478 number_of_super_particles_to_inject=number_of_super_particles_to_inject, 479 ) 480 self.attributes.reset_idx() 481 self.attributes.sanitize() 482 483 self.attributes.mark_updated("multiplicity") 484 for key in self.attributes.get_extensive_attribute_keys(): 485 self.attributes.mark_updated(key) 486 487 def deposition(self): 488 self.backend.deposition( 489 multiplicity=self.attributes["multiplicity"], 490 signed_water_mass=self.attributes["signed water mass"], 491 current_temperature=self.environment["T"], 492 current_total_pressure=self.environment["p"], 493 current_relative_humidity=self.environment["RH"], 494 current_water_activity=self.environment["a_w_ice"], 495 current_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 496 current_dry_air_density=self.environment["rhod"], 497 current_dry_potential_temperature=self.environment["thd"], 498 cell_volume=self.environment.mesh.dv, 499 time_step=self.dt, 500 cell_id=self.attributes["cell id"], 501 reynolds_number=self.attributes["Reynolds number"], 502 schmidt_number=self.environment["Schmidt number"], 503 predicted_vapour_mixing_ratio=self.environment.get_predicted( 504 "water_vapour_mixing_ratio" 505 ), 506 predicted_dry_potential_temperature=self.environment.get_predicted("thd"), 507 ) 508 self.attributes.mark_updated("signed water mass") 509 # TODO #1524 - should we update here? 510 # self.update_TpRH(only_if_not_last='VapourDepositionOnIce') 511 512 def immersion_freezing_time_dependent(self, *, thaw: bool, rand: Storage): 513 self.backend.freeze_time_dependent( 514 rand=rand, 515 attributes=TimeDependentAttributes( 516 immersed_surface_area=self.attributes["immersed surface area"], 517 signed_water_mass=self.attributes["signed water mass"], 518 ), 519 timestep=self.dt, 520 cell=self.attributes["cell id"], 521 a_w_ice=self.environment["a_w_ice"], 522 temperature=self.environment["T"], 523 relative_humidity=self.environment["RH"], 524 thaw=thaw, 525 ) 526 self.attributes.mark_updated("signed water mass") 527 528 def immersion_freezing_singular(self, *, thaw: bool): 529 self.backend.freeze_singular( 530 attributes=SingularAttributes( 531 freezing_temperature=self.attributes["freezing temperature"], 532 signed_water_mass=self.attributes["signed water mass"], 533 ), 534 temperature=self.environment["T"], 535 relative_humidity=self.environment["RH"], 536 cell=self.attributes["cell id"], 537 thaw=thaw, 538 ) 539 self.attributes.mark_updated("signed water mass")
class
Particulator:
20class Particulator: # pylint: disable=too-many-public-methods,too-many-instance-attributes 21 def __init__(self, n_sd, backend): 22 assert isinstance(backend, BackendMethods) 23 self.__n_sd = n_sd 24 25 self.backend = backend 26 self.formulae = backend.formulae 27 self.environment = None 28 self.attributes: (ParticleAttributes, None) = None 29 self.dynamics = {} 30 self.products = {} 31 self.observers = [] 32 33 self.n_steps = 0 34 35 self.sorting_scheme = "default" 36 self.condensation_solver = None 37 38 self.Index = make_Index(backend) # pylint: disable=invalid-name 39 self.PairIndicator = make_PairIndicator(backend) # pylint: disable=invalid-name 40 self.PairwiseStorage = make_PairwiseStorage( # pylint: disable=invalid-name 41 backend 42 ) 43 self.IndexedStorage = make_IndexedStorage( # pylint: disable=invalid-name 44 backend 45 ) 46 47 self.timers = {} 48 self.null = self.Storage.empty(0, dtype=float) 49 50 def run(self, steps): 51 for _ in range(steps): 52 for key, dynamic in self.dynamics.items(): 53 with self.timers[key]: 54 dynamic() 55 self.n_steps += 1 56 self._notify_observers() 57 58 def _notify_observers(self): 59 reversed_order_so_that_environment_is_last = reversed(self.observers) 60 for observer in reversed_order_so_that_environment_is_last: 61 observer.notify() 62 63 @property 64 def Storage(self): 65 return self.backend.Storage 66 67 @property 68 def Random(self): 69 return self.backend.Random 70 71 @property 72 def n_sd(self) -> int: 73 return self.__n_sd 74 75 @property 76 def dt(self) -> float: 77 if self.environment is not None: 78 return self.environment.dt 79 return None 80 81 @property 82 def mesh(self): 83 if self.environment is not None: 84 return self.environment.mesh 85 return None 86 87 def normalize(self, prob, norm_factor): 88 self.backend.normalize( 89 prob=prob, 90 cell_id=self.attributes["cell id"], 91 cell_idx=self.attributes.cell_idx, 92 cell_start=self.attributes.cell_start, 93 norm_factor=norm_factor, 94 timestep=self.dt, 95 dv=self.mesh.dv, 96 ) 97 98 def update_TpRH(self): 99 self.backend.temperature_pressure_rh( 100 # input 101 rhod=self.environment.get_predicted("rhod"), 102 thd=self.environment.get_predicted("thd"), 103 water_vapour_mixing_ratio=self.environment.get_predicted( 104 "water_vapour_mixing_ratio" 105 ), 106 # output 107 T=self.environment.get_predicted("T"), 108 p=self.environment.get_predicted("p"), 109 RH=self.environment.get_predicted("RH"), 110 ) 111 112 def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): 113 """Updates droplet volumes by simulating condensation driven by prior changes 114 in environment thermodynamic state, updates the environment state. 115 In the case of parcel environment, condensation is driven solely by changes in 116 the dry-air density (theta and water_vapour_mixing_ratio should not be changed 117 by other dynamics). 118 In the case of prescribed-flow/kinematic environments, the dry-air density is 119 constant in time throughout the simulation. 120 This function should only change environment's predicted `thd` and 121 `water_vapour_mixing_ratio` (and not `rhod`). 122 """ 123 self.backend.condensation( 124 solver=self.condensation_solver, 125 n_cell=self.mesh.n_cell, 126 cell_start_arg=self.attributes.cell_start, 127 water_mass=self.attributes["signed water mass"], 128 multiplicity=self.attributes["multiplicity"], 129 vdry=self.attributes["dry volume"], 130 idx=self.attributes._ParticleAttributes__idx, 131 rhod=self.environment["rhod"], 132 thd=self.environment["thd"], 133 water_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 134 dv=self.environment.dv, 135 prhod=self.environment.get_predicted("rhod"), 136 pthd=self.environment.get_predicted("thd"), 137 predicted_water_vapour_mixing_ratio=self.environment.get_predicted( 138 "water_vapour_mixing_ratio" 139 ), 140 kappa=self.attributes["kappa"], 141 f_org=self.attributes["dry volume organic fraction"], 142 rtol_x=rtol_x, 143 rtol_thd=rtol_thd, 144 v_cr=self.attributes["critical volume"], 145 timestep=self.dt, 146 counters=counters, 147 cell_order=cell_order, 148 RH_max=RH_max, 149 success=success, 150 cell_id=self.attributes["cell id"], 151 reynolds_number=self.attributes["Reynolds number"], 152 air_density=self.environment["air density"], 153 air_dynamic_viscosity=self.environment["air dynamic viscosity"], 154 ) 155 self.attributes.mark_updated("signed water mass") 156 157 def collision_coalescence_breakup( 158 self, 159 *, 160 enable_breakup, 161 gamma, 162 rand, 163 Ec, 164 Eb, 165 fragment_mass, 166 coalescence_rate, 167 breakup_rate, 168 breakup_rate_deficit, 169 is_first_in_pair, 170 warn_overflows, 171 max_multiplicity, 172 ): 173 # pylint: disable=too-many-locals 174 idx = self.attributes._ParticleAttributes__idx 175 healthy = self.attributes._ParticleAttributes__healthy_memory 176 cell_id = self.attributes["cell id"] 177 multiplicity = self.attributes["multiplicity"] 178 attributes = self.attributes.get_extensive_attribute_storage() 179 if enable_breakup: 180 self.backend.collision_coalescence_breakup( 181 multiplicity=multiplicity, 182 idx=idx, 183 attributes=attributes, 184 gamma=gamma, 185 rand=rand, 186 Ec=Ec, 187 Eb=Eb, 188 fragment_mass=fragment_mass, 189 healthy=healthy, 190 cell_id=cell_id, 191 coalescence_rate=coalescence_rate, 192 breakup_rate=breakup_rate, 193 breakup_rate_deficit=breakup_rate_deficit, 194 is_first_in_pair=is_first_in_pair, 195 warn_overflows=warn_overflows, 196 particle_mass=self.attributes["water mass"], 197 max_multiplicity=max_multiplicity, 198 ) 199 else: 200 self.backend.collision_coalescence( 201 multiplicity=multiplicity, 202 idx=idx, 203 attributes=attributes, 204 gamma=gamma, 205 healthy=healthy, 206 cell_id=cell_id, 207 coalescence_rate=coalescence_rate, 208 is_first_in_pair=is_first_in_pair, 209 ) 210 self.attributes.sanitize() 211 self.attributes.mark_updated("multiplicity") 212 for key in self.attributes.get_extensive_attribute_keys(): 213 self.attributes.mark_updated(key) 214 215 def oxidation( 216 self, 217 *, 218 kinetic_consts, 219 timestep, 220 equilibrium_consts, 221 dissociation_factors, 222 do_chemistry_flag, 223 ): 224 self.backend.oxidation( 225 n_sd=self.n_sd, 226 cell_ids=self.attributes["cell id"], 227 do_chemistry_flag=do_chemistry_flag, 228 k0=kinetic_consts["k0"], 229 k1=kinetic_consts["k1"], 230 k2=kinetic_consts["k2"], 231 k3=kinetic_consts["k3"], 232 K_SO2=equilibrium_consts["K_SO2"], 233 K_HSO3=equilibrium_consts["K_HSO3"], 234 dissociation_factor_SO2=dissociation_factors["SO2"], 235 timestep=timestep, 236 # input 237 droplet_volume=self.attributes["volume"], 238 pH=self.attributes["pH"], 239 # output 240 moles_O3=self.attributes["moles_O3"], 241 moles_H2O2=self.attributes["moles_H2O2"], 242 moles_S_IV=self.attributes["moles_S_IV"], 243 moles_S_VI=self.attributes["moles_S_VI"], 244 ) 245 for attr in ("moles_S_IV", "moles_S_VI", "moles_H2O2", "moles_O3"): 246 self.attributes.mark_updated(attr) 247 248 def dissolution( 249 self, 250 *, 251 gaseous_compounds, 252 system_type, 253 dissociation_factors, 254 timestep, 255 environment_mixing_ratios, 256 do_chemistry_flag, 257 ): 258 self.backend.dissolution( 259 n_cell=self.mesh.n_cell, 260 n_threads=1, 261 cell_order=np.arange(self.mesh.n_cell), 262 cell_start_arg=self.attributes.cell_start, 263 idx=self.attributes._ParticleAttributes__idx, 264 do_chemistry_flag=do_chemistry_flag, 265 mole_amounts={ 266 key: self.attributes["moles_" + key] for key in gaseous_compounds.keys() 267 }, 268 env_mixing_ratio=environment_mixing_ratios, 269 # note: assuming condensation was called 270 env_p=self.environment.get_predicted("p"), 271 env_T=self.environment.get_predicted("T"), 272 env_rho_d=self.environment.get_predicted("rhod"), 273 timestep=timestep, 274 dv=self.mesh.dv, 275 droplet_volume=self.attributes["volume"], 276 multiplicity=self.attributes["multiplicity"], 277 system_type=system_type, 278 dissociation_factors=dissociation_factors, 279 ) 280 for key in gaseous_compounds.keys(): 281 self.attributes.mark_updated(f"moles_{key}") 282 283 def chem_recalculate_cell_data(self, *, equilibrium_consts, kinetic_consts): 284 self.backend.chem_recalculate_cell_data( 285 equilibrium_consts=equilibrium_consts, 286 kinetic_consts=kinetic_consts, 287 temperature=self.environment.get_predicted("T"), 288 ) 289 290 def chem_recalculate_drop_data(self, *, dissociation_factors, equilibrium_consts): 291 self.backend.chem_recalculate_drop_data( 292 dissociation_factors=dissociation_factors, 293 equilibrium_consts=equilibrium_consts, 294 pH=self.attributes["pH"], 295 cell_id=self.attributes["cell id"], 296 ) 297 298 def recalculate_cell_id(self): 299 if not self.attributes.has_attribute("cell origin"): 300 return 301 self.backend.cell_id( 302 self.attributes["cell id"], 303 self.attributes["cell origin"], 304 self.backend.Storage.from_ndarray(self.environment.mesh.strides), 305 ) 306 self.attributes._ParticleAttributes__sorted = False 307 308 def sort_within_pair_by_attr(self, is_first_in_pair, attr_name): 309 self.backend.sort_within_pair_by_attr( 310 self.attributes._ParticleAttributes__idx, 311 is_first_in_pair, 312 self.attributes[attr_name], 313 ) 314 315 def moments( 316 self, 317 *, 318 moment_0, 319 moments, 320 specs: dict, 321 attr_name="signed water mass", 322 attr_range=(-np.inf, np.inf), 323 weighting_attribute="water mass", 324 weighting_rank=0, 325 skip_division_by_m0=False, 326 ): 327 """ 328 Writes to `moment_0` and `moment` the zero-th and the k-th statistical moments 329 of particle attributes computed filtering by value of the attribute `attr_name` 330 to fall within `attr_range`. The moment ranks are defined by `specs`. 331 332 Parameters: 333 specs: e.g., `specs={'volume': (1,2,3), 'kappa': (1)}` computes three moments 334 of volume and one moment of kappa 335 skip_division_by_m0: if set to `True`, the values written to `moments` are 336 multiplied by the 0-th moment (e.g., total volume instead of mean volume) 337 """ 338 if len(specs) == 0: 339 raise ValueError("empty specs passed") 340 attr_data, ranks = [], [] 341 for attr in specs: 342 for rank in specs[attr]: 343 attr_data.append(self.attributes[attr]) 344 ranks.append(rank) 345 assert len(set(attr_data)) <= 1 346 if len(attr_data) == 0: 347 attr_data = self.backend.Storage.empty((0,), dtype=float) 348 else: 349 attr_data = attr_data[0] 350 351 ranks = self.backend.Storage.from_ndarray(np.array(ranks, dtype=float)) 352 353 self.backend.moments( 354 moment_0=moment_0, 355 moments=moments, 356 multiplicity=self.attributes["multiplicity"], 357 attr_data=attr_data, 358 cell_id=self.attributes["cell id"], 359 idx=self.attributes._ParticleAttributes__idx, 360 length=self.attributes.super_droplet_count, 361 ranks=ranks, 362 min_x=attr_range[0], 363 max_x=attr_range[1], 364 x_attr=self.attributes[attr_name], 365 weighting_attribute=self.attributes[weighting_attribute], 366 weighting_rank=weighting_rank, 367 skip_division_by_m0=skip_division_by_m0, 368 ) 369 370 def spectrum_moments( 371 self, 372 *, 373 moment_0, 374 moments, 375 attr, 376 rank, 377 attr_bins, 378 attr_name="water mass", 379 weighting_attribute="water mass", 380 weighting_rank=0, 381 ): 382 attr_data = self.attributes[attr] 383 self.backend.spectrum_moments( 384 moment_0=moment_0, 385 moments=moments, 386 multiplicity=self.attributes["multiplicity"], 387 attr_data=attr_data, 388 cell_id=self.attributes["cell id"], 389 idx=self.attributes._ParticleAttributes__idx, 390 length=self.attributes.super_droplet_count, 391 rank=rank, 392 x_bins=attr_bins, 393 x_attr=self.attributes[attr_name], 394 weighting_attribute=self.attributes[weighting_attribute], 395 weighting_rank=weighting_rank, 396 ) 397 398 def adaptive_sdm_end(self, dt_left): 399 return self.backend.adaptive_sdm_end(dt_left, self.attributes.cell_start) 400 401 def remove_precipitated( 402 self, *, displacement, precipitation_counting_level_index 403 ) -> float: 404 rainfall_mass = self.backend.flag_precipitated( 405 cell_origin=self.attributes["cell origin"], 406 position_in_cell=self.attributes["position in cell"], 407 water_mass=self.attributes["water mass"], 408 multiplicity=self.attributes["multiplicity"], 409 idx=self.attributes._ParticleAttributes__idx, 410 length=self.attributes.super_droplet_count, 411 healthy=self.attributes._ParticleAttributes__healthy_memory, 412 precipitation_counting_level_index=precipitation_counting_level_index, 413 displacement=displacement, 414 ) 415 self.attributes.sanitize() 416 return rainfall_mass 417 418 def flag_out_of_column(self): 419 self.backend.flag_out_of_column( 420 cell_origin=self.attributes["cell origin"], 421 position_in_cell=self.attributes["position in cell"], 422 idx=self.attributes._ParticleAttributes__idx, 423 length=self.attributes.super_droplet_count, 424 healthy=self.attributes._ParticleAttributes__healthy_memory, 425 domain_top_level_index=self.mesh.grid[-1], 426 ) 427 self.attributes.sanitize() 428 429 def calculate_displacement( 430 self, *, displacement, courant, cell_origin, position_in_cell, n_substeps 431 ): 432 for dim in range(len(self.environment.mesh.grid)): 433 self.backend.calculate_displacement( 434 dim=dim, 435 displacement=displacement, 436 courant=courant[dim], 437 cell_origin=cell_origin, 438 position_in_cell=position_in_cell, 439 n_substeps=n_substeps, 440 ) 441 442 def isotopic_fractionation(self, heavy_isotopes: tuple): 443 self.backend.isotopic_fractionation() 444 for isotope in heavy_isotopes: 445 self.attributes.mark_updated(f"moles_{isotope}") 446 447 def seeding( 448 self, 449 *, 450 seeded_particle_index, 451 seeded_particle_multiplicity, 452 seeded_particle_extensive_attributes, 453 number_of_super_particles_to_inject, 454 ): 455 n_null = self.n_sd - self.attributes.super_droplet_count 456 if n_null == 0: 457 raise ValueError( 458 "No available seeds to inject. Please provide particles with nan filled attributes." 459 ) 460 461 if number_of_super_particles_to_inject > n_null: 462 raise ValueError( 463 "Trying to inject more super particles than space available." 464 ) 465 466 if number_of_super_particles_to_inject > len(seeded_particle_multiplicity): 467 raise ValueError( 468 "Trying to inject multiple super particles with the same attributes. \ 469 Instead increase multiplicity of injected particles." 470 ) 471 472 self.backend.seeding( 473 idx=self.attributes._ParticleAttributes__idx, 474 multiplicity=self.attributes["multiplicity"], 475 extensive_attributes=self.attributes.get_extensive_attribute_storage(), 476 seeded_particle_index=seeded_particle_index, 477 seeded_particle_multiplicity=seeded_particle_multiplicity, 478 seeded_particle_extensive_attributes=seeded_particle_extensive_attributes, 479 number_of_super_particles_to_inject=number_of_super_particles_to_inject, 480 ) 481 self.attributes.reset_idx() 482 self.attributes.sanitize() 483 484 self.attributes.mark_updated("multiplicity") 485 for key in self.attributes.get_extensive_attribute_keys(): 486 self.attributes.mark_updated(key) 487 488 def deposition(self): 489 self.backend.deposition( 490 multiplicity=self.attributes["multiplicity"], 491 signed_water_mass=self.attributes["signed water mass"], 492 current_temperature=self.environment["T"], 493 current_total_pressure=self.environment["p"], 494 current_relative_humidity=self.environment["RH"], 495 current_water_activity=self.environment["a_w_ice"], 496 current_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 497 current_dry_air_density=self.environment["rhod"], 498 current_dry_potential_temperature=self.environment["thd"], 499 cell_volume=self.environment.mesh.dv, 500 time_step=self.dt, 501 cell_id=self.attributes["cell id"], 502 reynolds_number=self.attributes["Reynolds number"], 503 schmidt_number=self.environment["Schmidt number"], 504 predicted_vapour_mixing_ratio=self.environment.get_predicted( 505 "water_vapour_mixing_ratio" 506 ), 507 predicted_dry_potential_temperature=self.environment.get_predicted("thd"), 508 ) 509 self.attributes.mark_updated("signed water mass") 510 # TODO #1524 - should we update here? 511 # self.update_TpRH(only_if_not_last='VapourDepositionOnIce') 512 513 def immersion_freezing_time_dependent(self, *, thaw: bool, rand: Storage): 514 self.backend.freeze_time_dependent( 515 rand=rand, 516 attributes=TimeDependentAttributes( 517 immersed_surface_area=self.attributes["immersed surface area"], 518 signed_water_mass=self.attributes["signed water mass"], 519 ), 520 timestep=self.dt, 521 cell=self.attributes["cell id"], 522 a_w_ice=self.environment["a_w_ice"], 523 temperature=self.environment["T"], 524 relative_humidity=self.environment["RH"], 525 thaw=thaw, 526 ) 527 self.attributes.mark_updated("signed water mass") 528 529 def immersion_freezing_singular(self, *, thaw: bool): 530 self.backend.freeze_singular( 531 attributes=SingularAttributes( 532 freezing_temperature=self.attributes["freezing temperature"], 533 signed_water_mass=self.attributes["signed water mass"], 534 ), 535 temperature=self.environment["T"], 536 relative_humidity=self.environment["RH"], 537 cell=self.attributes["cell id"], 538 thaw=thaw, 539 ) 540 self.attributes.mark_updated("signed water mass")
Particulator(n_sd, backend)
21 def __init__(self, n_sd, backend): 22 assert isinstance(backend, BackendMethods) 23 self.__n_sd = n_sd 24 25 self.backend = backend 26 self.formulae = backend.formulae 27 self.environment = None 28 self.attributes: (ParticleAttributes, None) = None 29 self.dynamics = {} 30 self.products = {} 31 self.observers = [] 32 33 self.n_steps = 0 34 35 self.sorting_scheme = "default" 36 self.condensation_solver = None 37 38 self.Index = make_Index(backend) # pylint: disable=invalid-name 39 self.PairIndicator = make_PairIndicator(backend) # pylint: disable=invalid-name 40 self.PairwiseStorage = make_PairwiseStorage( # pylint: disable=invalid-name 41 backend 42 ) 43 self.IndexedStorage = make_IndexedStorage( # pylint: disable=invalid-name 44 backend 45 ) 46 47 self.timers = {} 48 self.null = self.Storage.empty(0, dtype=float)
attributes: (<class 'PySDM.impl.particle_attributes.ParticleAttributes'>, None)
def
update_TpRH(self):
98 def update_TpRH(self): 99 self.backend.temperature_pressure_rh( 100 # input 101 rhod=self.environment.get_predicted("rhod"), 102 thd=self.environment.get_predicted("thd"), 103 water_vapour_mixing_ratio=self.environment.get_predicted( 104 "water_vapour_mixing_ratio" 105 ), 106 # output 107 T=self.environment.get_predicted("T"), 108 p=self.environment.get_predicted("p"), 109 RH=self.environment.get_predicted("RH"), 110 )
def
condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order):
112 def condensation(self, *, rtol_x, rtol_thd, counters, RH_max, success, cell_order): 113 """Updates droplet volumes by simulating condensation driven by prior changes 114 in environment thermodynamic state, updates the environment state. 115 In the case of parcel environment, condensation is driven solely by changes in 116 the dry-air density (theta and water_vapour_mixing_ratio should not be changed 117 by other dynamics). 118 In the case of prescribed-flow/kinematic environments, the dry-air density is 119 constant in time throughout the simulation. 120 This function should only change environment's predicted `thd` and 121 `water_vapour_mixing_ratio` (and not `rhod`). 122 """ 123 self.backend.condensation( 124 solver=self.condensation_solver, 125 n_cell=self.mesh.n_cell, 126 cell_start_arg=self.attributes.cell_start, 127 water_mass=self.attributes["signed water mass"], 128 multiplicity=self.attributes["multiplicity"], 129 vdry=self.attributes["dry volume"], 130 idx=self.attributes._ParticleAttributes__idx, 131 rhod=self.environment["rhod"], 132 thd=self.environment["thd"], 133 water_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 134 dv=self.environment.dv, 135 prhod=self.environment.get_predicted("rhod"), 136 pthd=self.environment.get_predicted("thd"), 137 predicted_water_vapour_mixing_ratio=self.environment.get_predicted( 138 "water_vapour_mixing_ratio" 139 ), 140 kappa=self.attributes["kappa"], 141 f_org=self.attributes["dry volume organic fraction"], 142 rtol_x=rtol_x, 143 rtol_thd=rtol_thd, 144 v_cr=self.attributes["critical volume"], 145 timestep=self.dt, 146 counters=counters, 147 cell_order=cell_order, 148 RH_max=RH_max, 149 success=success, 150 cell_id=self.attributes["cell id"], 151 reynolds_number=self.attributes["Reynolds number"], 152 air_density=self.environment["air density"], 153 air_dynamic_viscosity=self.environment["air dynamic viscosity"], 154 ) 155 self.attributes.mark_updated("signed water mass")
Updates droplet volumes by simulating condensation driven by prior changes
in environment thermodynamic state, updates the environment state.
In the case of parcel environment, condensation is driven solely by changes in
the dry-air density (theta and water_vapour_mixing_ratio should not be changed
by other dynamics).
In the case of prescribed-flow/kinematic environments, the dry-air density is
constant in time throughout the simulation.
This function should only change environment's predicted thd and
water_vapour_mixing_ratio (and not rhod).
def
collision_coalescence_breakup( self, *, enable_breakup, gamma, rand, Ec, Eb, fragment_mass, coalescence_rate, breakup_rate, breakup_rate_deficit, is_first_in_pair, warn_overflows, max_multiplicity):
157 def collision_coalescence_breakup( 158 self, 159 *, 160 enable_breakup, 161 gamma, 162 rand, 163 Ec, 164 Eb, 165 fragment_mass, 166 coalescence_rate, 167 breakup_rate, 168 breakup_rate_deficit, 169 is_first_in_pair, 170 warn_overflows, 171 max_multiplicity, 172 ): 173 # pylint: disable=too-many-locals 174 idx = self.attributes._ParticleAttributes__idx 175 healthy = self.attributes._ParticleAttributes__healthy_memory 176 cell_id = self.attributes["cell id"] 177 multiplicity = self.attributes["multiplicity"] 178 attributes = self.attributes.get_extensive_attribute_storage() 179 if enable_breakup: 180 self.backend.collision_coalescence_breakup( 181 multiplicity=multiplicity, 182 idx=idx, 183 attributes=attributes, 184 gamma=gamma, 185 rand=rand, 186 Ec=Ec, 187 Eb=Eb, 188 fragment_mass=fragment_mass, 189 healthy=healthy, 190 cell_id=cell_id, 191 coalescence_rate=coalescence_rate, 192 breakup_rate=breakup_rate, 193 breakup_rate_deficit=breakup_rate_deficit, 194 is_first_in_pair=is_first_in_pair, 195 warn_overflows=warn_overflows, 196 particle_mass=self.attributes["water mass"], 197 max_multiplicity=max_multiplicity, 198 ) 199 else: 200 self.backend.collision_coalescence( 201 multiplicity=multiplicity, 202 idx=idx, 203 attributes=attributes, 204 gamma=gamma, 205 healthy=healthy, 206 cell_id=cell_id, 207 coalescence_rate=coalescence_rate, 208 is_first_in_pair=is_first_in_pair, 209 ) 210 self.attributes.sanitize() 211 self.attributes.mark_updated("multiplicity") 212 for key in self.attributes.get_extensive_attribute_keys(): 213 self.attributes.mark_updated(key)
def
oxidation( self, *, kinetic_consts, timestep, equilibrium_consts, dissociation_factors, do_chemistry_flag):
215 def oxidation( 216 self, 217 *, 218 kinetic_consts, 219 timestep, 220 equilibrium_consts, 221 dissociation_factors, 222 do_chemistry_flag, 223 ): 224 self.backend.oxidation( 225 n_sd=self.n_sd, 226 cell_ids=self.attributes["cell id"], 227 do_chemistry_flag=do_chemistry_flag, 228 k0=kinetic_consts["k0"], 229 k1=kinetic_consts["k1"], 230 k2=kinetic_consts["k2"], 231 k3=kinetic_consts["k3"], 232 K_SO2=equilibrium_consts["K_SO2"], 233 K_HSO3=equilibrium_consts["K_HSO3"], 234 dissociation_factor_SO2=dissociation_factors["SO2"], 235 timestep=timestep, 236 # input 237 droplet_volume=self.attributes["volume"], 238 pH=self.attributes["pH"], 239 # output 240 moles_O3=self.attributes["moles_O3"], 241 moles_H2O2=self.attributes["moles_H2O2"], 242 moles_S_IV=self.attributes["moles_S_IV"], 243 moles_S_VI=self.attributes["moles_S_VI"], 244 ) 245 for attr in ("moles_S_IV", "moles_S_VI", "moles_H2O2", "moles_O3"): 246 self.attributes.mark_updated(attr)
def
dissolution( self, *, gaseous_compounds, system_type, dissociation_factors, timestep, environment_mixing_ratios, do_chemistry_flag):
248 def dissolution( 249 self, 250 *, 251 gaseous_compounds, 252 system_type, 253 dissociation_factors, 254 timestep, 255 environment_mixing_ratios, 256 do_chemistry_flag, 257 ): 258 self.backend.dissolution( 259 n_cell=self.mesh.n_cell, 260 n_threads=1, 261 cell_order=np.arange(self.mesh.n_cell), 262 cell_start_arg=self.attributes.cell_start, 263 idx=self.attributes._ParticleAttributes__idx, 264 do_chemistry_flag=do_chemistry_flag, 265 mole_amounts={ 266 key: self.attributes["moles_" + key] for key in gaseous_compounds.keys() 267 }, 268 env_mixing_ratio=environment_mixing_ratios, 269 # note: assuming condensation was called 270 env_p=self.environment.get_predicted("p"), 271 env_T=self.environment.get_predicted("T"), 272 env_rho_d=self.environment.get_predicted("rhod"), 273 timestep=timestep, 274 dv=self.mesh.dv, 275 droplet_volume=self.attributes["volume"], 276 multiplicity=self.attributes["multiplicity"], 277 system_type=system_type, 278 dissociation_factors=dissociation_factors, 279 ) 280 for key in gaseous_compounds.keys(): 281 self.attributes.mark_updated(f"moles_{key}")
def
recalculate_cell_id(self):
298 def recalculate_cell_id(self): 299 if not self.attributes.has_attribute("cell origin"): 300 return 301 self.backend.cell_id( 302 self.attributes["cell id"], 303 self.attributes["cell origin"], 304 self.backend.Storage.from_ndarray(self.environment.mesh.strides), 305 ) 306 self.attributes._ParticleAttributes__sorted = False
def
moments( self, *, moment_0, moments, specs: dict, attr_name='signed water mass', attr_range=(-inf, inf), weighting_attribute='water mass', weighting_rank=0, skip_division_by_m0=False):
315 def moments( 316 self, 317 *, 318 moment_0, 319 moments, 320 specs: dict, 321 attr_name="signed water mass", 322 attr_range=(-np.inf, np.inf), 323 weighting_attribute="water mass", 324 weighting_rank=0, 325 skip_division_by_m0=False, 326 ): 327 """ 328 Writes to `moment_0` and `moment` the zero-th and the k-th statistical moments 329 of particle attributes computed filtering by value of the attribute `attr_name` 330 to fall within `attr_range`. The moment ranks are defined by `specs`. 331 332 Parameters: 333 specs: e.g., `specs={'volume': (1,2,3), 'kappa': (1)}` computes three moments 334 of volume and one moment of kappa 335 skip_division_by_m0: if set to `True`, the values written to `moments` are 336 multiplied by the 0-th moment (e.g., total volume instead of mean volume) 337 """ 338 if len(specs) == 0: 339 raise ValueError("empty specs passed") 340 attr_data, ranks = [], [] 341 for attr in specs: 342 for rank in specs[attr]: 343 attr_data.append(self.attributes[attr]) 344 ranks.append(rank) 345 assert len(set(attr_data)) <= 1 346 if len(attr_data) == 0: 347 attr_data = self.backend.Storage.empty((0,), dtype=float) 348 else: 349 attr_data = attr_data[0] 350 351 ranks = self.backend.Storage.from_ndarray(np.array(ranks, dtype=float)) 352 353 self.backend.moments( 354 moment_0=moment_0, 355 moments=moments, 356 multiplicity=self.attributes["multiplicity"], 357 attr_data=attr_data, 358 cell_id=self.attributes["cell id"], 359 idx=self.attributes._ParticleAttributes__idx, 360 length=self.attributes.super_droplet_count, 361 ranks=ranks, 362 min_x=attr_range[0], 363 max_x=attr_range[1], 364 x_attr=self.attributes[attr_name], 365 weighting_attribute=self.attributes[weighting_attribute], 366 weighting_rank=weighting_rank, 367 skip_division_by_m0=skip_division_by_m0, 368 )
Writes to moment_0 and moment the zero-th and the k-th statistical moments
of particle attributes computed filtering by value of the attribute attr_name
to fall within attr_range. The moment ranks are defined by specs.
Parameters:
specs: e.g., specs={'volume': (1,2,3), 'kappa': (1)} computes three moments
of volume and one moment of kappa
skip_division_by_m0: if set to True, the values written to moments are
multiplied by the 0-th moment (e.g., total volume instead of mean volume)
def
spectrum_moments( self, *, moment_0, moments, attr, rank, attr_bins, attr_name='water mass', weighting_attribute='water mass', weighting_rank=0):
370 def spectrum_moments( 371 self, 372 *, 373 moment_0, 374 moments, 375 attr, 376 rank, 377 attr_bins, 378 attr_name="water mass", 379 weighting_attribute="water mass", 380 weighting_rank=0, 381 ): 382 attr_data = self.attributes[attr] 383 self.backend.spectrum_moments( 384 moment_0=moment_0, 385 moments=moments, 386 multiplicity=self.attributes["multiplicity"], 387 attr_data=attr_data, 388 cell_id=self.attributes["cell id"], 389 idx=self.attributes._ParticleAttributes__idx, 390 length=self.attributes.super_droplet_count, 391 rank=rank, 392 x_bins=attr_bins, 393 x_attr=self.attributes[attr_name], 394 weighting_attribute=self.attributes[weighting_attribute], 395 weighting_rank=weighting_rank, 396 )
def
remove_precipitated(self, *, displacement, precipitation_counting_level_index) -> float:
401 def remove_precipitated( 402 self, *, displacement, precipitation_counting_level_index 403 ) -> float: 404 rainfall_mass = self.backend.flag_precipitated( 405 cell_origin=self.attributes["cell origin"], 406 position_in_cell=self.attributes["position in cell"], 407 water_mass=self.attributes["water mass"], 408 multiplicity=self.attributes["multiplicity"], 409 idx=self.attributes._ParticleAttributes__idx, 410 length=self.attributes.super_droplet_count, 411 healthy=self.attributes._ParticleAttributes__healthy_memory, 412 precipitation_counting_level_index=precipitation_counting_level_index, 413 displacement=displacement, 414 ) 415 self.attributes.sanitize() 416 return rainfall_mass
def
flag_out_of_column(self):
418 def flag_out_of_column(self): 419 self.backend.flag_out_of_column( 420 cell_origin=self.attributes["cell origin"], 421 position_in_cell=self.attributes["position in cell"], 422 idx=self.attributes._ParticleAttributes__idx, 423 length=self.attributes.super_droplet_count, 424 healthy=self.attributes._ParticleAttributes__healthy_memory, 425 domain_top_level_index=self.mesh.grid[-1], 426 ) 427 self.attributes.sanitize()
def
calculate_displacement( self, *, displacement, courant, cell_origin, position_in_cell, n_substeps):
429 def calculate_displacement( 430 self, *, displacement, courant, cell_origin, position_in_cell, n_substeps 431 ): 432 for dim in range(len(self.environment.mesh.grid)): 433 self.backend.calculate_displacement( 434 dim=dim, 435 displacement=displacement, 436 courant=courant[dim], 437 cell_origin=cell_origin, 438 position_in_cell=position_in_cell, 439 n_substeps=n_substeps, 440 )
def
seeding( self, *, seeded_particle_index, seeded_particle_multiplicity, seeded_particle_extensive_attributes, number_of_super_particles_to_inject):
447 def seeding( 448 self, 449 *, 450 seeded_particle_index, 451 seeded_particle_multiplicity, 452 seeded_particle_extensive_attributes, 453 number_of_super_particles_to_inject, 454 ): 455 n_null = self.n_sd - self.attributes.super_droplet_count 456 if n_null == 0: 457 raise ValueError( 458 "No available seeds to inject. Please provide particles with nan filled attributes." 459 ) 460 461 if number_of_super_particles_to_inject > n_null: 462 raise ValueError( 463 "Trying to inject more super particles than space available." 464 ) 465 466 if number_of_super_particles_to_inject > len(seeded_particle_multiplicity): 467 raise ValueError( 468 "Trying to inject multiple super particles with the same attributes. \ 469 Instead increase multiplicity of injected particles." 470 ) 471 472 self.backend.seeding( 473 idx=self.attributes._ParticleAttributes__idx, 474 multiplicity=self.attributes["multiplicity"], 475 extensive_attributes=self.attributes.get_extensive_attribute_storage(), 476 seeded_particle_index=seeded_particle_index, 477 seeded_particle_multiplicity=seeded_particle_multiplicity, 478 seeded_particle_extensive_attributes=seeded_particle_extensive_attributes, 479 number_of_super_particles_to_inject=number_of_super_particles_to_inject, 480 ) 481 self.attributes.reset_idx() 482 self.attributes.sanitize() 483 484 self.attributes.mark_updated("multiplicity") 485 for key in self.attributes.get_extensive_attribute_keys(): 486 self.attributes.mark_updated(key)
def
deposition(self):
488 def deposition(self): 489 self.backend.deposition( 490 multiplicity=self.attributes["multiplicity"], 491 signed_water_mass=self.attributes["signed water mass"], 492 current_temperature=self.environment["T"], 493 current_total_pressure=self.environment["p"], 494 current_relative_humidity=self.environment["RH"], 495 current_water_activity=self.environment["a_w_ice"], 496 current_vapour_mixing_ratio=self.environment["water_vapour_mixing_ratio"], 497 current_dry_air_density=self.environment["rhod"], 498 current_dry_potential_temperature=self.environment["thd"], 499 cell_volume=self.environment.mesh.dv, 500 time_step=self.dt, 501 cell_id=self.attributes["cell id"], 502 reynolds_number=self.attributes["Reynolds number"], 503 schmidt_number=self.environment["Schmidt number"], 504 predicted_vapour_mixing_ratio=self.environment.get_predicted( 505 "water_vapour_mixing_ratio" 506 ), 507 predicted_dry_potential_temperature=self.environment.get_predicted("thd"), 508 ) 509 self.attributes.mark_updated("signed water mass") 510 # TODO #1524 - should we update here? 511 # self.update_TpRH(only_if_not_last='VapourDepositionOnIce')
def
immersion_freezing_time_dependent(self, *, thaw: bool, rand: <property object>):
513 def immersion_freezing_time_dependent(self, *, thaw: bool, rand: Storage): 514 self.backend.freeze_time_dependent( 515 rand=rand, 516 attributes=TimeDependentAttributes( 517 immersed_surface_area=self.attributes["immersed surface area"], 518 signed_water_mass=self.attributes["signed water mass"], 519 ), 520 timestep=self.dt, 521 cell=self.attributes["cell id"], 522 a_w_ice=self.environment["a_w_ice"], 523 temperature=self.environment["T"], 524 relative_humidity=self.environment["RH"], 525 thaw=thaw, 526 ) 527 self.attributes.mark_updated("signed water mass")
def
immersion_freezing_singular(self, *, thaw: bool):
529 def immersion_freezing_singular(self, *, thaw: bool): 530 self.backend.freeze_singular( 531 attributes=SingularAttributes( 532 freezing_temperature=self.attributes["freezing temperature"], 533 signed_water_mass=self.attributes["signed water mass"], 534 ), 535 temperature=self.environment["T"], 536 relative_humidity=self.environment["RH"], 537 cell=self.attributes["cell id"], 538 thaw=thaw, 539 ) 540 self.attributes.mark_updated("signed water mass")