PySDM_examples.Lowe_et_al_2019.simulation
1import numpy as np 2from PySDM_examples.utils import BasicSimulation 3 4import PySDM.products as PySDM_products 5from PySDM import Builder 6from PySDM.dynamics import AmbientThermodynamics, Condensation 7from PySDM.environments import Parcel 8from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii 9from PySDM.initialisation.spectra import Sum 10from PySDM.initialisation.sampling.spectral_sampling import ConstantMultiplicity 11 12 13class Simulation(BasicSimulation): 14 def __init__(self, settings, products=None): 15 n_sd = settings.n_sd_per_mode * len(settings.aerosol.modes) 16 builder = Builder( 17 n_sd=n_sd, 18 backend=settings.backend, 19 environment=Parcel( 20 dt=settings.dt, 21 mass_of_dry_air=settings.mass_of_dry_air, 22 p0=settings.p0, 23 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 24 T0=settings.T0, 25 w=settings.w, 26 ), 27 ) 28 29 attributes = { 30 "dry volume": np.empty(0), 31 "dry volume organic": np.empty(0), 32 "kappa times dry volume": np.empty(0), 33 "multiplicity": np.ndarray(0), 34 } 35 initial_volume = settings.mass_of_dry_air / settings.rho0 36 for mode in settings.aerosol.modes: 37 r_dry, n_in_dv = ConstantMultiplicity( 38 spectrum=mode["spectrum"] 39 ).sample_deterministic(settings.n_sd_per_mode) 40 v_dry = settings.formulae.trivia.volume(radius=r_dry) 41 attributes["multiplicity"] = np.append( 42 attributes["multiplicity"], n_in_dv * initial_volume 43 ) 44 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 45 attributes["dry volume organic"] = np.append( 46 attributes["dry volume organic"], mode["f_org"] * v_dry 47 ) 48 attributes["kappa times dry volume"] = np.append( 49 attributes["kappa times dry volume"], 50 v_dry * mode["kappa"][settings.model], 51 ) 52 for attribute in attributes.values(): 53 assert attribute.shape[0] == n_sd 54 55 np.testing.assert_approx_equal( 56 np.sum(attributes["multiplicity"]) / initial_volume, 57 Sum( 58 tuple( 59 settings.aerosol.modes[i]["spectrum"] 60 for i in range(len(settings.aerosol.modes)) 61 ) 62 ).norm_factor, 63 significant=5, 64 ) 65 r_wet = equilibrate_wet_radii( 66 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 67 environment=builder.particulator.environment, 68 kappa_times_dry_volume=attributes["kappa times dry volume"], 69 f_org=attributes["dry volume organic"] / attributes["dry volume"], 70 ) 71 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 72 73 if settings.model == "Constant": 74 del attributes["dry volume organic"] 75 76 builder.add_dynamic(AmbientThermodynamics()) 77 builder.add_dynamic(Condensation()) 78 79 products = products or ( 80 PySDM_products.ParcelDisplacement(name="z"), 81 PySDM_products.Time(name="t"), 82 PySDM_products.PeakSaturation(name="S_max"), 83 PySDM_products.AmbientRelativeHumidity(name="RH"), 84 PySDM_products.ActivatedParticleConcentration( 85 name="CDNC_cm3", 86 unit="cm^-3", 87 count_activated=True, 88 count_unactivated=False, 89 ), 90 PySDM_products.ParticleSizeSpectrumPerVolume( 91 radius_bins_edges=settings.wet_radius_bins_edges 92 ), 93 PySDM_products.ActivableFraction(name="Activated Fraction"), 94 PySDM_products.WaterMixingRatio(), 95 PySDM_products.AmbientDryAirDensity(name="rhod"), 96 PySDM_products.ActivatedEffectiveRadius( 97 name="reff", count_activated=True, count_unactivated=False 98 ), 99 PySDM_products.ParcelLiquidWaterPath( 100 name="lwp", count_activated=True, count_unactivated=False 101 ), 102 PySDM_products.CloudOpticalDepth(name="tau"), 103 PySDM_products.CloudAlbedo(name="albedo"), 104 ) 105 106 particulator = builder.build(attributes=attributes, products=products) 107 self.settings = settings 108 super().__init__(particulator=particulator) 109 110 def _save_scalars(self, output): 111 for k, v in self.particulator.products.items(): 112 if len(v.shape) > 1 or k in ("lwp", "Activated Fraction", "tau", "albedo"): 113 continue 114 value = v.get() 115 if isinstance(value, np.ndarray) and value.size == 1: 116 value = value[0] 117 output[k].append(value) 118 119 def _save_final_timestep_products(self, output): 120 output["spectrum"] = self.particulator.products[ 121 "particle size spectrum per volume" 122 ].get() 123 124 for name, args_call in { 125 "Activated Fraction": lambda: {"S_max": np.nanmax(output["S_max"])}, 126 "lwp": lambda: {}, 127 "tau": lambda: { 128 "effective_radius": output["reff"][-1], 129 "liquid_water_path": output["lwp"][0], 130 }, 131 "albedo": lambda: {"optical_depth": output["tau"]}, 132 }.items(): 133 output[name] = self.particulator.products[name].get(**args_call()) 134 135 def run(self): 136 output = {k: [] for k in self.particulator.products} 137 for step in self.settings.output_steps: 138 self.particulator.run(step - self.particulator.n_steps) 139 self._save_scalars(output) 140 self._save_final_timestep_products(output) 141 return output
14class Simulation(BasicSimulation): 15 def __init__(self, settings, products=None): 16 n_sd = settings.n_sd_per_mode * len(settings.aerosol.modes) 17 builder = Builder( 18 n_sd=n_sd, 19 backend=settings.backend, 20 environment=Parcel( 21 dt=settings.dt, 22 mass_of_dry_air=settings.mass_of_dry_air, 23 p0=settings.p0, 24 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 25 T0=settings.T0, 26 w=settings.w, 27 ), 28 ) 29 30 attributes = { 31 "dry volume": np.empty(0), 32 "dry volume organic": np.empty(0), 33 "kappa times dry volume": np.empty(0), 34 "multiplicity": np.ndarray(0), 35 } 36 initial_volume = settings.mass_of_dry_air / settings.rho0 37 for mode in settings.aerosol.modes: 38 r_dry, n_in_dv = ConstantMultiplicity( 39 spectrum=mode["spectrum"] 40 ).sample_deterministic(settings.n_sd_per_mode) 41 v_dry = settings.formulae.trivia.volume(radius=r_dry) 42 attributes["multiplicity"] = np.append( 43 attributes["multiplicity"], n_in_dv * initial_volume 44 ) 45 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 46 attributes["dry volume organic"] = np.append( 47 attributes["dry volume organic"], mode["f_org"] * v_dry 48 ) 49 attributes["kappa times dry volume"] = np.append( 50 attributes["kappa times dry volume"], 51 v_dry * mode["kappa"][settings.model], 52 ) 53 for attribute in attributes.values(): 54 assert attribute.shape[0] == n_sd 55 56 np.testing.assert_approx_equal( 57 np.sum(attributes["multiplicity"]) / initial_volume, 58 Sum( 59 tuple( 60 settings.aerosol.modes[i]["spectrum"] 61 for i in range(len(settings.aerosol.modes)) 62 ) 63 ).norm_factor, 64 significant=5, 65 ) 66 r_wet = equilibrate_wet_radii( 67 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 68 environment=builder.particulator.environment, 69 kappa_times_dry_volume=attributes["kappa times dry volume"], 70 f_org=attributes["dry volume organic"] / attributes["dry volume"], 71 ) 72 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 73 74 if settings.model == "Constant": 75 del attributes["dry volume organic"] 76 77 builder.add_dynamic(AmbientThermodynamics()) 78 builder.add_dynamic(Condensation()) 79 80 products = products or ( 81 PySDM_products.ParcelDisplacement(name="z"), 82 PySDM_products.Time(name="t"), 83 PySDM_products.PeakSaturation(name="S_max"), 84 PySDM_products.AmbientRelativeHumidity(name="RH"), 85 PySDM_products.ActivatedParticleConcentration( 86 name="CDNC_cm3", 87 unit="cm^-3", 88 count_activated=True, 89 count_unactivated=False, 90 ), 91 PySDM_products.ParticleSizeSpectrumPerVolume( 92 radius_bins_edges=settings.wet_radius_bins_edges 93 ), 94 PySDM_products.ActivableFraction(name="Activated Fraction"), 95 PySDM_products.WaterMixingRatio(), 96 PySDM_products.AmbientDryAirDensity(name="rhod"), 97 PySDM_products.ActivatedEffectiveRadius( 98 name="reff", count_activated=True, count_unactivated=False 99 ), 100 PySDM_products.ParcelLiquidWaterPath( 101 name="lwp", count_activated=True, count_unactivated=False 102 ), 103 PySDM_products.CloudOpticalDepth(name="tau"), 104 PySDM_products.CloudAlbedo(name="albedo"), 105 ) 106 107 particulator = builder.build(attributes=attributes, products=products) 108 self.settings = settings 109 super().__init__(particulator=particulator) 110 111 def _save_scalars(self, output): 112 for k, v in self.particulator.products.items(): 113 if len(v.shape) > 1 or k in ("lwp", "Activated Fraction", "tau", "albedo"): 114 continue 115 value = v.get() 116 if isinstance(value, np.ndarray) and value.size == 1: 117 value = value[0] 118 output[k].append(value) 119 120 def _save_final_timestep_products(self, output): 121 output["spectrum"] = self.particulator.products[ 122 "particle size spectrum per volume" 123 ].get() 124 125 for name, args_call in { 126 "Activated Fraction": lambda: {"S_max": np.nanmax(output["S_max"])}, 127 "lwp": lambda: {}, 128 "tau": lambda: { 129 "effective_radius": output["reff"][-1], 130 "liquid_water_path": output["lwp"][0], 131 }, 132 "albedo": lambda: {"optical_depth": output["tau"]}, 133 }.items(): 134 output[name] = self.particulator.products[name].get(**args_call()) 135 136 def run(self): 137 output = {k: [] for k in self.particulator.products} 138 for step in self.settings.output_steps: 139 self.particulator.run(step - self.particulator.n_steps) 140 self._save_scalars(output) 141 self._save_final_timestep_products(output) 142 return output
Simulation(settings, products=None)
15 def __init__(self, settings, products=None): 16 n_sd = settings.n_sd_per_mode * len(settings.aerosol.modes) 17 builder = Builder( 18 n_sd=n_sd, 19 backend=settings.backend, 20 environment=Parcel( 21 dt=settings.dt, 22 mass_of_dry_air=settings.mass_of_dry_air, 23 p0=settings.p0, 24 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 25 T0=settings.T0, 26 w=settings.w, 27 ), 28 ) 29 30 attributes = { 31 "dry volume": np.empty(0), 32 "dry volume organic": np.empty(0), 33 "kappa times dry volume": np.empty(0), 34 "multiplicity": np.ndarray(0), 35 } 36 initial_volume = settings.mass_of_dry_air / settings.rho0 37 for mode in settings.aerosol.modes: 38 r_dry, n_in_dv = ConstantMultiplicity( 39 spectrum=mode["spectrum"] 40 ).sample_deterministic(settings.n_sd_per_mode) 41 v_dry = settings.formulae.trivia.volume(radius=r_dry) 42 attributes["multiplicity"] = np.append( 43 attributes["multiplicity"], n_in_dv * initial_volume 44 ) 45 attributes["dry volume"] = np.append(attributes["dry volume"], v_dry) 46 attributes["dry volume organic"] = np.append( 47 attributes["dry volume organic"], mode["f_org"] * v_dry 48 ) 49 attributes["kappa times dry volume"] = np.append( 50 attributes["kappa times dry volume"], 51 v_dry * mode["kappa"][settings.model], 52 ) 53 for attribute in attributes.values(): 54 assert attribute.shape[0] == n_sd 55 56 np.testing.assert_approx_equal( 57 np.sum(attributes["multiplicity"]) / initial_volume, 58 Sum( 59 tuple( 60 settings.aerosol.modes[i]["spectrum"] 61 for i in range(len(settings.aerosol.modes)) 62 ) 63 ).norm_factor, 64 significant=5, 65 ) 66 r_wet = equilibrate_wet_radii( 67 r_dry=settings.formulae.trivia.radius(volume=attributes["dry volume"]), 68 environment=builder.particulator.environment, 69 kappa_times_dry_volume=attributes["kappa times dry volume"], 70 f_org=attributes["dry volume organic"] / attributes["dry volume"], 71 ) 72 attributes["volume"] = settings.formulae.trivia.volume(radius=r_wet) 73 74 if settings.model == "Constant": 75 del attributes["dry volume organic"] 76 77 builder.add_dynamic(AmbientThermodynamics()) 78 builder.add_dynamic(Condensation()) 79 80 products = products or ( 81 PySDM_products.ParcelDisplacement(name="z"), 82 PySDM_products.Time(name="t"), 83 PySDM_products.PeakSaturation(name="S_max"), 84 PySDM_products.AmbientRelativeHumidity(name="RH"), 85 PySDM_products.ActivatedParticleConcentration( 86 name="CDNC_cm3", 87 unit="cm^-3", 88 count_activated=True, 89 count_unactivated=False, 90 ), 91 PySDM_products.ParticleSizeSpectrumPerVolume( 92 radius_bins_edges=settings.wet_radius_bins_edges 93 ), 94 PySDM_products.ActivableFraction(name="Activated Fraction"), 95 PySDM_products.WaterMixingRatio(), 96 PySDM_products.AmbientDryAirDensity(name="rhod"), 97 PySDM_products.ActivatedEffectiveRadius( 98 name="reff", count_activated=True, count_unactivated=False 99 ), 100 PySDM_products.ParcelLiquidWaterPath( 101 name="lwp", count_activated=True, count_unactivated=False 102 ), 103 PySDM_products.CloudOpticalDepth(name="tau"), 104 PySDM_products.CloudAlbedo(name="albedo"), 105 ) 106 107 particulator = builder.build(attributes=attributes, products=products) 108 self.settings = settings 109 super().__init__(particulator=particulator)