Changeset 642 for codes/icosagcm/devel/Python/src

Timestamp:

12/19/17 15:26:51 (7 years ago)

Author:

dubos

Message:

devel/unstructured : bubble test case with Fortran time stepping

Location:

codes/icosagcm/devel/Python/src

Files:

• functions.h (modified) (1 diff)
• unstructured.pyx (modified) (10 diffs)

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

@@ -7 +7 @@
 void dynamico_init_params(void);
 
-void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 
+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,

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

@@ -18 +18 @@
 
 cdef extern from "functions.h":
-     cdef void dynamico_ARK_step(double *mass_col, double *rhodz, double *theta_rhodz, 
+     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,
@@ -78 +79 @@
 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):
     n=data.ndim
@@ -83 +85 @@
     if n==2 : return ptr2(data)
     if n==3 : return ptr3(data)
-
+    if n==4 : return ptr4(data)
+    if n>4: raise IndexError
+        
 cdef alloc(dbl_ptr *p, allocator, n=1):
     data=allocator(n)
@@ -93 +97 @@
         
 cdef class Caldyn_step:
+    # number of time steps to do at each invocation of advance()
+    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, 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
     # allocated arrays, must remain referenced or segfault
-    cdef readonly np.ndarray mass, theta_rhodz, u, geopot, w
+    cdef readonly np.ndarray mass, theta_rhodz, u, geopot, W
     cdef readonly np.ndarray mass_col, dmass_col, ps, theta, pk, hflux, qv
     cdef readonly np.ndarray drhodz, dtheta_rhodz, du_fast, du_slow
     cdef readonly np.ndarray dPhi_fast, dPhi_slow, dW_fast, dW_slow
 
-    def __init__(self,mesh,time_scheme):
+    def __init__(self,mesh,time_scheme, nstep):
+        self.nstep=nstep
         #        self.mesh=mesh
         fps, ftheta, fmass = mesh.field_ps, mesh.field_theta, mesh.field_mass
@@ -113 +120 @@
         cdef double[:,:] cflj_ = time_scheme.cfjl
         ns = time_scheme.nstage
-        nb_stage[0]= ns
+        nb_stage[0]=ns
+
         cdef int i,j
         for i in range(ns):
@@ -123 +131 @@
         #    prognostic/diagnostic
         self.ps                       = alloc(&self.p_ps, fps) 
-        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps), 
+        self.mass_col, self.dmass_col = alloc(&self.p_mass_col, fps), alloc(&self.p_dmass_col, fps,ns), 
         self.mass, self.theta_rhodz   = alloc(&self.p_mass, fmass), alloc(&self.p_theta_rhodz, fmass),
         self.theta, self.pk           = alloc(&self.p_theta, fmass), alloc(&self.p_pk, fmass), 
-        self.geopot, self.w           = alloc(&self.p_geopot, fw), alloc(&self.p_w, fw),
+        self.geopot, self.W           = alloc(&self.p_geopot, fw), alloc(&self.p_W, fw),
         self.hflux, self.u            = alloc(&self.p_hflux, fu), alloc(&self.p_u, fu) 
         self.qv                       = alloc(&self.p_qv,fz)
@@ -137 +145 @@
         #        global elapsed
         # time1=time.time()
-        check_ptr('mass', self.p_mass, self.mass)
-        check_ptr('u', self.p_u, self.u)
-        dynamico_ARK_step(self.p_mass_col, self.p_mass, self.p_theta_rhodz,
-                          self.p_u, self.p_geopot, self.p_w,
+        dynamico_ARK_step(self.nstep,
+                          self.p_mass_col, self.p_mass, self.p_theta_rhodz,
+                          self.p_u, self.p_geopot, self.p_W,
                           self.p_theta, self.p_ps, self.p_pk, self.p_hflux, self.p_qv,
                           self.p_dmass_col, self.p_drhodz, self.p_dtheta_rhodz,
@@ -148 +155 @@
         #time2=time.time()
         #if time2>time1: elapsed=elapsed+time2-time1
-    def setup_TRSW(self):
-        setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
-                (True,thermo_boussinesq,eta_lag))
-        dynamico_init_params()
-    def setup_HPE(self, caldyn_thermo, caldyn_eta, g, BC,thermo):
-        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
-                 'g','ptop','Rd','cpp','preff','Treff'),
-                (True,caldyn_thermo,caldyn_eta,
-                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0))
-        dynamico_init_params()
-    def setup_NH(self, caldyn_thermo, caldyn_eta, g, BC,thermo):
-        setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
-                 'g','ptop','Rd','cpp','preff','Treff',
-                 'pbot','rho_bot'),
-                (False,caldyn_thermo,caldyn_eta,
-                 g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0,
-                 BC.pbot.max(), BC.rho_bot.max()))
-        dynamico_init_params()
-
+
+def caldyn_step_TRSW(mesh,time_scheme,nstep):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta'),
+            (True,thermo_boussinesq,eta_lag))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_HPE(mesh,time_scheme,nstep, caldyn_thermo,caldyn_eta, thermo,BC,g):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+             'g','ptop','Rd','cpp','preff','Treff'),
+             (True,caldyn_thermo,caldyn_eta,
+              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+def caldyn_step_NH(mesh,time_scheme,nstep, caldyn_thermo, caldyn_eta, thermo,BC,g):
+    setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
+             'g','ptop','Rd','cpp','preff','Treff','pbot','rho_bot'),
+             (False,caldyn_thermo,caldyn_eta,
+              g,BC.ptop,thermo.Rd,thermo.Cpd,thermo.p0,thermo.T0,
+              BC.pbot.max(), BC.rho_bot.max()))
+    dynamico_init_params()
+    return Caldyn_step(mesh,time_scheme, nstep)
+    
 #----------------------------- Base class for dynamics ------------------------
 
@@ -212 +222 @@
 
 class Caldyn_HPE(Caldyn):
-    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
         Caldyn.__init__(self,mesh)
         setvars(('hydrostatic','caldyn_thermo','caldyn_eta',
@@ -231 +241 @@
 
 class Caldyn_NH(Caldyn):
-    def __init__(self,caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g):
+    def __init__(self,caldyn_thermo,caldyn_eta, mesh,thermo,BC,g):
         Caldyn.__init__(self,mesh)
         setvars(('hydrostatic','caldyn_thermo','caldyn_eta',