PyMPDATA_examples.Shipway_and_Hill_2012.settings
1from typing import Optional 2 3import numpy as np 4from PyMPDATA_examples.Olesik_et_al_2022.settings import ksi_1 as default_ksi_1 5from pystrict import strict 6from scipy.integrate import solve_ivp 7from scipy.interpolate import interp1d 8from scipy.misc import derivative 9 10from . import formulae 11from .arakawa_c import arakawa_c 12from .formulae import const, si 13 14 15@strict 16class Settings: 17 def __init__( 18 self, 19 dt: float, 20 dz: float, 21 rhod_w_const: float, 22 t_max: float = 15 * si.minutes, 23 nr: int = 1, 24 r_min: float = np.nan, 25 r_max: float = np.nan, 26 p0: Optional[float] = None, 27 ksi_1: float = default_ksi_1.to_base_units().magnitude, 28 z_max: float = 3000 * si.metres, 29 apprx_drhod_dz: bool = True, 30 ): 31 self.dt = dt 32 self.dz = dz 33 34 self.nr = nr 35 self.ksi_1 = ksi_1 36 37 self.z_max = z_max 38 self.t_max = t_max 39 40 self.qv = interp1d((0, 740, 3260), (0.015, 0.0138, 0.0024)) 41 self._th = interp1d((0, 740, 3260), (297.9, 297.9, 312.66)) 42 43 # note: not in the paper, 44 # https://github.com/BShipway/KiD/tree/master/src/physconst.f90#L43 45 p0 = p0 or 1000 * si.hPa 46 47 self.rhod0 = formulae.rho_d(p0, self.qv(0), self._th(0)) 48 self.thd = lambda z: formulae.th_dry(self._th(z), self.qv(z)) 49 50 def drhod_dz(z, rhod): 51 T = formulae.temperature(rhod[0], self.thd(z)) 52 p = formulae.pressure(rhod[0], T, self.qv(z)) 53 drhod_dz = formulae.drho_dz(const.g, p, T, self.qv(z), const.lv) 54 if not apprx_drhod_dz: # to resolve issue #335 55 qv = self.qv(z) 56 dqv_dz = derivative(self.qv, z) 57 drhod_dz = drhod_dz / (1 + qv) - rhod * dqv_dz / (1 + qv) 58 return drhod_dz 59 60 z_points = np.arange(0, self.z_max + self.dz / 2, self.dz / 2) 61 rhod_solution = solve_ivp( 62 fun=drhod_dz, 63 t_span=(0, self.z_max), 64 y0=np.asarray((self.rhod0,)), 65 t_eval=z_points, 66 ) 67 assert rhod_solution.success 68 69 self.rhod = interp1d(z_points, rhod_solution.y[0]) 70 71 self.t_1 = 600 * si.s 72 self.rhod_w = lambda t: ( 73 rhod_w_const * np.sin(np.pi * t / self.t_1) if t < self.t_1 else 0 74 ) 75 76 self.r_min = r_min 77 self.r_max = r_max 78 self.bin_boundaries, self.dr = np.linspace( 79 self.r_min, self.r_max, self.nr + 1, retstep=True 80 ) 81 82 self.dr_power = {} 83 for k in (1, 2, 3, 4): 84 self.dr_power[k] = ( 85 self.bin_boundaries[1:] ** k - self.bin_boundaries[:-1] ** k 86 ) 87 self.dr_power[k] = self.dr_power[k].reshape(1, -1).T 88 89 self.z_vec = self.dz * arakawa_c.z_vector_coord((self.nz,)) 90 91 @property 92 def nz(self): 93 nz = self.z_max / self.dz 94 assert nz == int(nz) 95 return int(nz) 96 97 @property 98 def nt(self): 99 nt = self.t_max / self.dt 100 assert nt == int(nt) 101 return int(nt)
@strict
class
Settings:
16@strict 17class Settings: 18 def __init__( 19 self, 20 dt: float, 21 dz: float, 22 rhod_w_const: float, 23 t_max: float = 15 * si.minutes, 24 nr: int = 1, 25 r_min: float = np.nan, 26 r_max: float = np.nan, 27 p0: Optional[float] = None, 28 ksi_1: float = default_ksi_1.to_base_units().magnitude, 29 z_max: float = 3000 * si.metres, 30 apprx_drhod_dz: bool = True, 31 ): 32 self.dt = dt 33 self.dz = dz 34 35 self.nr = nr 36 self.ksi_1 = ksi_1 37 38 self.z_max = z_max 39 self.t_max = t_max 40 41 self.qv = interp1d((0, 740, 3260), (0.015, 0.0138, 0.0024)) 42 self._th = interp1d((0, 740, 3260), (297.9, 297.9, 312.66)) 43 44 # note: not in the paper, 45 # https://github.com/BShipway/KiD/tree/master/src/physconst.f90#L43 46 p0 = p0 or 1000 * si.hPa 47 48 self.rhod0 = formulae.rho_d(p0, self.qv(0), self._th(0)) 49 self.thd = lambda z: formulae.th_dry(self._th(z), self.qv(z)) 50 51 def drhod_dz(z, rhod): 52 T = formulae.temperature(rhod[0], self.thd(z)) 53 p = formulae.pressure(rhod[0], T, self.qv(z)) 54 drhod_dz = formulae.drho_dz(const.g, p, T, self.qv(z), const.lv) 55 if not apprx_drhod_dz: # to resolve issue #335 56 qv = self.qv(z) 57 dqv_dz = derivative(self.qv, z) 58 drhod_dz = drhod_dz / (1 + qv) - rhod * dqv_dz / (1 + qv) 59 return drhod_dz 60 61 z_points = np.arange(0, self.z_max + self.dz / 2, self.dz / 2) 62 rhod_solution = solve_ivp( 63 fun=drhod_dz, 64 t_span=(0, self.z_max), 65 y0=np.asarray((self.rhod0,)), 66 t_eval=z_points, 67 ) 68 assert rhod_solution.success 69 70 self.rhod = interp1d(z_points, rhod_solution.y[0]) 71 72 self.t_1 = 600 * si.s 73 self.rhod_w = lambda t: ( 74 rhod_w_const * np.sin(np.pi * t / self.t_1) if t < self.t_1 else 0 75 ) 76 77 self.r_min = r_min 78 self.r_max = r_max 79 self.bin_boundaries, self.dr = np.linspace( 80 self.r_min, self.r_max, self.nr + 1, retstep=True 81 ) 82 83 self.dr_power = {} 84 for k in (1, 2, 3, 4): 85 self.dr_power[k] = ( 86 self.bin_boundaries[1:] ** k - self.bin_boundaries[:-1] ** k 87 ) 88 self.dr_power[k] = self.dr_power[k].reshape(1, -1).T 89 90 self.z_vec = self.dz * arakawa_c.z_vector_coord((self.nz,)) 91 92 @property 93 def nz(self): 94 nz = self.z_max / self.dz 95 assert nz == int(nz) 96 return int(nz) 97 98 @property 99 def nt(self): 100 nt = self.t_max / self.dt 101 assert nt == int(nt) 102 return int(nt)
Settings( dt: float, dz: float, rhod_w_const: float, t_max: float = 900.0, nr: int = 1, r_min: float = nan, r_max: float = nan, p0: Optional[float] = None, ksi_1: float = 1e-10, z_max: float = 3000.0, apprx_drhod_dz: bool = True)
18 def __init__( 19 self, 20 dt: float, 21 dz: float, 22 rhod_w_const: float, 23 t_max: float = 15 * si.minutes, 24 nr: int = 1, 25 r_min: float = np.nan, 26 r_max: float = np.nan, 27 p0: Optional[float] = None, 28 ksi_1: float = default_ksi_1.to_base_units().magnitude, 29 z_max: float = 3000 * si.metres, 30 apprx_drhod_dz: bool = True, 31 ): 32 self.dt = dt 33 self.dz = dz 34 35 self.nr = nr 36 self.ksi_1 = ksi_1 37 38 self.z_max = z_max 39 self.t_max = t_max 40 41 self.qv = interp1d((0, 740, 3260), (0.015, 0.0138, 0.0024)) 42 self._th = interp1d((0, 740, 3260), (297.9, 297.9, 312.66)) 43 44 # note: not in the paper, 45 # https://github.com/BShipway/KiD/tree/master/src/physconst.f90#L43 46 p0 = p0 or 1000 * si.hPa 47 48 self.rhod0 = formulae.rho_d(p0, self.qv(0), self._th(0)) 49 self.thd = lambda z: formulae.th_dry(self._th(z), self.qv(z)) 50 51 def drhod_dz(z, rhod): 52 T = formulae.temperature(rhod[0], self.thd(z)) 53 p = formulae.pressure(rhod[0], T, self.qv(z)) 54 drhod_dz = formulae.drho_dz(const.g, p, T, self.qv(z), const.lv) 55 if not apprx_drhod_dz: # to resolve issue #335 56 qv = self.qv(z) 57 dqv_dz = derivative(self.qv, z) 58 drhod_dz = drhod_dz / (1 + qv) - rhod * dqv_dz / (1 + qv) 59 return drhod_dz 60 61 z_points = np.arange(0, self.z_max + self.dz / 2, self.dz / 2) 62 rhod_solution = solve_ivp( 63 fun=drhod_dz, 64 t_span=(0, self.z_max), 65 y0=np.asarray((self.rhod0,)), 66 t_eval=z_points, 67 ) 68 assert rhod_solution.success 69 70 self.rhod = interp1d(z_points, rhod_solution.y[0]) 71 72 self.t_1 = 600 * si.s 73 self.rhod_w = lambda t: ( 74 rhod_w_const * np.sin(np.pi * t / self.t_1) if t < self.t_1 else 0 75 ) 76 77 self.r_min = r_min 78 self.r_max = r_max 79 self.bin_boundaries, self.dr = np.linspace( 80 self.r_min, self.r_max, self.nr + 1, retstep=True 81 ) 82 83 self.dr_power = {} 84 for k in (1, 2, 3, 4): 85 self.dr_power[k] = ( 86 self.bin_boundaries[1:] ** k - self.bin_boundaries[:-1] ** k 87 ) 88 self.dr_power[k] = self.dr_power[k].reshape(1, -1).T 89 90 self.z_vec = self.dz * arakawa_c.z_vector_coord((self.nz,))