PySDM_examples.sandbox.seeding.simulation
1import numpy as np 2 3from PySDM_examples.sandbox.seeding.settings import Settings 4 5from PySDM import Particulator 6from PySDM.backends import CPU 7from PySDM.environments import Parcel 8from PySDM.dynamics import Condensation, AmbientThermodynamics, Coalescence, Seeding 9from PySDM.dynamics.collisions.collision_kernels import Geometric 10from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity 11from PySDM import products 12from PySDM.physics import si 13 14 15class Simulation: 16 def __init__(self, settings: Settings): 17 environment = Parcel( 18 dt=settings.timestep, 19 mass_of_dry_air=settings.mass_of_dry_air, 20 w=settings.updraft, 21 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 22 p0=settings.initial_total_pressure, 23 T0=settings.initial_temperature, 24 backend=CPU(formulae=settings.formulae), 25 ) 26 r_dry, n_in_dv = ConstantMultiplicity( 27 settings.initial_aerosol_dry_radii 28 ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend) 29 attributes = environment.init_attributes( 30 n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry 31 ) 32 self.particulator = Particulator( 33 environment=environment, 34 n_sd=settings.n_sd_seeding + settings.n_sd_initial, 35 dynamics=[ 36 AmbientThermodynamics(), 37 Condensation(), 38 ] 39 + ( 40 [] 41 if not settings.enable_collisions 42 else [ 43 Coalescence(collision_kernel=Geometric()), 44 ] 45 ) 46 + [ 47 Seeding( 48 **{ 49 k: getattr(settings, k) 50 for k in ( 51 "super_droplet_injection_rate", 52 "seeded_particle_multiplicity", 53 "seeded_particle_extensive_attributes", 54 ) 55 } 56 ), 57 ], 58 attributes={ 59 k: np.pad( 60 array=v, 61 pad_width=(0, settings.n_sd_seeding), 62 mode="constant", 63 constant_values=np.nan if k == "multiplicity" else 0, 64 ) 65 for k, v in attributes.items() 66 }, 67 products=( 68 products.SuperDropletCountPerGridbox(name="sd_count"), 69 products.Time(), 70 products.WaterMixingRatio( 71 radius_range=(settings.rain_water_radius_threshold, np.inf), 72 name="rain water mixing ratio", 73 ), 74 products.EffectiveRadius( 75 name="r_eff", 76 unit="um", 77 radius_range=(0.5 * si.um, 25 * si.um), 78 ), 79 products.ParticleConcentration( 80 name="n_drop", 81 unit="cm^-3", 82 radius_range=(0.5 * si.um, 25 * si.um), 83 ), 84 ), 85 ) 86 self.n_steps = int(settings.t_max // settings.timestep) 87 88 def run(self): 89 output = { 90 "attributes": {"water mass": []}, 91 "products": {key: [] for key in self.particulator.products}, 92 } 93 for step in range(self.n_steps + 1): 94 if step != 0: 95 self.particulator.advance(n_steps=1) 96 for key, attr in output["attributes"].items(): 97 data = self.particulator.attributes[key].to_ndarray(raw=True) 98 data[data == 0] = np.nan 99 attr.append(data) 100 for key, prod in output["products"].items(): 101 value = self.particulator.products[key].get() 102 if not isinstance(value, float): 103 (value,) = value 104 prod.append(float(value)) 105 for out in ("attributes", "products"): 106 for key, val in output[out].items(): 107 output[out][key] = np.array(val) 108 return output
class
Simulation:
16class Simulation: 17 def __init__(self, settings: Settings): 18 environment = Parcel( 19 dt=settings.timestep, 20 mass_of_dry_air=settings.mass_of_dry_air, 21 w=settings.updraft, 22 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 23 p0=settings.initial_total_pressure, 24 T0=settings.initial_temperature, 25 backend=CPU(formulae=settings.formulae), 26 ) 27 r_dry, n_in_dv = ConstantMultiplicity( 28 settings.initial_aerosol_dry_radii 29 ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend) 30 attributes = environment.init_attributes( 31 n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry 32 ) 33 self.particulator = Particulator( 34 environment=environment, 35 n_sd=settings.n_sd_seeding + settings.n_sd_initial, 36 dynamics=[ 37 AmbientThermodynamics(), 38 Condensation(), 39 ] 40 + ( 41 [] 42 if not settings.enable_collisions 43 else [ 44 Coalescence(collision_kernel=Geometric()), 45 ] 46 ) 47 + [ 48 Seeding( 49 **{ 50 k: getattr(settings, k) 51 for k in ( 52 "super_droplet_injection_rate", 53 "seeded_particle_multiplicity", 54 "seeded_particle_extensive_attributes", 55 ) 56 } 57 ), 58 ], 59 attributes={ 60 k: np.pad( 61 array=v, 62 pad_width=(0, settings.n_sd_seeding), 63 mode="constant", 64 constant_values=np.nan if k == "multiplicity" else 0, 65 ) 66 for k, v in attributes.items() 67 }, 68 products=( 69 products.SuperDropletCountPerGridbox(name="sd_count"), 70 products.Time(), 71 products.WaterMixingRatio( 72 radius_range=(settings.rain_water_radius_threshold, np.inf), 73 name="rain water mixing ratio", 74 ), 75 products.EffectiveRadius( 76 name="r_eff", 77 unit="um", 78 radius_range=(0.5 * si.um, 25 * si.um), 79 ), 80 products.ParticleConcentration( 81 name="n_drop", 82 unit="cm^-3", 83 radius_range=(0.5 * si.um, 25 * si.um), 84 ), 85 ), 86 ) 87 self.n_steps = int(settings.t_max // settings.timestep) 88 89 def run(self): 90 output = { 91 "attributes": {"water mass": []}, 92 "products": {key: [] for key in self.particulator.products}, 93 } 94 for step in range(self.n_steps + 1): 95 if step != 0: 96 self.particulator.advance(n_steps=1) 97 for key, attr in output["attributes"].items(): 98 data = self.particulator.attributes[key].to_ndarray(raw=True) 99 data[data == 0] = np.nan 100 attr.append(data) 101 for key, prod in output["products"].items(): 102 value = self.particulator.products[key].get() 103 if not isinstance(value, float): 104 (value,) = value 105 prod.append(float(value)) 106 for out in ("attributes", "products"): 107 for key, val in output[out].items(): 108 output[out][key] = np.array(val) 109 return output
Simulation(settings: PySDM_examples.sandbox.seeding.settings.Settings)
17 def __init__(self, settings: Settings): 18 environment = Parcel( 19 dt=settings.timestep, 20 mass_of_dry_air=settings.mass_of_dry_air, 21 w=settings.updraft, 22 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 23 p0=settings.initial_total_pressure, 24 T0=settings.initial_temperature, 25 backend=CPU(formulae=settings.formulae), 26 ) 27 r_dry, n_in_dv = ConstantMultiplicity( 28 settings.initial_aerosol_dry_radii 29 ).sample_deterministic(n_sd=settings.n_sd_initial, backend=environment.backend) 30 attributes = environment.init_attributes( 31 n_in_dv=n_in_dv, kappa=settings.initial_aerosol_kappa, r_dry=r_dry 32 ) 33 self.particulator = Particulator( 34 environment=environment, 35 n_sd=settings.n_sd_seeding + settings.n_sd_initial, 36 dynamics=[ 37 AmbientThermodynamics(), 38 Condensation(), 39 ] 40 + ( 41 [] 42 if not settings.enable_collisions 43 else [ 44 Coalescence(collision_kernel=Geometric()), 45 ] 46 ) 47 + [ 48 Seeding( 49 **{ 50 k: getattr(settings, k) 51 for k in ( 52 "super_droplet_injection_rate", 53 "seeded_particle_multiplicity", 54 "seeded_particle_extensive_attributes", 55 ) 56 } 57 ), 58 ], 59 attributes={ 60 k: np.pad( 61 array=v, 62 pad_width=(0, settings.n_sd_seeding), 63 mode="constant", 64 constant_values=np.nan if k == "multiplicity" else 0, 65 ) 66 for k, v in attributes.items() 67 }, 68 products=( 69 products.SuperDropletCountPerGridbox(name="sd_count"), 70 products.Time(), 71 products.WaterMixingRatio( 72 radius_range=(settings.rain_water_radius_threshold, np.inf), 73 name="rain water mixing ratio", 74 ), 75 products.EffectiveRadius( 76 name="r_eff", 77 unit="um", 78 radius_range=(0.5 * si.um, 25 * si.um), 79 ), 80 products.ParticleConcentration( 81 name="n_drop", 82 unit="cm^-3", 83 radius_range=(0.5 * si.um, 25 * si.um), 84 ), 85 ), 86 ) 87 self.n_steps = int(settings.t_max // settings.timestep)
def
run(self):
89 def run(self): 90 output = { 91 "attributes": {"water mass": []}, 92 "products": {key: [] for key in self.particulator.products}, 93 } 94 for step in range(self.n_steps + 1): 95 if step != 0: 96 self.particulator.advance(n_steps=1) 97 for key, attr in output["attributes"].items(): 98 data = self.particulator.attributes[key].to_ndarray(raw=True) 99 data[data == 0] = np.nan 100 attr.append(data) 101 for key, prod in output["products"].items(): 102 value = self.particulator.products[key].get() 103 if not isinstance(value, float): 104 (value,) = value 105 prod.append(float(value)) 106 for out in ("attributes", "products"): 107 for key, val in output[out].items(): 108 output[out][key] = np.array(val) 109 return output