Changeset 640 for codes/icosagcm/devel
- Timestamp:
- 12/18/17 12:43:01 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/src/unstructured.pyx
r631 r640 5 5 from ctypes import c_void_p, c_int, c_double, c_bool 6 6 7 from libc.time cimport time as ctime, time_t 8 cimport numpy as np 9 7 10 #------------- direct Cython interface to DYNAMICO routines -------------# 8 11 12 cdef enum: max_nb_stage=5 13 cdef extern : 14 cdef double tauj[max_nb_stage] 15 cdef double cslj[max_nb_stage][max_nb_stage] 16 cdef double cflj[max_nb_stage][max_nb_stage] 17 cdef int nb_stage[1] 18 9 19 cdef extern from "functions.h": 20 cdef void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 21 double *u, double *geopot, double *w, 22 double *theta, double *ps, double *pk, double *hflux, double *qv, 23 double *dmass_col, double *drhodz, double *dtheta_rhodz, 24 double *du_fast, double *du_slow, 25 double *dPhi_fast, double *dPhi_slow, 26 double *dW_fast, double *dW_slow) 10 27 cdef void dynamico_init_params() 11 28 cpdef void dynamico_setup_xios() … … 46 63 47 64 kernels.addvars( 48 c_bool,'hydrostatic','debug_hevi_solver', 'rigid',65 c_bool,'hydrostatic','debug_hevi_solver', 49 66 c_int,'llm','nqdyn','primal_num','max_primal_deg', 50 67 'dual_num','max_dual_deg','edge_num','max_trisk_deg', … … 54 71 55 72 elapsed=0. 73 74 #------------------------ Extension type performing a full ARK time step ---------------------- 75 76 ctypedef double *dbl_ptr 77 cdef dbl_ptr ptr1(double[:] data) : return &data[0] 78 cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0] 79 cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0] 80 cdef dbl_ptr ptr(data): 81 n=data.ndim 82 if n==1 : return ptr1(data) 83 if n==2 : return ptr2(data) 84 if n==3 : return ptr3(data) 85 86 cdef alloc(dbl_ptr *p, allocator, n=1): 87 data=allocator(n) 88 p[0]=ptr(data) 89 return data 90 91 cdef check_ptr(name, dbl_ptr p, np.ndarray data): 92 if p != ptr(data) : print name, 'p <> ptr(data) !!' 93 94 cdef class Caldyn_step: 95 # pointer to allocated arrays 96 cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_w # prognostic 97 cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic 98 cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow # tendencies 99 cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow # tendencies 100 # allocated arrays, must remain referenced or segfault 101 cdef readonly np.ndarray mass, theta_rhodz, u, geopot, w 102 cdef readonly np.ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv 103 cdef readonly np.ndarray drhodz, dtheta_rhodz, du_fast, du_slow 104 cdef readonly np.ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow 105 106 def __init__(self,mesh,time_scheme): 107 # self.mesh=mesh 108 fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass 109 fw, fu, fz = mesh.field_w, mesh.field_u, mesh.field_z 110 # collect coefficients of time scheme 111 cdef double[:] tauj_ = time_scheme.tauj 112 cdef double[:,:] cslj_ = time_scheme.csjl 113 cdef double[:,:] cflj_ = time_scheme.cfjl 114 ns = time_scheme.nstage 115 nb_stage[0]= ns 116 cdef int i,j 117 for i in range(ns): 118 tauj[i]=tauj_[i] 119 for j in range(ns): 120 cslj[i][j]=cslj_[i,j] 121 cflj[i][j]=cflj_[i,j] 122 # allocate arrays, store pointers to avoid overhead when calling dynamico 123 # prognostic/diagnostic 124 self.ps = alloc(&self.p_ps, fps) 125 self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps), 126 self.mass, self.theta_rhodz = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass), 127 self.theta, self.pk = alloc(&self.p_theta, fmass), alloc(&self.p_pk, fmass), 128 self.geopot, self.w = alloc(&self.p_geopot, fw), alloc(&self.p_w, fw), 129 self.hflux, self.u = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu) 130 self.qv = alloc(&self.p_qv,fz) 131 # tendencies 132 self.drhodz, self.dtheta_rhodz = alloc(&self.p_drhodz,fmass,ns), alloc(&self.p_dtheta_rhodz,fmass,ns) 133 self.du_fast, self.du_slow = alloc(&self.p_du_fast,fu,ns), alloc(&self.p_du_slow,fu,ns) 134 self.dPhi_fast, self.dPhi_slow = alloc(&self.p_dPhi_fast,fw,ns), alloc(&self.p_dPhi_slow,fw,ns) 135 self.dW_fast, self.dW_slow = alloc(&self.p_dW_fast,fw,ns), alloc(&self.p_dW_slow,fw,ns) 136 def next(self): 137 # global elapsed 138 # time1=time.time() 139 check_ptr('mass', self.p_mass, self.mass) 140 check_ptr('u', self.p_u, self.u) 141 dynamico_ARK_step(self.p_mass_col, self.p_mass, self.p_theta_rhodz, 142 self.p_u, self.p_geopot, self.p_w, 143 self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv, 144 self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz, 145 self.p_du_fast, self.p_du_slow, 146 self.p_dPhi_fast, self.p_dPhi_slow, 147 self.p_dW_fast, self.p_dW_slow) 148 #time2=time.time() 149 #if time2>time1: elapsed=elapsed+time2-time1 150 def setup_TRSW(self): 151 setvars(('hydrostatic','caldyn_thermo','caldyn_eta'), 152 (True,thermo_boussinesq,eta_lag)) 153 dynamico_init_params() 154 def setup_HPE(self, caldyn_thermo, caldyn_eta, g, BC,thermo): 155 setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 156 'g','ptop','Rd','cpp','preff','Treff'), 157 (True,caldyn_thermo,caldyn_eta, 158 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0)) 159 dynamico_init_params() 160 def setup_NH(self, caldyn_thermo, caldyn_eta, g, BC,thermo): 161 setvars(('hydrostatic','caldyn_thermo','caldyn_eta', 162 'g','ptop','Rd','cpp','preff','Treff', 163 'pbot','rho_bot'), 164 (False,caldyn_thermo,caldyn_eta, 165 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0, 166 BC.pbot.max(), BC.rho_bot.max())) 167 dynamico_init_params() 56 168 57 169 #----------------------------- Base class for dynamics ------------------------
Note: See TracChangeset
for help on using the changeset viewer.