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