PySDM_examples.Grabowski_and_Pawlowska_2023.simulation
1import numpy as np 2from PySDM_examples.utils import BasicSimulation 3 4from PySDM import Builder 5from PySDM.backends import CPU 6from PySDM.backends.impl_numba.test_helpers import scipy_ode_condensation_solver 7from PySDM.dynamics import AmbientThermodynamics, Condensation 8from PySDM.environments import Parcel 9from PySDM.initialisation import equilibrate_wet_radii 10from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity 11from PySDM.physics import si 12 13 14class Simulation(BasicSimulation): 15 def __init__( 16 self, 17 settings, 18 products=None, 19 scipy_solver=False, 20 sampling_class=ConstantMultiplicity, 21 ): 22 builder = Builder( 23 n_sd=settings.n_sd, 24 backend=CPU( 25 formulae=settings.formulae, override_jit_flags={"parallel": False} 26 ), 27 environment=Parcel( 28 dt=settings.timestep, 29 p0=settings.initial_pressure, 30 initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio, 31 T0=settings.initial_temperature, 32 w=settings.vertical_velocity, 33 mass_of_dry_air=44 * si.kg, 34 ), 35 ) 36 builder.add_dynamic(AmbientThermodynamics()) 37 builder.add_dynamic( 38 Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x) 39 ) 40 for attribute in ( 41 "critical saturation", 42 "equilibrium saturation", 43 "critical volume", 44 ): 45 builder.request_attribute(attribute) 46 47 env = builder.particulator.environment 48 volume = env.mass_of_dry_air / settings.initial_air_density 49 attributes = { 50 k: np.empty(0) 51 for k in ("dry volume", "kappa times dry volume", "multiplicity") 52 } 53 54 assert len(settings.aerosol_modes_by_kappa.keys()) == 1 55 kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0] 56 spectrum = settings.aerosol_modes_by_kappa[kappa] 57 58 r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd) 59 v_dry = settings.formulae.trivia.volume(radius=r_dry) 60 attributes["multiplicity"] = np.append( 61 attributes["multiplicity"], n_per_volume * volume 62 ) 63 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 64 attributes["kappa times dry volume"] = np.append( 65 attributes["kappa times dry volume"], v_dry * kappa 66 ) 67 r_wet = equilibrate_wet_radii( 68 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 69 environment=env, 70 kappa_times_dry_volume=attributes["kappa times dry volume"], 71 ) 72 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 73 74 super().__init__( 75 particulator=builder.build(attributes=attributes, products=products) 76 ) 77 if scipy_solver: 78 scipy_ode_condensation_solver.patch_particulator(self.particulator) 79 80 self.output_attributes = { 81 "volume": tuple([] for _ in range(self.particulator.n_sd)), 82 "dry volume": tuple([] for _ in range(self.particulator.n_sd)), 83 "critical saturation": tuple([] for _ in range(self.particulator.n_sd)), 84 "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)), 85 "critical volume": tuple([] for _ in range(self.particulator.n_sd)), 86 "multiplicity": tuple([] for _ in range(self.particulator.n_sd)), 87 } 88 self.settings = settings 89 90 self.__sanity_checks(attributes, volume) 91 92 def __sanity_checks(self, attributes, volume): 93 for attribute in attributes.values(): 94 assert attribute.shape[0] == self.particulator.n_sd 95 np.testing.assert_approx_equal( 96 sum(attributes["multiplicity"]) / volume, 97 sum( 98 mode.norm_factor 99 for mode in self.settings.aerosol_modes_by_kappa.values() 100 ), 101 significant=4, 102 ) 103 104 def _save(self, output): 105 for key, attr in self.output_attributes.items(): 106 attr_data = self.particulator.attributes[key].to_ndarray() 107 for drop_id in range(self.particulator.n_sd): 108 attr[drop_id].append(attr_data[drop_id]) 109 super()._save(output) 110 111 def run(self): 112 output_products = super()._run( 113 self.settings.nt, self.settings.steps_per_output_interval 114 ) 115 return {"products": output_products, "attributes": self.output_attributes}
15class Simulation(BasicSimulation): 16 def __init__( 17 self, 18 settings, 19 products=None, 20 scipy_solver=False, 21 sampling_class=ConstantMultiplicity, 22 ): 23 builder = Builder( 24 n_sd=settings.n_sd, 25 backend=CPU( 26 formulae=settings.formulae, override_jit_flags={"parallel": False} 27 ), 28 environment=Parcel( 29 dt=settings.timestep, 30 p0=settings.initial_pressure, 31 initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio, 32 T0=settings.initial_temperature, 33 w=settings.vertical_velocity, 34 mass_of_dry_air=44 * si.kg, 35 ), 36 ) 37 builder.add_dynamic(AmbientThermodynamics()) 38 builder.add_dynamic( 39 Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x) 40 ) 41 for attribute in ( 42 "critical saturation", 43 "equilibrium saturation", 44 "critical volume", 45 ): 46 builder.request_attribute(attribute) 47 48 env = builder.particulator.environment 49 volume = env.mass_of_dry_air / settings.initial_air_density 50 attributes = { 51 k: np.empty(0) 52 for k in ("dry volume", "kappa times dry volume", "multiplicity") 53 } 54 55 assert len(settings.aerosol_modes_by_kappa.keys()) == 1 56 kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0] 57 spectrum = settings.aerosol_modes_by_kappa[kappa] 58 59 r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd) 60 v_dry = settings.formulae.trivia.volume(radius=r_dry) 61 attributes["multiplicity"] = np.append( 62 attributes["multiplicity"], n_per_volume * volume 63 ) 64 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 65 attributes["kappa times dry volume"] = np.append( 66 attributes["kappa times dry volume"], v_dry * kappa 67 ) 68 r_wet = equilibrate_wet_radii( 69 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 70 environment=env, 71 kappa_times_dry_volume=attributes["kappa times dry volume"], 72 ) 73 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 74 75 super().__init__( 76 particulator=builder.build(attributes=attributes, products=products) 77 ) 78 if scipy_solver: 79 scipy_ode_condensation_solver.patch_particulator(self.particulator) 80 81 self.output_attributes = { 82 "volume": tuple([] for _ in range(self.particulator.n_sd)), 83 "dry volume": tuple([] for _ in range(self.particulator.n_sd)), 84 "critical saturation": tuple([] for _ in range(self.particulator.n_sd)), 85 "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)), 86 "critical volume": tuple([] for _ in range(self.particulator.n_sd)), 87 "multiplicity": tuple([] for _ in range(self.particulator.n_sd)), 88 } 89 self.settings = settings 90 91 self.__sanity_checks(attributes, volume) 92 93 def __sanity_checks(self, attributes, volume): 94 for attribute in attributes.values(): 95 assert attribute.shape[0] == self.particulator.n_sd 96 np.testing.assert_approx_equal( 97 sum(attributes["multiplicity"]) / volume, 98 sum( 99 mode.norm_factor 100 for mode in self.settings.aerosol_modes_by_kappa.values() 101 ), 102 significant=4, 103 ) 104 105 def _save(self, output): 106 for key, attr in self.output_attributes.items(): 107 attr_data = self.particulator.attributes[key].to_ndarray() 108 for drop_id in range(self.particulator.n_sd): 109 attr[drop_id].append(attr_data[drop_id]) 110 super()._save(output) 111 112 def run(self): 113 output_products = super()._run( 114 self.settings.nt, self.settings.steps_per_output_interval 115 ) 116 return {"products": output_products, "attributes": self.output_attributes}
Simulation( settings, products=None, scipy_solver=False, sampling_class=<class 'PySDM.initialisation.sampling.spectral_sampling.ConstantMultiplicity'>)
16 def __init__( 17 self, 18 settings, 19 products=None, 20 scipy_solver=False, 21 sampling_class=ConstantMultiplicity, 22 ): 23 builder = Builder( 24 n_sd=settings.n_sd, 25 backend=CPU( 26 formulae=settings.formulae, override_jit_flags={"parallel": False} 27 ), 28 environment=Parcel( 29 dt=settings.timestep, 30 p0=settings.initial_pressure, 31 initial_water_vapour_mixing_ratio=settings.initial_vapour_mixing_ratio, 32 T0=settings.initial_temperature, 33 w=settings.vertical_velocity, 34 mass_of_dry_air=44 * si.kg, 35 ), 36 ) 37 builder.add_dynamic(AmbientThermodynamics()) 38 builder.add_dynamic( 39 Condensation(rtol_thd=settings.rtol_thd, rtol_x=settings.rtol_x) 40 ) 41 for attribute in ( 42 "critical saturation", 43 "equilibrium saturation", 44 "critical volume", 45 ): 46 builder.request_attribute(attribute) 47 48 env = builder.particulator.environment 49 volume = env.mass_of_dry_air / settings.initial_air_density 50 attributes = { 51 k: np.empty(0) 52 for k in ("dry volume", "kappa times dry volume", "multiplicity") 53 } 54 55 assert len(settings.aerosol_modes_by_kappa.keys()) == 1 56 kappa = tuple(settings.aerosol_modes_by_kappa.keys())[0] 57 spectrum = settings.aerosol_modes_by_kappa[kappa] 58 59 r_dry, n_per_volume = sampling_class(spectrum).sample(settings.n_sd) 60 v_dry = settings.formulae.trivia.volume(radius=r_dry) 61 attributes["multiplicity"] = np.append( 62 attributes["multiplicity"], n_per_volume * volume 63 ) 64 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 65 attributes["kappa times dry volume"] = np.append( 66 attributes["kappa times dry volume"], v_dry * kappa 67 ) 68 r_wet = equilibrate_wet_radii( 69 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 70 environment=env, 71 kappa_times_dry_volume=attributes["kappa times dry volume"], 72 ) 73 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 74 75 super().__init__( 76 particulator=builder.build(attributes=attributes, products=products) 77 ) 78 if scipy_solver: 79 scipy_ode_condensation_solver.patch_particulator(self.particulator) 80 81 self.output_attributes = { 82 "volume": tuple([] for _ in range(self.particulator.n_sd)), 83 "dry volume": tuple([] for _ in range(self.particulator.n_sd)), 84 "critical saturation": tuple([] for _ in range(self.particulator.n_sd)), 85 "equilibrium saturation": tuple([] for _ in range(self.particulator.n_sd)), 86 "critical volume": tuple([] for _ in range(self.particulator.n_sd)), 87 "multiplicity": tuple([] for _ in range(self.particulator.n_sd)), 88 } 89 self.settings = settings 90 91 self.__sanity_checks(attributes, volume)