PySDM_examples.featured.Jensen_and_Nugent_2017.simulation
1import numpy as np 2from PySDM_examples.utils import BasicSimulation 3from PySDM_examples.featured.Jensen_and_Nugent_2017.settings import Settings 4from PySDM_examples.featured.Jensen_and_Nugent_2017 import table_3 5from PySDM import Particulator 6from PySDM.physics import si 7from PySDM.backends import CPU 8from PySDM.products import ( 9 PeakSaturation, 10 ParcelDisplacement, 11 Time, 12 ActivatedMeanRadius, 13 RadiusStandardDeviation, 14) 15from PySDM.environments import Parcel 16from PySDM.dynamics import Condensation, AmbientThermodynamics, Coalescence 17from PySDM.dynamics.collisions.collision_kernels import Geometric 18from PySDM.initialisation.sampling.spectral_sampling import Logarithmic 19 20# note: 100 in caption of Table 1 21N_SD_NON_GCCN = 100 22 23 24class Simulation(BasicSimulation): 25 def __init__( 26 self, 27 settings: Settings, 28 gccn: bool = False, 29 gravitational_coalsecence: bool = False, 30 ): 31 32 n_gccn = np.count_nonzero(table_3.NA) if gccn else 0 33 environment = Parcel( 34 dt=settings.dt, 35 mass_of_dry_air=666 * si.kg, 36 p0=settings.p0, 37 initial_relative_humidity=settings.RH0, 38 T0=settings.T0, 39 w=settings.vertical_velocity, 40 z0=settings.z0, 41 backend=CPU(formulae=settings.formulae), 42 ) 43 44 self.r_dry, n_in_unit_volume = Logarithmic( 45 spectrum=settings.dry_radii_spectrum, 46 ).sample_deterministic(N_SD_NON_GCCN) 47 48 if gccn: 49 nonzero_concentration_mask = np.nonzero(table_3.NA) 50 self.r_dry = np.concatenate( 51 [self.r_dry, table_3.RD[nonzero_concentration_mask]] 52 ) 53 n_in_unit_volume = np.concatenate( 54 [n_in_unit_volume, table_3.NA[nonzero_concentration_mask]] 55 ) # TODO #1266: check which temp, pres, RH assumed in the paper for NA??? 56 57 pd0 = settings.formulae.trivia.p_d( 58 settings.p0, 59 settings.formulae.trivia.water_vapour_mixing_ratio( 60 settings.p0, 61 settings.RH0, 62 settings.formulae.saturation_vapour_pressure.pvs_water(settings.T0), 63 ), 64 ) 65 rhod0 = settings.formulae.state_variable_triplet.rhod_of_pd_T(pd0, settings.T0) 66 67 attributes = environment.init_attributes( 68 n_in_dv=n_in_unit_volume * environment.mass_of_dry_air / rhod0, 69 kappa=settings.kappa, 70 r_dry=self.r_dry, 71 ) 72 73 super().__init__( 74 Particulator( 75 n_sd=N_SD_NON_GCCN + n_gccn, 76 environment=environment, 77 dynamics=[ 78 # TODO #1266: order matters here, but error message is not saying it! 79 AmbientThermodynamics(), 80 Condensation(), 81 ] 82 + ( 83 [] 84 if not gravitational_coalsecence 85 else [Coalescence(collision_kernel=Geometric())] 86 ), 87 attributes=attributes, 88 products=( 89 PeakSaturation(name="S_max"), 90 ParcelDisplacement(name="z"), 91 Time(name="t"), 92 ActivatedMeanRadius( 93 name="r_mean_act", count_activated=True, count_unactivated=False 94 ), 95 RadiusStandardDeviation( 96 name="r_std_act", count_activated=True, count_unactivated=False 97 ), 98 ), 99 requested_attributes=( 100 additional_derived_attributes := ( 101 "radius", 102 "equilibrium saturation", 103 ) 104 ), 105 ) 106 ) 107 108 # TODO #1266: copied from G & P 2023 109 self.output_attributes = { 110 attr: tuple([] for _ in range(self.particulator.n_sd)) 111 for attr in additional_derived_attributes 112 } 113 114 def run( 115 self, *, n_steps: int = 2250, steps_per_output_interval: int = 10 116 ): # TODO #1266: essentially copied from G & P 2023 117 output_products = super()._run( 118 nt=n_steps, steps_per_output_interval=steps_per_output_interval 119 ) 120 return {"products": output_products, "attributes": self.output_attributes} 121 122 def _save(self, output): # TODO #1266: copied from G&P 2023 123 for key, attr in self.output_attributes.items(): 124 attr_data = self.particulator.attributes[key].to_ndarray() 125 for drop_id in range(self.particulator.n_sd): 126 attr[drop_id].append(attr_data[drop_id]) 127 super()._save(output)
N_SD_NON_GCCN =
100
25class Simulation(BasicSimulation): 26 def __init__( 27 self, 28 settings: Settings, 29 gccn: bool = False, 30 gravitational_coalsecence: bool = False, 31 ): 32 33 n_gccn = np.count_nonzero(table_3.NA) if gccn else 0 34 environment = Parcel( 35 dt=settings.dt, 36 mass_of_dry_air=666 * si.kg, 37 p0=settings.p0, 38 initial_relative_humidity=settings.RH0, 39 T0=settings.T0, 40 w=settings.vertical_velocity, 41 z0=settings.z0, 42 backend=CPU(formulae=settings.formulae), 43 ) 44 45 self.r_dry, n_in_unit_volume = Logarithmic( 46 spectrum=settings.dry_radii_spectrum, 47 ).sample_deterministic(N_SD_NON_GCCN) 48 49 if gccn: 50 nonzero_concentration_mask = np.nonzero(table_3.NA) 51 self.r_dry = np.concatenate( 52 [self.r_dry, table_3.RD[nonzero_concentration_mask]] 53 ) 54 n_in_unit_volume = np.concatenate( 55 [n_in_unit_volume, table_3.NA[nonzero_concentration_mask]] 56 ) # TODO #1266: check which temp, pres, RH assumed in the paper for NA??? 57 58 pd0 = settings.formulae.trivia.p_d( 59 settings.p0, 60 settings.formulae.trivia.water_vapour_mixing_ratio( 61 settings.p0, 62 settings.RH0, 63 settings.formulae.saturation_vapour_pressure.pvs_water(settings.T0), 64 ), 65 ) 66 rhod0 = settings.formulae.state_variable_triplet.rhod_of_pd_T(pd0, settings.T0) 67 68 attributes = environment.init_attributes( 69 n_in_dv=n_in_unit_volume * environment.mass_of_dry_air / rhod0, 70 kappa=settings.kappa, 71 r_dry=self.r_dry, 72 ) 73 74 super().__init__( 75 Particulator( 76 n_sd=N_SD_NON_GCCN + n_gccn, 77 environment=environment, 78 dynamics=[ 79 # TODO #1266: order matters here, but error message is not saying it! 80 AmbientThermodynamics(), 81 Condensation(), 82 ] 83 + ( 84 [] 85 if not gravitational_coalsecence 86 else [Coalescence(collision_kernel=Geometric())] 87 ), 88 attributes=attributes, 89 products=( 90 PeakSaturation(name="S_max"), 91 ParcelDisplacement(name="z"), 92 Time(name="t"), 93 ActivatedMeanRadius( 94 name="r_mean_act", count_activated=True, count_unactivated=False 95 ), 96 RadiusStandardDeviation( 97 name="r_std_act", count_activated=True, count_unactivated=False 98 ), 99 ), 100 requested_attributes=( 101 additional_derived_attributes := ( 102 "radius", 103 "equilibrium saturation", 104 ) 105 ), 106 ) 107 ) 108 109 # TODO #1266: copied from G & P 2023 110 self.output_attributes = { 111 attr: tuple([] for _ in range(self.particulator.n_sd)) 112 for attr in additional_derived_attributes 113 } 114 115 def run( 116 self, *, n_steps: int = 2250, steps_per_output_interval: int = 10 117 ): # TODO #1266: essentially copied from G & P 2023 118 output_products = super()._run( 119 nt=n_steps, steps_per_output_interval=steps_per_output_interval 120 ) 121 return {"products": output_products, "attributes": self.output_attributes} 122 123 def _save(self, output): # TODO #1266: copied from G&P 2023 124 for key, attr in self.output_attributes.items(): 125 attr_data = self.particulator.attributes[key].to_ndarray() 126 for drop_id in range(self.particulator.n_sd): 127 attr[drop_id].append(attr_data[drop_id]) 128 super()._save(output)
Simulation( settings: PySDM_examples.featured.Jensen_and_Nugent_2017.settings.Settings, gccn: bool = False, gravitational_coalsecence: bool = False)
26 def __init__( 27 self, 28 settings: Settings, 29 gccn: bool = False, 30 gravitational_coalsecence: bool = False, 31 ): 32 33 n_gccn = np.count_nonzero(table_3.NA) if gccn else 0 34 environment = Parcel( 35 dt=settings.dt, 36 mass_of_dry_air=666 * si.kg, 37 p0=settings.p0, 38 initial_relative_humidity=settings.RH0, 39 T0=settings.T0, 40 w=settings.vertical_velocity, 41 z0=settings.z0, 42 backend=CPU(formulae=settings.formulae), 43 ) 44 45 self.r_dry, n_in_unit_volume = Logarithmic( 46 spectrum=settings.dry_radii_spectrum, 47 ).sample_deterministic(N_SD_NON_GCCN) 48 49 if gccn: 50 nonzero_concentration_mask = np.nonzero(table_3.NA) 51 self.r_dry = np.concatenate( 52 [self.r_dry, table_3.RD[nonzero_concentration_mask]] 53 ) 54 n_in_unit_volume = np.concatenate( 55 [n_in_unit_volume, table_3.NA[nonzero_concentration_mask]] 56 ) # TODO #1266: check which temp, pres, RH assumed in the paper for NA??? 57 58 pd0 = settings.formulae.trivia.p_d( 59 settings.p0, 60 settings.formulae.trivia.water_vapour_mixing_ratio( 61 settings.p0, 62 settings.RH0, 63 settings.formulae.saturation_vapour_pressure.pvs_water(settings.T0), 64 ), 65 ) 66 rhod0 = settings.formulae.state_variable_triplet.rhod_of_pd_T(pd0, settings.T0) 67 68 attributes = environment.init_attributes( 69 n_in_dv=n_in_unit_volume * environment.mass_of_dry_air / rhod0, 70 kappa=settings.kappa, 71 r_dry=self.r_dry, 72 ) 73 74 super().__init__( 75 Particulator( 76 n_sd=N_SD_NON_GCCN + n_gccn, 77 environment=environment, 78 dynamics=[ 79 # TODO #1266: order matters here, but error message is not saying it! 80 AmbientThermodynamics(), 81 Condensation(), 82 ] 83 + ( 84 [] 85 if not gravitational_coalsecence 86 else [Coalescence(collision_kernel=Geometric())] 87 ), 88 attributes=attributes, 89 products=( 90 PeakSaturation(name="S_max"), 91 ParcelDisplacement(name="z"), 92 Time(name="t"), 93 ActivatedMeanRadius( 94 name="r_mean_act", count_activated=True, count_unactivated=False 95 ), 96 RadiusStandardDeviation( 97 name="r_std_act", count_activated=True, count_unactivated=False 98 ), 99 ), 100 requested_attributes=( 101 additional_derived_attributes := ( 102 "radius", 103 "equilibrium saturation", 104 ) 105 ), 106 ) 107 ) 108 109 # TODO #1266: copied from G & P 2023 110 self.output_attributes = { 111 attr: tuple([] for _ in range(self.particulator.n_sd)) 112 for attr in additional_derived_attributes 113 }
def
run(self, *, n_steps: int = 2250, steps_per_output_interval: int = 10):
115 def run( 116 self, *, n_steps: int = 2250, steps_per_output_interval: int = 10 117 ): # TODO #1266: essentially copied from G & P 2023 118 output_products = super()._run( 119 nt=n_steps, steps_per_output_interval=steps_per_output_interval 120 ) 121 return {"products": output_products, "attributes": self.output_attributes}