Changeset 689 for codes/icosagcm/devel
- Timestamp:
- 04/09/18 15:24:15 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r687 r689 13 13 14 14 radian=180/math.pi 15 16 def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=unst.np_num) 15 17 16 18 #------------------- Hybrid mass-based coordinate ------------- … … 68 70 69 71 #----------------------- Cartesian mesh ----------------------- 70 71 def squeeze(dims): return np.zeros([n for n in dims if n>1])72 72 73 73 # arrays is a list of arrays … … 177 177 Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee) 178 178 179 def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.ny,self.nx,self.llm))180 def field_mass(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))181 def field_z(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))182 def field_w(self,n=1): return squeeze((n,self.ny,self.nx,self.llm+1))179 def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm)) 180 def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm)) 181 def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm)) 182 def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1)) 183 183 def field_u(self,n=1): 184 184 if n==1: … … 186 186 else: 187 187 return np.zeros((n,self.ny,2*self.nx,self.llm)) 188 def field_ps(self,n=1): return squeeze((n,self.ny,self.nx))188 def field_ps(self,n=1): return zeros((n,self.ny,self.nx)) 189 189 def ucomp(self,u): return u[:,range(0,2*self.nx,2),:] 190 190 def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u … … 279 279 280 280 class Abstract_Mesh(Base_class): 281 def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.primal_num,self.llm))282 def field_mass(self,n=1): return squeeze((n,self.primal_num,self.llm))283 def field_z(self,n=1): return squeeze((n,self.dual_num,self.llm))284 def field_w(self,n=1): return squeeze((n,self.primal_num,self.llm+1))285 def field_u(self,n=1): return squeeze((n,self.edge_num,self.llm))286 def field_ps(self,n=1): return squeeze((n,self.primal_num,))281 def field_theta(self,n=1): return zeros((n,self.nqdyn,self.primal_num,self.llm)) 282 def field_mass(self,n=1): return zeros((n,self.primal_num,self.llm)) 283 def field_z(self,n=1): return zeros((n,self.dual_num,self.llm)) 284 def field_w(self,n=1): return zeros((n,self.primal_num,self.llm+1)) 285 def field_u(self,n=1): return zeros((n,self.edge_num,self.llm)) 286 def field_ps(self,n=1): return zeros((n,self.primal_num,)) 287 287 def ucov2D(self, ulon, ulat): 288 288 return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e)) -
codes/icosagcm/devel/Python/src/functions.h
r642 r689 1 1 /* KERNELS */ 2 2 3 #undef PREC_DOUBLE 4 5 #ifdef PREC_DOUBLE 6 typedef double number; 7 #else 8 typedef float number; 9 #endif 10 3 11 enum {max_nb_stage=5}; 4 extern doubletauj[max_nb_stage];5 extern doublecslj[max_nb_stage][max_nb_stage], cflj[max_nb_stage][max_nb_stage];12 extern number tauj[max_nb_stage]; 13 extern number cslj[max_nb_stage][max_nb_stage], cflj[max_nb_stage][max_nb_stage]; 6 14 7 15 void dynamico_init_params(void); 8 16 9 17 void dynamico_ARK_step(int nstep, 10 double *mass_col, double *rhodz, double *theta_rhodz, 11 double *u, double *geopot, double *w, 12 double *theta, double *ps, double *pk, double *hflux, double *qv, 13 double *dmass_col, double *drhodz, double *dtheta_rhodz, 14 double *du_fast, double *du_slow, 15 double *dPhi_fast, double *dPhi_slow, 16 double *dW_fast, double *dW_slow); 18 number *mass_col, number *rhodz, number *theta_rhodz, 19 number *u, number *geopot, number *w, 20 number *theta, number *ps, number *pk, number *hflux, number *qv, 21 number *dmass_col, number *drhodz, number *dtheta_rhodz, 22 number *du_fast, number *du_slow, 23 number *dPhi_fast, number *dPhi_slow, 24 number *dW_fast, number *dW_slow); 25 26 void dynamico_remap(number *rhodz, number *theta_rhodz, number *u); 17 27 18 28 /* KERNELS -> XIOS */ -
codes/icosagcm/devel/Python/src/unstructured.pyx
r687 r689 5 5 from util import list_stencil 6 6 7 from ctypes import c_void_p, c_int, c_double, c_bool 8 9 #from libc.time cimport time as ctime, time_t 7 from ctypes import c_void_p, c_int, c_double, c_float, c_bool 10 8 from numpy cimport ndarray 11 9 cimport numpy as np 10 import numpy as np 11 12 #-------------- choose precision of kernel computations ------------# 13 14 DEF prec_double=False 15 16 IF prec_double: 17 c_num=c_double 18 ctypedef double num 19 np_num=np.float64 20 ELSE: 21 c_num=c_float 22 ctypedef float num 23 np_num=np.float32 24 25 ctypedef num *num_ptr 12 26 13 27 #------------- direct Cython interface to DYNAMICO routines -------------# 28 14 29 15 30 cdef enum: max_nb_stage=5 16 31 cdef extern : 17 cdef doubletauj[max_nb_stage]18 cdef doublecslj[max_nb_stage][max_nb_stage]19 cdef doublecflj[max_nb_stage][max_nb_stage]32 cdef num tauj[max_nb_stage] 33 cdef num cslj[max_nb_stage][max_nb_stage] 34 cdef num cflj[max_nb_stage][max_nb_stage] 20 35 cdef int nb_stage[1] 21 36 22 37 cdef extern from "functions.h": 23 38 cdef void dynamico_ARK_step(int nstep, 24 double *mass_col, double *rhodz, double*theta_rhodz,25 double *u, double *geopot, double*w,26 double *theta, double *ps, double *pk, double *hflux, double*qv,27 double *dmass_col, double *drhodz, double*dtheta_rhodz,28 double *du_fast, double*du_slow,29 double *dPhi_fast, double*dPhi_slow,30 double *dW_fast, double*dW_slow)31 cdef void dynamico_remap( double *rhodz, double *theta_rhodz, double*u)39 num *mass_col, num *rhodz, num *theta_rhodz, 40 num *u, num *geopot, num *w, 41 num *theta, num *ps, num *pk, num *hflux, num *qv, 42 num *dmass_col, num *drhodz, num *dtheta_rhodz, 43 num *du_fast, num *du_slow, 44 num *dPhi_fast, num *dPhi_slow, 45 num *dW_fast, num *dW_slow) 46 cdef void dynamico_remap(num *rhodz, num *theta_rhodz, num *u) 32 47 cdef void dynamico_init_params() 33 48 cpdef void dynamico_setup_xios() … … 59 74 ['dynamico_init_metric', c_void_p,6], 60 75 ['dynamico_init_hybrid', c_void_p,3], 61 ['dynamico_caldyn_unstructured', c_ double, c_void_p,20],76 ['dynamico_caldyn_unstructured', c_num, c_void_p,20], 62 77 ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p], 63 78 ['dynamico_init_transfer', c_int, c_int,2,c_void_p,3, c_int,2,c_void_p,3], … … 75 90 'dual_num','max_dual_deg','edge_num','max_trisk_deg', 76 91 'caldyn_thermo','caldyn_eta','nb_threads', 77 c_double,'elapsed','g', 'ptop', 'cpp', 'cppv', 92 c_double,'elapsed', 93 c_num, 'g', 'ptop', 'cpp', 'cppv', 78 94 'Rd', 'Rv', 'preff', 'Treff', 'pbot', 'rho_bot', 'Phi_bot') 79 95 … … 82 98 #------------------------ Extension type performing a full ARK time step ---------------------- 83 99 84 ctypedef double *dbl_ptr 85 cdef dbl_ptr ptr1(double[:] data) : return &data[0] 86 cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0] 87 cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0] 88 cdef dbl_ptr ptr4(double[:,:,:,:] data) : return &data[0,0,0,0] 89 cdef dbl_ptr ptr(data): 100 cdef num_ptr ptr1(num[:] data) except *: return &data[0] 101 cdef num_ptr ptr2(num[:,:] data) except *: return &data[0,0] 102 cdef num_ptr ptr3(num[:,:,:] data) except *: return &data[0,0,0] 103 cdef num_ptr ptr4(num[:,:,:,:] data) except *: return &data[0,0,0,0] 104 cdef num_ptr ptr(data) except * : 90 105 n=data.ndim 91 106 if n==1 : return ptr1(data) … … 95 110 if n>4: raise IndexError 96 111 97 cdef alloc( dbl_ptr *p, allocator, n=1):112 cdef alloc(num_ptr *p, allocator, n=1): 98 113 data=allocator(n) 99 114 p[0]=ptr(data) 100 115 return data 101 116 102 cdef check_ptr(name, dbl_ptr p, ndarray data):117 cdef check_ptr(name, num_ptr p, ndarray data): 103 118 if p != ptr(data) : print name, 'p <> ptr(data) !!' 104 119 … … 107 122 cdef int nstep 108 123 # pointer to allocated arrays 109 cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W # prognostic110 cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic111 cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow # tendencies112 cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow # tendencies124 cdef num_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W # prognostic 125 cdef num_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic 126 cdef num_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow # tendencies 127 cdef num_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow # tendencies 113 128 # allocated arrays, must remain referenced or segfault 114 129 cdef readonly ndarray mass, theta_rhodz, u, geopot, W … … 298 313 # i.e. index[start:end] with start=vtxdist[cell], end=vtxdist[cell+1] lists the edges of cell 'cell' 299 314 300 import numpy as np301 302 315 def loc_stencil(degree, stencil): 303 316 loc=0
Note: See TracChangeset
for help on using the changeset viewer.