PySDM_examples.Shipway_and_Hill_2012.settings
1from typing import Iterable, Optional 2 3import numpy as np 4from numdifftools import Derivative 5from scipy.integrate import solve_ivp 6from scipy.interpolate import interp1d 7 8from PySDM import Formulae 9from PySDM.dynamics import condensation 10from PySDM.initialisation import spectra 11from PySDM.physics import si 12from PySDM.dynamics.collisions.collision_kernels import Geometric 13 14 15class Settings: 16 def __dir__(self) -> Iterable[str]: 17 return ( 18 "n_sd_per_gridbox", 19 "p0", 20 "kappa", 21 "rho_times_w_1", 22 "particles_per_volume_STP", 23 "dt", 24 "dz", 25 "precip", 26 "z_max", 27 "t_max", 28 "cloud_water_radius_range", 29 "rain_water_radius_range", 30 "r_bins_edges_dry", 31 "r_bins_edges", 32 ) 33 34 def __init__( 35 self, 36 *, 37 n_sd_per_gridbox: int, 38 p0: float = 1007 * si.hPa, # as used in Olesik et al. 2022 (GMD) 39 kappa: float = 1, 40 rho_times_w_1: float = 2 * si.m / si.s * si.kg / si.m**3, 41 particles_per_volume_STP: int = 50 / si.cm**3, 42 dt: float = 1 * si.s, 43 dz: float = 25 * si.m, 44 z_max: float = 3000 * si.m, 45 z_part: Optional[tuple] = None, 46 t_max: float = 60 * si.minutes, 47 precip: bool = True, 48 enable_condensation: bool = True, 49 formulae: Formulae = None, 50 save_spec_and_attr_times=(), 51 collision_kernel=None 52 ): 53 self.formulae = formulae or Formulae() 54 self.n_sd_per_gridbox = n_sd_per_gridbox 55 self.p0 = p0 56 self.kappa = kappa 57 self.rho_times_w_1 = rho_times_w_1 58 self.particles_per_volume_STP = particles_per_volume_STP 59 self.dt = dt 60 self.dz = dz 61 self.precip = precip 62 self.enable_condensation = enable_condensation 63 self.z_part = z_part 64 self.z_max = z_max 65 self.t_max = t_max 66 self.collision_kernel = collision_kernel or Geometric(collection_efficiency=1) 67 68 t_1 = 600 * si.s 69 self.rho_times_w = lambda t: ( 70 rho_times_w_1 * np.sin(np.pi * t / t_1) if t < t_1 else 0 71 ) 72 apprx_w1 = rho_times_w_1 / self.formulae.constants.rho_STP 73 self.particle_reservoir_depth = ( 74 (2 * apprx_w1 * t_1 / np.pi) // self.dz + 1 75 ) * self.dz 76 77 self.wet_radius_spectrum_per_mass_of_dry_air = spectra.Lognormal( 78 norm_factor=particles_per_volume_STP / self.formulae.constants.rho_STP, 79 m_mode=0.08 / 2 * si.um, 80 s_geom=1.4, 81 ) 82 83 self._th = interp1d( 84 (0.0 * si.m, 740.0 * si.m, 3260.00 * si.m), 85 (297.9 * si.K, 297.9 * si.K, 312.66 * si.K), 86 fill_value="extrapolate", 87 ) 88 89 self.water_vapour_mixing_ratio = interp1d( 90 (-max(self.particle_reservoir_depth, 1), 0, 740, 3260), 91 (0.015, 0.015, 0.0138, 0.0024), 92 fill_value="extrapolate", 93 ) 94 95 self.thd = ( 96 lambda z_above_reservoir: self.formulae.state_variable_triplet.th_dry( 97 self._th(z_above_reservoir), 98 self.water_vapour_mixing_ratio(z_above_reservoir), 99 ) 100 ) 101 102 g = self.formulae.constants.g_std 103 self.rhod0 = self.formulae.state_variable_triplet.rho_d( 104 p=p0, 105 water_vapour_mixing_ratio=self.water_vapour_mixing_ratio(0 * si.m), 106 theta_std=self._th(0 * si.m), 107 ) 108 109 def drhod_dz(z_above_reservoir, rhod): 110 if z_above_reservoir < 0: 111 return 0 112 water_vapour_mixing_ratio = self.water_vapour_mixing_ratio( 113 z_above_reservoir 114 ) 115 d_water_vapour_mixing_ratio__dz = Derivative( 116 self.water_vapour_mixing_ratio 117 )(z_above_reservoir) 118 T = self.formulae.state_variable_triplet.T( 119 rhod[0], self.thd(z_above_reservoir) 120 ) 121 p = self.formulae.state_variable_triplet.p( 122 rhod[0], T, water_vapour_mixing_ratio 123 ) 124 lv = self.formulae.latent_heat.lv(T) 125 return self.formulae.hydrostatics.drho_dz( 126 g, p, T, water_vapour_mixing_ratio, lv 127 ) / ( 128 1 + water_vapour_mixing_ratio 129 ) - rhod * d_water_vapour_mixing_ratio__dz / ( 130 1 + water_vapour_mixing_ratio 131 ) 132 133 z_span = (-self.particle_reservoir_depth, self.z_max) 134 z_points = np.linspace(*z_span, 2 * self.nz + 1) 135 rhod_solution = solve_ivp( 136 fun=drhod_dz, 137 t_span=z_span, 138 y0=np.asarray((self.rhod0,)), 139 t_eval=z_points, 140 max_step=dz / 2, 141 ) 142 assert rhod_solution.success 143 self.rhod = interp1d(z_points, rhod_solution.y[0]) 144 145 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 146 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 147 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 148 self.condensation_adaptive = True 149 self.condensation_update_thd = False 150 self.coalescence_adaptive = True 151 152 self.number_of_bins = 100 153 self.r_bins_edges_dry = np.logspace( 154 np.log10(0.001 * si.um), 155 np.log10(1 * si.um), 156 self.number_of_bins + 1, 157 endpoint=True, 158 ) 159 self.r_bins_edges = np.logspace( 160 np.log10(0.001 * si.um), 161 np.log10(10 * si.mm), 162 self.number_of_bins + 1, 163 endpoint=True, 164 ) 165 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 166 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 167 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 168 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 169 self.save_spec_and_attr_times = save_spec_and_attr_times 170 171 @property 172 def n_sd(self): 173 return self.nz * self.n_sd_per_gridbox 174 175 @property 176 def nz(self): 177 assert ( 178 self.particle_reservoir_depth / self.dz 179 == self.particle_reservoir_depth // self.dz 180 ) 181 nz = (self.z_max + self.particle_reservoir_depth) / self.dz 182 assert nz == int(nz) 183 return int(nz) 184 185 @property 186 def nt(self): 187 nt = self.t_max / self.dt 188 assert nt == int(nt) 189 return int(nt)
class
Settings:
16class Settings: 17 def __dir__(self) -> Iterable[str]: 18 return ( 19 "n_sd_per_gridbox", 20 "p0", 21 "kappa", 22 "rho_times_w_1", 23 "particles_per_volume_STP", 24 "dt", 25 "dz", 26 "precip", 27 "z_max", 28 "t_max", 29 "cloud_water_radius_range", 30 "rain_water_radius_range", 31 "r_bins_edges_dry", 32 "r_bins_edges", 33 ) 34 35 def __init__( 36 self, 37 *, 38 n_sd_per_gridbox: int, 39 p0: float = 1007 * si.hPa, # as used in Olesik et al. 2022 (GMD) 40 kappa: float = 1, 41 rho_times_w_1: float = 2 * si.m / si.s * si.kg / si.m**3, 42 particles_per_volume_STP: int = 50 / si.cm**3, 43 dt: float = 1 * si.s, 44 dz: float = 25 * si.m, 45 z_max: float = 3000 * si.m, 46 z_part: Optional[tuple] = None, 47 t_max: float = 60 * si.minutes, 48 precip: bool = True, 49 enable_condensation: bool = True, 50 formulae: Formulae = None, 51 save_spec_and_attr_times=(), 52 collision_kernel=None 53 ): 54 self.formulae = formulae or Formulae() 55 self.n_sd_per_gridbox = n_sd_per_gridbox 56 self.p0 = p0 57 self.kappa = kappa 58 self.rho_times_w_1 = rho_times_w_1 59 self.particles_per_volume_STP = particles_per_volume_STP 60 self.dt = dt 61 self.dz = dz 62 self.precip = precip 63 self.enable_condensation = enable_condensation 64 self.z_part = z_part 65 self.z_max = z_max 66 self.t_max = t_max 67 self.collision_kernel = collision_kernel or Geometric(collection_efficiency=1) 68 69 t_1 = 600 * si.s 70 self.rho_times_w = lambda t: ( 71 rho_times_w_1 * np.sin(np.pi * t / t_1) if t < t_1 else 0 72 ) 73 apprx_w1 = rho_times_w_1 / self.formulae.constants.rho_STP 74 self.particle_reservoir_depth = ( 75 (2 * apprx_w1 * t_1 / np.pi) // self.dz + 1 76 ) * self.dz 77 78 self.wet_radius_spectrum_per_mass_of_dry_air = spectra.Lognormal( 79 norm_factor=particles_per_volume_STP / self.formulae.constants.rho_STP, 80 m_mode=0.08 / 2 * si.um, 81 s_geom=1.4, 82 ) 83 84 self._th = interp1d( 85 (0.0 * si.m, 740.0 * si.m, 3260.00 * si.m), 86 (297.9 * si.K, 297.9 * si.K, 312.66 * si.K), 87 fill_value="extrapolate", 88 ) 89 90 self.water_vapour_mixing_ratio = interp1d( 91 (-max(self.particle_reservoir_depth, 1), 0, 740, 3260), 92 (0.015, 0.015, 0.0138, 0.0024), 93 fill_value="extrapolate", 94 ) 95 96 self.thd = ( 97 lambda z_above_reservoir: self.formulae.state_variable_triplet.th_dry( 98 self._th(z_above_reservoir), 99 self.water_vapour_mixing_ratio(z_above_reservoir), 100 ) 101 ) 102 103 g = self.formulae.constants.g_std 104 self.rhod0 = self.formulae.state_variable_triplet.rho_d( 105 p=p0, 106 water_vapour_mixing_ratio=self.water_vapour_mixing_ratio(0 * si.m), 107 theta_std=self._th(0 * si.m), 108 ) 109 110 def drhod_dz(z_above_reservoir, rhod): 111 if z_above_reservoir < 0: 112 return 0 113 water_vapour_mixing_ratio = self.water_vapour_mixing_ratio( 114 z_above_reservoir 115 ) 116 d_water_vapour_mixing_ratio__dz = Derivative( 117 self.water_vapour_mixing_ratio 118 )(z_above_reservoir) 119 T = self.formulae.state_variable_triplet.T( 120 rhod[0], self.thd(z_above_reservoir) 121 ) 122 p = self.formulae.state_variable_triplet.p( 123 rhod[0], T, water_vapour_mixing_ratio 124 ) 125 lv = self.formulae.latent_heat.lv(T) 126 return self.formulae.hydrostatics.drho_dz( 127 g, p, T, water_vapour_mixing_ratio, lv 128 ) / ( 129 1 + water_vapour_mixing_ratio 130 ) - rhod * d_water_vapour_mixing_ratio__dz / ( 131 1 + water_vapour_mixing_ratio 132 ) 133 134 z_span = (-self.particle_reservoir_depth, self.z_max) 135 z_points = np.linspace(*z_span, 2 * self.nz + 1) 136 rhod_solution = solve_ivp( 137 fun=drhod_dz, 138 t_span=z_span, 139 y0=np.asarray((self.rhod0,)), 140 t_eval=z_points, 141 max_step=dz / 2, 142 ) 143 assert rhod_solution.success 144 self.rhod = interp1d(z_points, rhod_solution.y[0]) 145 146 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 147 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 148 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 149 self.condensation_adaptive = True 150 self.condensation_update_thd = False 151 self.coalescence_adaptive = True 152 153 self.number_of_bins = 100 154 self.r_bins_edges_dry = np.logspace( 155 np.log10(0.001 * si.um), 156 np.log10(1 * si.um), 157 self.number_of_bins + 1, 158 endpoint=True, 159 ) 160 self.r_bins_edges = np.logspace( 161 np.log10(0.001 * si.um), 162 np.log10(10 * si.mm), 163 self.number_of_bins + 1, 164 endpoint=True, 165 ) 166 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 167 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 168 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 169 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 170 self.save_spec_and_attr_times = save_spec_and_attr_times 171 172 @property 173 def n_sd(self): 174 return self.nz * self.n_sd_per_gridbox 175 176 @property 177 def nz(self): 178 assert ( 179 self.particle_reservoir_depth / self.dz 180 == self.particle_reservoir_depth // self.dz 181 ) 182 nz = (self.z_max + self.particle_reservoir_depth) / self.dz 183 assert nz == int(nz) 184 return int(nz) 185 186 @property 187 def nt(self): 188 nt = self.t_max / self.dt 189 assert nt == int(nt) 190 return int(nt)
Settings( *, n_sd_per_gridbox: int, p0: float = 100700.0, kappa: float = 1, rho_times_w_1: float = 2.0, particles_per_volume_STP: int = 49999999.99999999, dt: float = 1.0, dz: float = 25.0, z_max: float = 3000.0, z_part: Optional[tuple] = None, t_max: float = 3600.0, precip: bool = True, enable_condensation: bool = True, formulae: PySDM.formulae.Formulae = None, save_spec_and_attr_times=(), collision_kernel=None)
35 def __init__( 36 self, 37 *, 38 n_sd_per_gridbox: int, 39 p0: float = 1007 * si.hPa, # as used in Olesik et al. 2022 (GMD) 40 kappa: float = 1, 41 rho_times_w_1: float = 2 * si.m / si.s * si.kg / si.m**3, 42 particles_per_volume_STP: int = 50 / si.cm**3, 43 dt: float = 1 * si.s, 44 dz: float = 25 * si.m, 45 z_max: float = 3000 * si.m, 46 z_part: Optional[tuple] = None, 47 t_max: float = 60 * si.minutes, 48 precip: bool = True, 49 enable_condensation: bool = True, 50 formulae: Formulae = None, 51 save_spec_and_attr_times=(), 52 collision_kernel=None 53 ): 54 self.formulae = formulae or Formulae() 55 self.n_sd_per_gridbox = n_sd_per_gridbox 56 self.p0 = p0 57 self.kappa = kappa 58 self.rho_times_w_1 = rho_times_w_1 59 self.particles_per_volume_STP = particles_per_volume_STP 60 self.dt = dt 61 self.dz = dz 62 self.precip = precip 63 self.enable_condensation = enable_condensation 64 self.z_part = z_part 65 self.z_max = z_max 66 self.t_max = t_max 67 self.collision_kernel = collision_kernel or Geometric(collection_efficiency=1) 68 69 t_1 = 600 * si.s 70 self.rho_times_w = lambda t: ( 71 rho_times_w_1 * np.sin(np.pi * t / t_1) if t < t_1 else 0 72 ) 73 apprx_w1 = rho_times_w_1 / self.formulae.constants.rho_STP 74 self.particle_reservoir_depth = ( 75 (2 * apprx_w1 * t_1 / np.pi) // self.dz + 1 76 ) * self.dz 77 78 self.wet_radius_spectrum_per_mass_of_dry_air = spectra.Lognormal( 79 norm_factor=particles_per_volume_STP / self.formulae.constants.rho_STP, 80 m_mode=0.08 / 2 * si.um, 81 s_geom=1.4, 82 ) 83 84 self._th = interp1d( 85 (0.0 * si.m, 740.0 * si.m, 3260.00 * si.m), 86 (297.9 * si.K, 297.9 * si.K, 312.66 * si.K), 87 fill_value="extrapolate", 88 ) 89 90 self.water_vapour_mixing_ratio = interp1d( 91 (-max(self.particle_reservoir_depth, 1), 0, 740, 3260), 92 (0.015, 0.015, 0.0138, 0.0024), 93 fill_value="extrapolate", 94 ) 95 96 self.thd = ( 97 lambda z_above_reservoir: self.formulae.state_variable_triplet.th_dry( 98 self._th(z_above_reservoir), 99 self.water_vapour_mixing_ratio(z_above_reservoir), 100 ) 101 ) 102 103 g = self.formulae.constants.g_std 104 self.rhod0 = self.formulae.state_variable_triplet.rho_d( 105 p=p0, 106 water_vapour_mixing_ratio=self.water_vapour_mixing_ratio(0 * si.m), 107 theta_std=self._th(0 * si.m), 108 ) 109 110 def drhod_dz(z_above_reservoir, rhod): 111 if z_above_reservoir < 0: 112 return 0 113 water_vapour_mixing_ratio = self.water_vapour_mixing_ratio( 114 z_above_reservoir 115 ) 116 d_water_vapour_mixing_ratio__dz = Derivative( 117 self.water_vapour_mixing_ratio 118 )(z_above_reservoir) 119 T = self.formulae.state_variable_triplet.T( 120 rhod[0], self.thd(z_above_reservoir) 121 ) 122 p = self.formulae.state_variable_triplet.p( 123 rhod[0], T, water_vapour_mixing_ratio 124 ) 125 lv = self.formulae.latent_heat.lv(T) 126 return self.formulae.hydrostatics.drho_dz( 127 g, p, T, water_vapour_mixing_ratio, lv 128 ) / ( 129 1 + water_vapour_mixing_ratio 130 ) - rhod * d_water_vapour_mixing_ratio__dz / ( 131 1 + water_vapour_mixing_ratio 132 ) 133 134 z_span = (-self.particle_reservoir_depth, self.z_max) 135 z_points = np.linspace(*z_span, 2 * self.nz + 1) 136 rhod_solution = solve_ivp( 137 fun=drhod_dz, 138 t_span=z_span, 139 y0=np.asarray((self.rhod0,)), 140 t_eval=z_points, 141 max_step=dz / 2, 142 ) 143 assert rhod_solution.success 144 self.rhod = interp1d(z_points, rhod_solution.y[0]) 145 146 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 147 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 148 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 149 self.condensation_adaptive = True 150 self.condensation_update_thd = False 151 self.coalescence_adaptive = True 152 153 self.number_of_bins = 100 154 self.r_bins_edges_dry = np.logspace( 155 np.log10(0.001 * si.um), 156 np.log10(1 * si.um), 157 self.number_of_bins + 1, 158 endpoint=True, 159 ) 160 self.r_bins_edges = np.logspace( 161 np.log10(0.001 * si.um), 162 np.log10(10 * si.mm), 163 self.number_of_bins + 1, 164 endpoint=True, 165 ) 166 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 167 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 168 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 169 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 170 self.save_spec_and_attr_times = save_spec_and_attr_times