PySDM_examples.Shipway_and_Hill_2012.mpdata_1d
1from functools import cached_property 2 3import numpy as np 4 5from PyMPDATA import Options, ScalarField, Solver, Stepper, VectorField 6from PyMPDATA.boundary_conditions import Extrapolated 7 8from PySDM.impl import arakawa_c 9 10 11class MPDATA_1D: 12 def __init__( 13 self, 14 nz, 15 dt, 16 advector_of_t, 17 advectee_of_zZ_at_t0, 18 g_factor_of_zZ, 19 mpdata_settings, 20 ): 21 self.__t = 0 22 self.dt = dt 23 self.advector_of_t = advector_of_t 24 25 self._grid = (nz,) 26 self._options = Options( 27 n_iters=mpdata_settings["n_iters"], 28 infinite_gauge=mpdata_settings["iga"], 29 nonoscillatory=mpdata_settings["fct"], 30 third_order_terms=mpdata_settings["tot"], 31 ) 32 bcs = (Extrapolated(),) 33 zZ_scalar = arakawa_c.z_scalar_coord(self._grid) / nz 34 self._g_factor = ScalarField( 35 data=g_factor_of_zZ(zZ_scalar), 36 halo=self._options.n_halo, 37 boundary_conditions=bcs, 38 ) 39 self._advector = VectorField( 40 data=(np.full(nz + 1, advector_of_t(0)),), 41 halo=self._options.n_halo, 42 boundary_conditions=bcs, 43 ) 44 self._advectee = ScalarField( 45 data=advectee_of_zZ_at_t0(zZ_scalar), 46 halo=self._options.n_halo, 47 boundary_conditions=bcs, 48 ) 49 50 @cached_property 51 def solver(self): 52 return Solver( 53 stepper=Stepper( 54 options=self._options, grid=self._grid, non_unit_g_factor=True 55 ), 56 advectee=self._advectee, 57 advector=self._advector, 58 g_factor=self._g_factor, 59 ) 60 61 @property 62 def solver_cached(self): 63 return "solver" in self.__dict__ 64 65 @property 66 def advectee(self): 67 return (self.solver.advectee if self.solver_cached else self._advectee).get() 68 69 @property 70 def advector(self): 71 return ( 72 self.solver.advector if self.solver_cached else self._advector 73 ).get_component(0) 74 75 def update_advector_field(self): 76 self.__t += 0.5 * self.dt 77 self.advector[:] = self.advector_of_t(self.__t) 78 np.testing.assert_array_less(np.abs(self.advector), 1) 79 self.__t += 0.5 * self.dt 80 81 def __call__(self, _): 82 self.solver.advance(1)
class
MPDATA_1D:
12class MPDATA_1D: 13 def __init__( 14 self, 15 nz, 16 dt, 17 advector_of_t, 18 advectee_of_zZ_at_t0, 19 g_factor_of_zZ, 20 mpdata_settings, 21 ): 22 self.__t = 0 23 self.dt = dt 24 self.advector_of_t = advector_of_t 25 26 self._grid = (nz,) 27 self._options = Options( 28 n_iters=mpdata_settings["n_iters"], 29 infinite_gauge=mpdata_settings["iga"], 30 nonoscillatory=mpdata_settings["fct"], 31 third_order_terms=mpdata_settings["tot"], 32 ) 33 bcs = (Extrapolated(),) 34 zZ_scalar = arakawa_c.z_scalar_coord(self._grid) / nz 35 self._g_factor = ScalarField( 36 data=g_factor_of_zZ(zZ_scalar), 37 halo=self._options.n_halo, 38 boundary_conditions=bcs, 39 ) 40 self._advector = VectorField( 41 data=(np.full(nz + 1, advector_of_t(0)),), 42 halo=self._options.n_halo, 43 boundary_conditions=bcs, 44 ) 45 self._advectee = ScalarField( 46 data=advectee_of_zZ_at_t0(zZ_scalar), 47 halo=self._options.n_halo, 48 boundary_conditions=bcs, 49 ) 50 51 @cached_property 52 def solver(self): 53 return Solver( 54 stepper=Stepper( 55 options=self._options, grid=self._grid, non_unit_g_factor=True 56 ), 57 advectee=self._advectee, 58 advector=self._advector, 59 g_factor=self._g_factor, 60 ) 61 62 @property 63 def solver_cached(self): 64 return "solver" in self.__dict__ 65 66 @property 67 def advectee(self): 68 return (self.solver.advectee if self.solver_cached else self._advectee).get() 69 70 @property 71 def advector(self): 72 return ( 73 self.solver.advector if self.solver_cached else self._advector 74 ).get_component(0) 75 76 def update_advector_field(self): 77 self.__t += 0.5 * self.dt 78 self.advector[:] = self.advector_of_t(self.__t) 79 np.testing.assert_array_less(np.abs(self.advector), 1) 80 self.__t += 0.5 * self.dt 81 82 def __call__(self, _): 83 self.solver.advance(1)
MPDATA_1D( nz, dt, advector_of_t, advectee_of_zZ_at_t0, g_factor_of_zZ, mpdata_settings)
13 def __init__( 14 self, 15 nz, 16 dt, 17 advector_of_t, 18 advectee_of_zZ_at_t0, 19 g_factor_of_zZ, 20 mpdata_settings, 21 ): 22 self.__t = 0 23 self.dt = dt 24 self.advector_of_t = advector_of_t 25 26 self._grid = (nz,) 27 self._options = Options( 28 n_iters=mpdata_settings["n_iters"], 29 infinite_gauge=mpdata_settings["iga"], 30 nonoscillatory=mpdata_settings["fct"], 31 third_order_terms=mpdata_settings["tot"], 32 ) 33 bcs = (Extrapolated(),) 34 zZ_scalar = arakawa_c.z_scalar_coord(self._grid) / nz 35 self._g_factor = ScalarField( 36 data=g_factor_of_zZ(zZ_scalar), 37 halo=self._options.n_halo, 38 boundary_conditions=bcs, 39 ) 40 self._advector = VectorField( 41 data=(np.full(nz + 1, advector_of_t(0)),), 42 halo=self._options.n_halo, 43 boundary_conditions=bcs, 44 ) 45 self._advectee = ScalarField( 46 data=advectee_of_zZ_at_t0(zZ_scalar), 47 halo=self._options.n_halo, 48 boundary_conditions=bcs, 49 )