Changeset 689 for codes/icosagcm/devel

Timestamp:

04/09/18 15:24:15 (6 years ago)

Author:

dubos

Message:

devel/unstructured : select double or single precision for physical quantities

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/meshes.py (modified) (5 diffs)
• src/functions.h (modified) (1 diff)
• src/unstructured.pyx (modified) (7 diffs)

codes/icosagcm/devel/Python/dynamico/meshes.py

@@ -13 +13 @@
 
 radian=180/math.pi
+
+def zeros(dims): return np.zeros([n for n in dims if n>1], dtype=unst.np_num)
 
 #------------------- Hybrid mass-based coordinate -------------
@@ -68 +70 @@
 
 #----------------------- Cartesian mesh -----------------------
-
-def squeeze(dims): return np.zeros([n for n in dims if n>1])
 
 # arrays is a list of arrays
@@ -177 +177 @@
                   Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee)
 
-    def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.ny,self.nx,self.llm))
-    def field_mass(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))
-    def field_z(self,n=1): return squeeze((n,self.ny,self.nx,self.llm))
-    def field_w(self,n=1): return squeeze((n,self.ny,self.nx,self.llm+1))
+    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
+    def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
+    def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
+    def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))
     def field_u(self,n=1): 
         if n==1:
@@ -186 +186 @@
         else:
             return np.zeros((n,self.ny,2*self.nx,self.llm))
-    def field_ps(self,n=1): return squeeze((n,self.ny,self.nx))
+    def field_ps(self,n=1): return zeros((n,self.ny,self.nx))
     def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]
     def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u
@@ -279 +279 @@
 
 class Abstract_Mesh(Base_class):
-    def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.primal_num,self.llm))
-    def field_mass(self,n=1):  return squeeze((n,self.primal_num,self.llm))
-    def field_z(self,n=1):     return squeeze((n,self.dual_num,self.llm))
-    def field_w(self,n=1):     return squeeze((n,self.primal_num,self.llm+1))
-    def field_u(self,n=1):     return squeeze((n,self.edge_num,self.llm))
-    def field_ps(self,n=1):    return squeeze((n,self.primal_num,))
+    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.primal_num,self.llm))
+    def field_mass(self,n=1):  return zeros((n,self.primal_num,self.llm))
+    def field_z(self,n=1):     return zeros((n,self.dual_num,self.llm))
+    def field_w(self,n=1):     return zeros((n,self.primal_num,self.llm+1))
+    def field_u(self,n=1):     return zeros((n,self.edge_num,self.llm))
+    def field_ps(self,n=1):    return zeros((n,self.primal_num,))
     def ucov2D(self, ulon, ulat): 
         return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e))

codes/icosagcm/devel/Python/src/functions.h

@@ -1 +1 @@
 /* KERNELS */
 
+#undef PREC_DOUBLE
+
+#ifdef PREC_DOUBLE
+typedef double number;
+#else
+typedef float number;
+#endif
+
 enum {max_nb_stage=5};
-extern double tauj[max_nb_stage];
-extern double cslj[max_nb_stage][max_nb_stage], cflj[max_nb_stage][max_nb_stage];
+extern number tauj[max_nb_stage];
+extern number cslj[max_nb_stage][max_nb_stage], cflj[max_nb_stage][max_nb_stage];
 
 void dynamico_init_params(void);
 
 void dynamico_ARK_step(int nstep,
-                       double *mass_col, double *rhodz, double *theta_rhodz, 
-                       double *u, double *geopot, double *w,
-                       double *theta, double *ps, double *pk, double *hflux, double *qv,
-                       double *dmass_col, double *drhodz, double *dtheta_rhodz,
-                       double *du_fast, double *du_slow,
-                       double *dPhi_fast, double *dPhi_slow, 
-                       double *dW_fast, double *dW_slow);
+                       number *mass_col, number *rhodz, number *theta_rhodz, 
+                       number *u, number *geopot, number *w,
+                       number *theta, number *ps, number *pk, number *hflux, number *qv,
+                       number *dmass_col, number *drhodz, number *dtheta_rhodz,
+                       number *du_fast, number *du_slow,
+                       number *dPhi_fast, number *dPhi_slow, 
+                       number *dW_fast, number *dW_slow);
+
+void dynamico_remap(number *rhodz, number *theta_rhodz, number *u);
 
 /* KERNELS -> XIOS */

codes/icosagcm/devel/Python/src/unstructured.pyx

@@ -5 +5 @@
 from util import list_stencil
 
-from ctypes import c_void_p, c_int, c_double, c_bool
-
-#from libc.time cimport time as ctime, time_t
+from ctypes import c_void_p, c_int, c_double, c_float, c_bool
 from numpy cimport ndarray
 cimport numpy as np
+import numpy as np
+
+#--------------   choose precision of kernel computations  ------------#
+
+DEF prec_double=False
+
+IF prec_double:
+   c_num=c_double
+   ctypedef double num
+   np_num=np.float64
+ELSE:
+   c_num=c_float
+   ctypedef float num
+   np_num=np.float32
+
+ctypedef num *num_ptr
 
 #------------- direct Cython interface to DYNAMICO routines -------------#
