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