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