PySDM_examples.Spichtinger_et_al_2023.simulation
1import numpy as np 2 3from PySDM_examples.utils import BasicSimulation 4 5import PySDM.products as PySDM_products 6from PySDM.backends import Numba 7from PySDM import Particulator 8from PySDM.dynamics import ( 9 AmbientThermodynamics, 10 Condensation, 11 Freezing, 12 VapourDepositionOnIce, 13) 14from PySDM.environments import Parcel 15from PySDM.initialisation import discretise_multiplicities 16from PySDM.initialisation.hygroscopic_equilibrium import equilibrate_wet_radii 17 18 19class Simulation(BasicSimulation): 20 def __init__(self, settings, backend=Numba): 21 22 dt = settings.dt 23 24 formulae = settings.formulae 25 26 env = Parcel( 27 mixed_phase=True, 28 dt=dt, 29 mass_of_dry_air=settings.mass_of_dry_air, 30 p0=settings.initial_pressure, 31 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 32 T0=settings.initial_temperature, 33 w=settings.w_updraft, 34 backend=backend( 35 formulae=settings.formulae, 36 **( 37 {"override_jit_flags": {"parallel": False}} 38 if backend is Numba 39 else {} 40 ), 41 ), 42 ) 43 44 self.n_sd = settings.n_sd 45 self.multiplicities = discretise_multiplicities( 46 settings.specific_concentration * env.mass_of_dry_air 47 ) 48 self.r_dry = settings.r_dry 49 v_dry = settings.formulae.trivia.volume(radius=self.r_dry) 50 kappa = settings.kappa 51 52 self.r_wet = equilibrate_wet_radii( 53 r_dry=self.r_dry, 54 environment=env, 55 kappa_times_dry_volume=kappa * v_dry, 56 ) 57 58 attributes = { 59 "multiplicity": self.multiplicities, 60 "dry volume": v_dry, 61 "kappa times dry volume": kappa * v_dry, 62 "signed water mass": formulae.particle_shape_and_density.radius_to_mass( 63 self.r_wet 64 ), 65 } 66 67 products = [ 68 PySDM_products.Time(name="t"), 69 PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"), 70 PySDM_products.ParticleConcentration( 71 name="n_i", unit="1/m**3", radius_range=(-np.inf, 0) 72 ), 73 ] 74 75 self.n_output = settings.n_output 76 self.n_substeps = int(settings.t_duration / dt / self.n_output) 77 super().__init__( 78 Particulator( 79 attributes=attributes, 80 products=products, 81 n_sd=settings.n_sd, 82 environment=env, 83 dynamics=( 84 AmbientThermodynamics(), 85 Condensation(), 86 VapourDepositionOnIce(), 87 Freezing( 88 homogeneous_freezing="time-dependent", immersion_freezing=None 89 ), 90 ), 91 ) 92 ) 93 94 def save(self, output): 95 cell_id = 0 96 output["t"].append(self.particulator.products["t"].get()) 97 output["ni"].append(self.particulator.products["n_i"].get()[cell_id]) 98 output["RHi"].append(self.particulator.products["RH_ice"].get()[cell_id]) 99 100 def run(self): 101 output = { 102 "t": [], 103 "ni": [], 104 "RHi": [], 105 } 106 107 self.save(output) 108 109 RHi_old = self.particulator.products["RH_ice"].get()[0].copy() 110 for _ in range(self.n_output): 111 112 self.particulator.run(self.n_substeps) 113 114 self.save(output) 115 116 RHi = self.particulator.products["RH_ice"].get()[0].copy() 117 dRHi = (RHi_old - RHi) / RHi_old 118 if dRHi > 0.0 and RHi < 130.0: 119 print("break") 120 break 121 RHi_old = RHi 122 123 return output["ni"][-1]
20class Simulation(BasicSimulation): 21 def __init__(self, settings, backend=Numba): 22 23 dt = settings.dt 24 25 formulae = settings.formulae 26 27 env = Parcel( 28 mixed_phase=True, 29 dt=dt, 30 mass_of_dry_air=settings.mass_of_dry_air, 31 p0=settings.initial_pressure, 32 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 33 T0=settings.initial_temperature, 34 w=settings.w_updraft, 35 backend=backend( 36 formulae=settings.formulae, 37 **( 38 {"override_jit_flags": {"parallel": False}} 39 if backend is Numba 40 else {} 41 ), 42 ), 43 ) 44 45 self.n_sd = settings.n_sd 46 self.multiplicities = discretise_multiplicities( 47 settings.specific_concentration * env.mass_of_dry_air 48 ) 49 self.r_dry = settings.r_dry 50 v_dry = settings.formulae.trivia.volume(radius=self.r_dry) 51 kappa = settings.kappa 52 53 self.r_wet = equilibrate_wet_radii( 54 r_dry=self.r_dry, 55 environment=env, 56 kappa_times_dry_volume=kappa * v_dry, 57 ) 58 59 attributes = { 60 "multiplicity": self.multiplicities, 61 "dry volume": v_dry, 62 "kappa times dry volume": kappa * v_dry, 63 "signed water mass": formulae.particle_shape_and_density.radius_to_mass( 64 self.r_wet 65 ), 66 } 67 68 products = [ 69 PySDM_products.Time(name="t"), 70 PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"), 71 PySDM_products.ParticleConcentration( 72 name="n_i", unit="1/m**3", radius_range=(-np.inf, 0) 73 ), 74 ] 75 76 self.n_output = settings.n_output 77 self.n_substeps = int(settings.t_duration / dt / self.n_output) 78 super().__init__( 79 Particulator( 80 attributes=attributes, 81 products=products, 82 n_sd=settings.n_sd, 83 environment=env, 84 dynamics=( 85 AmbientThermodynamics(), 86 Condensation(), 87 VapourDepositionOnIce(), 88 Freezing( 89 homogeneous_freezing="time-dependent", immersion_freezing=None 90 ), 91 ), 92 ) 93 ) 94 95 def save(self, output): 96 cell_id = 0 97 output["t"].append(self.particulator.products["t"].get()) 98 output["ni"].append(self.particulator.products["n_i"].get()[cell_id]) 99 output["RHi"].append(self.particulator.products["RH_ice"].get()[cell_id]) 100 101 def run(self): 102 output = { 103 "t": [], 104 "ni": [], 105 "RHi": [], 106 } 107 108 self.save(output) 109 110 RHi_old = self.particulator.products["RH_ice"].get()[0].copy() 111 for _ in range(self.n_output): 112 113 self.particulator.run(self.n_substeps) 114 115 self.save(output) 116 117 RHi = self.particulator.products["RH_ice"].get()[0].copy() 118 dRHi = (RHi_old - RHi) / RHi_old 119 if dRHi > 0.0 and RHi < 130.0: 120 print("break") 121 break 122 RHi_old = RHi 123 124 return output["ni"][-1]
Simulation(settings, backend=<class 'PySDM.backends.Numba'>)
21 def __init__(self, settings, backend=Numba): 22 23 dt = settings.dt 24 25 formulae = settings.formulae 26 27 env = Parcel( 28 mixed_phase=True, 29 dt=dt, 30 mass_of_dry_air=settings.mass_of_dry_air, 31 p0=settings.initial_pressure, 32 initial_water_vapour_mixing_ratio=settings.initial_water_vapour_mixing_ratio, 33 T0=settings.initial_temperature, 34 w=settings.w_updraft, 35 backend=backend( 36 formulae=settings.formulae, 37 **( 38 {"override_jit_flags": {"parallel": False}} 39 if backend is Numba 40 else {} 41 ), 42 ), 43 ) 44 45 self.n_sd = settings.n_sd 46 self.multiplicities = discretise_multiplicities( 47 settings.specific_concentration * env.mass_of_dry_air 48 ) 49 self.r_dry = settings.r_dry 50 v_dry = settings.formulae.trivia.volume(radius=self.r_dry) 51 kappa = settings.kappa 52 53 self.r_wet = equilibrate_wet_radii( 54 r_dry=self.r_dry, 55 environment=env, 56 kappa_times_dry_volume=kappa * v_dry, 57 ) 58 59 attributes = { 60 "multiplicity": self.multiplicities, 61 "dry volume": v_dry, 62 "kappa times dry volume": kappa * v_dry, 63 "signed water mass": formulae.particle_shape_and_density.radius_to_mass( 64 self.r_wet 65 ), 66 } 67 68 products = [ 69 PySDM_products.Time(name="t"), 70 PySDM_products.AmbientRelativeHumidity(name="RH_ice", unit="%"), 71 PySDM_products.ParticleConcentration( 72 name="n_i", unit="1/m**3", radius_range=(-np.inf, 0) 73 ), 74 ] 75 76 self.n_output = settings.n_output 77 self.n_substeps = int(settings.t_duration / dt / self.n_output) 78 super().__init__( 79 Particulator( 80 attributes=attributes, 81 products=products, 82 n_sd=settings.n_sd, 83 environment=env, 84 dynamics=( 85 AmbientThermodynamics(), 86 Condensation(), 87 VapourDepositionOnIce(), 88 Freezing( 89 homogeneous_freezing="time-dependent", immersion_freezing=None 90 ), 91 ), 92 ) 93 )
def
run(self):
101 def run(self): 102 output = { 103 "t": [], 104 "ni": [], 105 "RHi": [], 106 } 107 108 self.save(output) 109 110 RHi_old = self.particulator.products["RH_ice"].get()[0].copy() 111 for _ in range(self.n_output): 112 113 self.particulator.run(self.n_substeps) 114 115 self.save(output) 116 117 RHi = self.particulator.products["RH_ice"].get()[0].copy() 118 dRHi = (RHi_old - RHi) / RHi_old 119 if dRHi > 0.0 and RHi < 130.0: 120 print("break") 121 break 122 RHi_old = RHi 123 124 return output["ni"][-1]