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 old_buggy_density_formula=False, 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 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_vapourisation.lv(T) 125 return self.formulae.hydrostatics.drho_dz( 126 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 2 if not old_buggy_density_formula else 1 133 ) 134 135 z_span = (-self.particle_reservoir_depth, self.z_max) 136 z_points = np.linspace(*z_span, 2 * self.nz + 1) 137 rhod_solution = solve_ivp( 138 fun=drhod_dz, 139 t_span=z_span, 140 y0=np.asarray((self.rhod0,)), 141 t_eval=z_points, 142 max_step=dz / 2, 143 ) 144 assert rhod_solution.success 145 self.rhod = interp1d(z_points, rhod_solution.y[0]) 146 147 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 148 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 149 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 150 self.condensation_adaptive = True 151 self.condensation_update_thd = False 152 self.coalescence_adaptive = True 153 154 self.number_of_bins = 100 155 self.r_bins_edges_dry = np.logspace( 156 np.log10(0.001 * si.um), 157 np.log10(1 * si.um), 158 self.number_of_bins + 1, 159 endpoint=True, 160 ) 161 self.r_bins_edges = np.logspace( 162 np.log10(0.001 * si.um), 163 np.log10(10 * si.mm), 164 self.number_of_bins + 1, 165 endpoint=True, 166 ) 167 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 168 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 169 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 170 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 171 self.save_spec_and_attr_times = save_spec_and_attr_times 172 173 @property 174 def n_sd(self): 175 return self.nz * self.n_sd_per_gridbox 176 177 @property 178 def nz(self): 179 assert ( 180 self.particle_reservoir_depth / self.dz 181 == self.particle_reservoir_depth // self.dz 182 ) 183 nz = (self.z_max + self.particle_reservoir_depth) / self.dz 184 assert nz == int(nz) 185 return int(nz) 186 187 @property 188 def nt(self): 189 nt = self.t_max / self.dt 190 assert nt == int(nt) 191 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 old_buggy_density_formula=False, 54 ): 55 self.formulae = formulae or Formulae() 56 self.n_sd_per_gridbox = n_sd_per_gridbox 57 self.p0 = p0 58 self.kappa = kappa 59 self.rho_times_w_1 = rho_times_w_1 60 self.particles_per_volume_STP = particles_per_volume_STP 61 self.dt = dt 62 self.dz = dz 63 self.precip = precip 64 self.enable_condensation = enable_condensation 65 self.z_part = z_part 66 self.z_max = z_max 67 self.t_max = t_max 68 self.collision_kernel = collision_kernel or Geometric(collection_efficiency=1) 69 70 t_1 = 600 * si.s 71 self.rho_times_w = lambda t: ( 72 rho_times_w_1 * np.sin(np.pi * t / t_1) if t < t_1 else 0 73 ) 74 apprx_w1 = rho_times_w_1 / self.formulae.constants.rho_STP 75 self.particle_reservoir_depth = ( 76 (2 * apprx_w1 * t_1 / np.pi) // self.dz + 1 77 ) * self.dz 78 79 self.wet_radius_spectrum_per_mass_of_dry_air = spectra.Lognormal( 80 norm_factor=particles_per_volume_STP / self.formulae.constants.rho_STP, 81 m_mode=0.08 / 2 * si.um, 82 s_geom=1.4, 83 ) 84 85 self._th = interp1d( 86 (0.0 * si.m, 740.0 * si.m, 3260.00 * si.m), 87 (297.9 * si.K, 297.9 * si.K, 312.66 * si.K), 88 fill_value="extrapolate", 89 ) 90 91 self.water_vapour_mixing_ratio = interp1d( 92 (-max(self.particle_reservoir_depth, 1), 0, 740, 3260), 93 (0.015, 0.015, 0.0138, 0.0024), 94 fill_value="extrapolate", 95 ) 96 97 self.thd = ( 98 lambda z_above_reservoir: self.formulae.state_variable_triplet.th_dry( 99 self._th(z_above_reservoir), 100 self.water_vapour_mixing_ratio(z_above_reservoir), 101 ) 102 ) 103 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_vapourisation.lv(T) 126 return self.formulae.hydrostatics.drho_dz( 127 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 2 if not old_buggy_density_formula else 1 134 ) 135 136 z_span = (-self.particle_reservoir_depth, self.z_max) 137 z_points = np.linspace(*z_span, 2 * self.nz + 1) 138 rhod_solution = solve_ivp( 139 fun=drhod_dz, 140 t_span=z_span, 141 y0=np.asarray((self.rhod0,)), 142 t_eval=z_points, 143 max_step=dz / 2, 144 ) 145 assert rhod_solution.success 146 self.rhod = interp1d(z_points, rhod_solution.y[0]) 147 148 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 149 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 150 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 151 self.condensation_adaptive = True 152 self.condensation_update_thd = False 153 self.coalescence_adaptive = True 154 155 self.number_of_bins = 100 156 self.r_bins_edges_dry = np.logspace( 157 np.log10(0.001 * si.um), 158 np.log10(1 * si.um), 159 self.number_of_bins + 1, 160 endpoint=True, 161 ) 162 self.r_bins_edges = np.logspace( 163 np.log10(0.001 * si.um), 164 np.log10(10 * si.mm), 165 self.number_of_bins + 1, 166 endpoint=True, 167 ) 168 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 169 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 170 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 171 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 172 self.save_spec_and_attr_times = save_spec_and_attr_times 173 174 @property 175 def n_sd(self): 176 return self.nz * self.n_sd_per_gridbox 177 178 @property 179 def nz(self): 180 assert ( 181 self.particle_reservoir_depth / self.dz 182 == self.particle_reservoir_depth // self.dz 183 ) 184 nz = (self.z_max + self.particle_reservoir_depth) / self.dz 185 assert nz == int(nz) 186 return int(nz) 187 188 @property 189 def nt(self): 190 nt = self.t_max / self.dt 191 assert nt == int(nt) 192 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, old_buggy_density_formula=False)
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 old_buggy_density_formula=False, 54 ): 55 self.formulae = formulae or Formulae() 56 self.n_sd_per_gridbox = n_sd_per_gridbox 57 self.p0 = p0 58 self.kappa = kappa 59 self.rho_times_w_1 = rho_times_w_1 60 self.particles_per_volume_STP = particles_per_volume_STP 61 self.dt = dt 62 self.dz = dz 63 self.precip = precip 64 self.enable_condensation = enable_condensation 65 self.z_part = z_part 66 self.z_max = z_max 67 self.t_max = t_max 68 self.collision_kernel = collision_kernel or Geometric(collection_efficiency=1) 69 70 t_1 = 600 * si.s 71 self.rho_times_w = lambda t: ( 72 rho_times_w_1 * np.sin(np.pi * t / t_1) if t < t_1 else 0 73 ) 74 apprx_w1 = rho_times_w_1 / self.formulae.constants.rho_STP 75 self.particle_reservoir_depth = ( 76 (2 * apprx_w1 * t_1 / np.pi) // self.dz + 1 77 ) * self.dz 78 79 self.wet_radius_spectrum_per_mass_of_dry_air = spectra.Lognormal( 80 norm_factor=particles_per_volume_STP / self.formulae.constants.rho_STP, 81 m_mode=0.08 / 2 * si.um, 82 s_geom=1.4, 83 ) 84 85 self._th = interp1d( 86 (0.0 * si.m, 740.0 * si.m, 3260.00 * si.m), 87 (297.9 * si.K, 297.9 * si.K, 312.66 * si.K), 88 fill_value="extrapolate", 89 ) 90 91 self.water_vapour_mixing_ratio = interp1d( 92 (-max(self.particle_reservoir_depth, 1), 0, 740, 3260), 93 (0.015, 0.015, 0.0138, 0.0024), 94 fill_value="extrapolate", 95 ) 96 97 self.thd = ( 98 lambda z_above_reservoir: self.formulae.state_variable_triplet.th_dry( 99 self._th(z_above_reservoir), 100 self.water_vapour_mixing_ratio(z_above_reservoir), 101 ) 102 ) 103 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_vapourisation.lv(T) 126 return self.formulae.hydrostatics.drho_dz( 127 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 2 if not old_buggy_density_formula else 1 134 ) 135 136 z_span = (-self.particle_reservoir_depth, self.z_max) 137 z_points = np.linspace(*z_span, 2 * self.nz + 1) 138 rhod_solution = solve_ivp( 139 fun=drhod_dz, 140 t_span=z_span, 141 y0=np.asarray((self.rhod0,)), 142 t_eval=z_points, 143 max_step=dz / 2, 144 ) 145 assert rhod_solution.success 146 self.rhod = interp1d(z_points, rhod_solution.y[0]) 147 148 self.mpdata_settings = {"n_iters": 3, "iga": True, "fct": True, "tot": True} 149 self.condensation_rtol_x = condensation.DEFAULTS.rtol_x 150 self.condensation_rtol_thd = condensation.DEFAULTS.rtol_thd 151 self.condensation_adaptive = True 152 self.condensation_update_thd = False 153 self.coalescence_adaptive = True 154 155 self.number_of_bins = 100 156 self.r_bins_edges_dry = np.logspace( 157 np.log10(0.001 * si.um), 158 np.log10(1 * si.um), 159 self.number_of_bins + 1, 160 endpoint=True, 161 ) 162 self.r_bins_edges = np.logspace( 163 np.log10(0.001 * si.um), 164 np.log10(10 * si.mm), 165 self.number_of_bins + 1, 166 endpoint=True, 167 ) 168 self.cloud_water_radius_range = [1 * si.um, 50 * si.um] 169 self.cloud_water_radius_range_igel = [1 * si.um, 25 * si.um] 170 self.rain_water_radius_range = [50 * si.um, np.inf * si.um] 171 self.rain_water_radius_range_igel = [25 * si.um, np.inf * si.um] 172 self.save_spec_and_attr_times = save_spec_and_attr_times