PyMPDATA_examples.Olesik_et_al_2022.simulation
1import math 2 3import numpy as np 4from PyMPDATA_examples.utils.discretisation import discretised_analytical_solution 5 6from PyMPDATA import Options, ScalarField, Solver, Stepper, VectorField 7from PyMPDATA.boundary_conditions import Constant, Extrapolated 8 9 10class Simulation: 11 @staticmethod 12 def make_condensational_growth_solver( 13 nr, 14 r_min, 15 r_max, 16 GC_max, 17 grid_layout, 18 psi_coord, 19 pdf_of_r, 20 drdt_of_r, 21 opts: Options, 22 ): 23 # psi = psi(p) 24 dp_dr = psi_coord.dx_dr 25 dx_dr = grid_layout.dx_dr 26 27 xh, dx = np.linspace( 28 grid_layout.x(r_min), grid_layout.x(r_max), nr + 1, retstep=True 29 ) 30 rh = grid_layout.r(xh) 31 32 x = np.linspace(xh[0] + dx / 2, xh[-1] - dx / 2, nr) 33 r = grid_layout.r(x) 34 35 def pdf_of_r_over_psi(r): 36 return pdf_of_r(r) / psi_coord.dx_dr(r) 37 38 psi = discretised_analytical_solution( 39 rh, pdf_of_r_over_psi, midpoint_value=True, r=r 40 ) 41 42 dp_dt = drdt_of_r(rh) * dp_dr(rh) 43 G = dp_dr(r) / dx_dr(r) 44 45 # C = dr_dt * dt / dr 46 # GC = dp_dr / dx_dr * dr_dt * dt / dr = 47 # \ \_____ / _..____/ 48 # \_____.._____/ \_ dt/dx 49 # | 50 # dp_dt 51 52 dt = GC_max * dx / np.amax(dp_dt) 53 GCh = dp_dt * dt / dx 54 55 # CFL condition 56 np.testing.assert_array_less(np.abs(GCh), 1) 57 58 g_factor = ScalarField( 59 G.astype(dtype=opts.dtype), 60 halo=opts.n_halo, 61 boundary_conditions=(Extrapolated(),), 62 ) 63 state = ScalarField( 64 psi.astype(dtype=opts.dtype), 65 halo=opts.n_halo, 66 boundary_conditions=(Constant(0),), 67 ) 68 GC_field = VectorField( 69 [GCh.astype(dtype=opts.dtype)], 70 halo=opts.n_halo, 71 boundary_conditions=(Constant(0),), 72 ) 73 stepper = Stepper(options=opts, n_dims=1, non_unit_g_factor=True) 74 return ( 75 Solver( 76 stepper=stepper, g_factor=g_factor, advectee=state, advector=GC_field 77 ), 78 r, 79 rh, 80 dx, 81 dt, 82 g_factor, 83 ) 84 85 @staticmethod 86 def __mgn(quantity, unit): 87 return quantity.to(unit).magnitude 88 89 def __init__(self, settings, grid_layout, psi_coord, opts, GC_max): 90 self.settings = settings 91 self.psi_coord = psi_coord 92 self.grid_layout = grid_layout 93 94 # units of calculation 95 self.__t_unit = self.settings.si.seconds 96 self.__r_unit = self.settings.si.micrometre 97 self.__p_unit = psi_coord.x(self.__r_unit) 98 self.__n_of_r_unit = ( 99 self.settings.si.centimetres**-3 / self.settings.si.micrometre 100 ) 101 102 ( 103 self.solver, 104 self.__r, 105 self.__rh, 106 self.dx, 107 dt, 108 self._g_factor, 109 ) = Simulation.make_condensational_growth_solver( 110 self.settings.nr, 111 self.__mgn(self.settings.r_min, self.__r_unit), 112 self.__mgn(self.settings.r_max, self.__r_unit), 113 GC_max, 114 grid_layout, 115 psi_coord, 116 lambda r: self.__mgn( 117 self.settings.pdf(r * self.__r_unit), self.__n_of_r_unit 118 ), 119 lambda r: self.__mgn( 120 self.settings.drdt(r * self.__r_unit), self.__r_unit / self.__t_unit 121 ), 122 opts, 123 ) 124 125 self.out_steps = tuple(math.ceil(t / dt) for t in settings.out_times) 126 self.dt = dt * self.__t_unit 127 128 def step(self, nt): 129 return self.solver.advance(nt) 130 131 @property 132 def r(self): 133 return self.__r * self.__r_unit 134 135 @property 136 def rh(self): 137 return self.__rh * self.__r_unit 138 139 @property 140 def n_of_r(self): 141 psi = self.solver.advectee.get() 142 n = psi * self.psi_coord.dx_dr(self.__r) 143 return n * self.__n_of_r_unit 144 145 @property 146 def psi(self): 147 psi_unit = self.__n_of_r_unit * self.__r_unit / self.__p_unit 148 return self.solver.advectee.get() * psi_unit 149 150 @property 151 def g_factor(self): 152 return self._g_factor.get() 153 154 @property 155 def dp_dr(self): 156 return self.psi_coord.dx_dr(self.__r)
class
Simulation:
11class Simulation: 12 @staticmethod 13 def make_condensational_growth_solver( 14 nr, 15 r_min, 16 r_max, 17 GC_max, 18 grid_layout, 19 psi_coord, 20 pdf_of_r, 21 drdt_of_r, 22 opts: Options, 23 ): 24 # psi = psi(p) 25 dp_dr = psi_coord.dx_dr 26 dx_dr = grid_layout.dx_dr 27 28 xh, dx = np.linspace( 29 grid_layout.x(r_min), grid_layout.x(r_max), nr + 1, retstep=True 30 ) 31 rh = grid_layout.r(xh) 32 33 x = np.linspace(xh[0] + dx / 2, xh[-1] - dx / 2, nr) 34 r = grid_layout.r(x) 35 36 def pdf_of_r_over_psi(r): 37 return pdf_of_r(r) / psi_coord.dx_dr(r) 38 39 psi = discretised_analytical_solution( 40 rh, pdf_of_r_over_psi, midpoint_value=True, r=r 41 ) 42 43 dp_dt = drdt_of_r(rh) * dp_dr(rh) 44 G = dp_dr(r) / dx_dr(r) 45 46 # C = dr_dt * dt / dr 47 # GC = dp_dr / dx_dr * dr_dt * dt / dr = 48 # \ \_____ / _..____/ 49 # \_____.._____/ \_ dt/dx 50 # | 51 # dp_dt 52 53 dt = GC_max * dx / np.amax(dp_dt) 54 GCh = dp_dt * dt / dx 55 56 # CFL condition 57 np.testing.assert_array_less(np.abs(GCh), 1) 58 59 g_factor = ScalarField( 60 G.astype(dtype=opts.dtype), 61 halo=opts.n_halo, 62 boundary_conditions=(Extrapolated(),), 63 ) 64 state = ScalarField( 65 psi.astype(dtype=opts.dtype), 66 halo=opts.n_halo, 67 boundary_conditions=(Constant(0),), 68 ) 69 GC_field = VectorField( 70 [GCh.astype(dtype=opts.dtype)], 71 halo=opts.n_halo, 72 boundary_conditions=(Constant(0),), 73 ) 74 stepper = Stepper(options=opts, n_dims=1, non_unit_g_factor=True) 75 return ( 76 Solver( 77 stepper=stepper, g_factor=g_factor, advectee=state, advector=GC_field 78 ), 79 r, 80 rh, 81 dx, 82 dt, 83 g_factor, 84 ) 85 86 @staticmethod 87 def __mgn(quantity, unit): 88 return quantity.to(unit).magnitude 89 90 def __init__(self, settings, grid_layout, psi_coord, opts, GC_max): 91 self.settings = settings 92 self.psi_coord = psi_coord 93 self.grid_layout = grid_layout 94 95 # units of calculation 96 self.__t_unit = self.settings.si.seconds 97 self.__r_unit = self.settings.si.micrometre 98 self.__p_unit = psi_coord.x(self.__r_unit) 99 self.__n_of_r_unit = ( 100 self.settings.si.centimetres**-3 / self.settings.si.micrometre 101 ) 102 103 ( 104 self.solver, 105 self.__r, 106 self.__rh, 107 self.dx, 108 dt, 109 self._g_factor, 110 ) = Simulation.make_condensational_growth_solver( 111 self.settings.nr, 112 self.__mgn(self.settings.r_min, self.__r_unit), 113 self.__mgn(self.settings.r_max, self.__r_unit), 114 GC_max, 115 grid_layout, 116 psi_coord, 117 lambda r: self.__mgn( 118 self.settings.pdf(r * self.__r_unit), self.__n_of_r_unit 119 ), 120 lambda r: self.__mgn( 121 self.settings.drdt(r * self.__r_unit), self.__r_unit / self.__t_unit 122 ), 123 opts, 124 ) 125 126 self.out_steps = tuple(math.ceil(t / dt) for t in settings.out_times) 127 self.dt = dt * self.__t_unit 128 129 def step(self, nt): 130 return self.solver.advance(nt) 131 132 @property 133 def r(self): 134 return self.__r * self.__r_unit 135 136 @property 137 def rh(self): 138 return self.__rh * self.__r_unit 139 140 @property 141 def n_of_r(self): 142 psi = self.solver.advectee.get() 143 n = psi * self.psi_coord.dx_dr(self.__r) 144 return n * self.__n_of_r_unit 145 146 @property 147 def psi(self): 148 psi_unit = self.__n_of_r_unit * self.__r_unit / self.__p_unit 149 return self.solver.advectee.get() * psi_unit 150 151 @property 152 def g_factor(self): 153 return self._g_factor.get() 154 155 @property 156 def dp_dr(self): 157 return self.psi_coord.dx_dr(self.__r)
Simulation(settings, grid_layout, psi_coord, opts, GC_max)
90 def __init__(self, settings, grid_layout, psi_coord, opts, GC_max): 91 self.settings = settings 92 self.psi_coord = psi_coord 93 self.grid_layout = grid_layout 94 95 # units of calculation 96 self.__t_unit = self.settings.si.seconds 97 self.__r_unit = self.settings.si.micrometre 98 self.__p_unit = psi_coord.x(self.__r_unit) 99 self.__n_of_r_unit = ( 100 self.settings.si.centimetres**-3 / self.settings.si.micrometre 101 ) 102 103 ( 104 self.solver, 105 self.__r, 106 self.__rh, 107 self.dx, 108 dt, 109 self._g_factor, 110 ) = Simulation.make_condensational_growth_solver( 111 self.settings.nr, 112 self.__mgn(self.settings.r_min, self.__r_unit), 113 self.__mgn(self.settings.r_max, self.__r_unit), 114 GC_max, 115 grid_layout, 116 psi_coord, 117 lambda r: self.__mgn( 118 self.settings.pdf(r * self.__r_unit), self.__n_of_r_unit 119 ), 120 lambda r: self.__mgn( 121 self.settings.drdt(r * self.__r_unit), self.__r_unit / self.__t_unit 122 ), 123 opts, 124 ) 125 126 self.out_steps = tuple(math.ceil(t / dt) for t in settings.out_times) 127 self.dt = dt * self.__t_unit
@staticmethod
def
make_condensational_growth_solver( nr, r_min, r_max, GC_max, grid_layout, psi_coord, pdf_of_r, drdt_of_r, opts: PyMPDATA.options.Options):
12 @staticmethod 13 def make_condensational_growth_solver( 14 nr, 15 r_min, 16 r_max, 17 GC_max, 18 grid_layout, 19 psi_coord, 20 pdf_of_r, 21 drdt_of_r, 22 opts: Options, 23 ): 24 # psi = psi(p) 25 dp_dr = psi_coord.dx_dr 26 dx_dr = grid_layout.dx_dr 27 28 xh, dx = np.linspace( 29 grid_layout.x(r_min), grid_layout.x(r_max), nr + 1, retstep=True 30 ) 31 rh = grid_layout.r(xh) 32 33 x = np.linspace(xh[0] + dx / 2, xh[-1] - dx / 2, nr) 34 r = grid_layout.r(x) 35 36 def pdf_of_r_over_psi(r): 37 return pdf_of_r(r) / psi_coord.dx_dr(r) 38 39 psi = discretised_analytical_solution( 40 rh, pdf_of_r_over_psi, midpoint_value=True, r=r 41 ) 42 43 dp_dt = drdt_of_r(rh) * dp_dr(rh) 44 G = dp_dr(r) / dx_dr(r) 45 46 # C = dr_dt * dt / dr 47 # GC = dp_dr / dx_dr * dr_dt * dt / dr = 48 # \ \_____ / _..____/ 49 # \_____.._____/ \_ dt/dx 50 # | 51 # dp_dt 52 53 dt = GC_max * dx / np.amax(dp_dt) 54 GCh = dp_dt * dt / dx 55 56 # CFL condition 57 np.testing.assert_array_less(np.abs(GCh), 1) 58 59 g_factor = ScalarField( 60 G.astype(dtype=opts.dtype), 61 halo=opts.n_halo, 62 boundary_conditions=(Extrapolated(),), 63 ) 64 state = ScalarField( 65 psi.astype(dtype=opts.dtype), 66 halo=opts.n_halo, 67 boundary_conditions=(Constant(0),), 68 ) 69 GC_field = VectorField( 70 [GCh.astype(dtype=opts.dtype)], 71 halo=opts.n_halo, 72 boundary_conditions=(Constant(0),), 73 ) 74 stepper = Stepper(options=opts, n_dims=1, non_unit_g_factor=True) 75 return ( 76 Solver( 77 stepper=stepper, g_factor=g_factor, advectee=state, advector=GC_field 78 ), 79 r, 80 rh, 81 dx, 82 dt, 83 g_factor, 84 )