+
 
 cdef enum: max_nb_stage=5
 cdef extern :
-    cdef double tauj[max_nb_stage]
-    cdef double cslj[max_nb_stage][max_nb_stage]
-    cdef double cflj[max_nb_stage][max_nb_stage]
+    cdef num tauj[max_nb_stage]
+    cdef num cslj[max_nb_stage][max_nb_stage]
+    cdef num cflj[max_nb_stage][max_nb_stage]
     cdef int nb_stage[1]
 
 cdef extern from "functions.h":
      cdef void dynamico_ARK_step(int nstep,
-                                 double *mass_col, double *rhodz, double *theta_rhodz, 
-                                 double *u, double *geopot, double *w,
-                                 double *theta, double *ps, double *pk, double *hflux, double *qv,
-                                 double *dmass_col, double *drhodz, double *dtheta_rhodz,
-                                 double *du_fast, double *du_slow,
-                                 double *dPhi_fast, double *dPhi_slow, 
-                                 double *dW_fast, double *dW_slow)
-     cdef void dynamico_remap(double *rhodz, double *theta_rhodz, double *u)
+                                 num *mass_col, num *rhodz, num *theta_rhodz, 
+                                 num *u, num *geopot, num *w,
+                                 num *theta, num *ps, num *pk, num *hflux, num *qv,
+                                 num *dmass_col, num *drhodz, num *dtheta_rhodz,
+                                 num *du_fast, num *du_slow,
+                                 num *dPhi_fast, num *dPhi_slow, 
+                                 num *dW_fast, num *dW_slow)
+     cdef void dynamico_remap(num *rhodz, num *theta_rhodz, num *u)
      cdef void dynamico_init_params()
      cpdef void dynamico_setup_xios()
@@ -59 +74 @@
                      ['dynamico_init_metric', c_void_p,6],
                      ['dynamico_init_hybrid', c_void_p,3],
-                     ['dynamico_caldyn_unstructured', c_double, c_void_p,20],
+                     ['dynamico_caldyn_unstructured', c_num, c_void_p,20],
                      ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p],
                      ['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',
     'caldyn_thermo','caldyn_eta','nb_threads',
-    c_double,'elapsed','g', 'ptop', 'cpp', 'cppv',
+    c_double,'elapsed',
+    c_num, 'g', 'ptop', 'cpp', 'cppv',
     'Rd', 'Rv', 'preff', 'Treff', 'pbot', 'rho_bot', 'Phi_bot')
 
@@ -82 +98 @@
 #------------------------ Extension type performing a full ARK time step ----------------------
 
-ctypedef double *dbl_ptr
-cdef dbl_ptr ptr1(double[:] data) : return &data[0]
-cdef dbl_ptr ptr2(double[:,:] data) : return &data[0,0]
-cdef dbl_ptr ptr3(double[:,:,:] data) : return &data[0,0,0]
-cdef dbl_ptr ptr4(double[:,:,:,:] data) : return &data[0,0,0,0]
-cdef dbl_ptr ptr(data):
+cdef num_ptr ptr1(num[:] data) except *: return &data[0]
+cdef num_ptr ptr2(num[:,:] data) except *: return &data[0,0]
+cdef num_ptr ptr3(num[:,:,:] data) except *: return &data[0,0,0]
+cdef num_ptr ptr4(num[:,:,:,:] data) except *: return &data[0,0,0,0]
+cdef num_ptr ptr(data) except * :
     n=data.ndim
     if n==1 : return ptr1(data)
@@ -95 +110 @@
     if n>4: raise IndexError
         
-cdef alloc(dbl_ptr *p, allocator, n=1):
+cdef alloc(num_ptr *p, allocator, n=1):
     data=allocator(n)
     p[0]=ptr(data)
     return data
 
-cdef check_ptr(name, dbl_ptr p, ndarray data):
+cdef check_ptr(name, num_ptr p, ndarray data):
     if p != ptr(data) : print name, 'p <> ptr(data) !!'
         
@@ -107 +122 @@
     cdef int nstep
     # pointer to allocated arrays
-    cdef dbl_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W                   # prognostic
-    cdef dbl_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic
-    cdef dbl_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow                # tendencies
-    cdef dbl_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow                # tendencies
+    cdef num_ptr p_mass, p_theta_rhodz, p_u, p_geopot, p_W                   # prognostic
+    cdef num_ptr p_mass_col, p_dmass_col, p_ps, p_theta, p_pk, p_hflux, p_qv # diagnostic
+    cdef num_ptr p_drhodz, p_dtheta_rhodz, p_du_fast, p_du_slow                # tendencies
+    cdef num_ptr p_dPhi_fast, p_dPhi_slow, p_dW_fast, p_dW_slow                # tendencies
     # allocated arrays, must remain referenced or segfault
     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'
 
-import numpy as np
-
 def loc_stencil(degree, stencil):
     loc=0