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