Changeset 790 for codes/icosagcm

Timestamp:

12/03/18 17:58:11 (5 years ago)

Author:

dubos

Message:

devel/Python : centralize logging and command-line argument parsing + update some test cases accordingly

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/getargs.py (added)
• dynamico/wrap.py (modified) (3 diffs)
• dynamico/xios.py (modified) (5 diffs)
• test/py/Baroclinic_3D_ullrich.py (modified) (14 diffs)
• test/py/test_halo.py (added)
• test/py/write_Cartesian_mesh.py (modified) (2 diffs)
• test/python.sh (modified) (1 diff)

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

@@ -1 +1 @@
 #------------------- wrap C/Fortran routines -------------------#                                                                                           
-
+import getargs
 import numpy as np 
 from ctypes import *
@@ -55 +55 @@
                 name = prefix+name
                 soname = self.prefix_so+name
-                print 'wrap.Sharedlib : importing', soname
+                if verbose : print 'wrap.Sharedlib : importing', soname
                 self.vardict[self.prefix_py+name] = Wrap(self.lib[soname], self.check_args, proto)
     def addvars(self,*types_and_names):
@@ -72 +72 @@
         for name,val in zip(names,vals):
             self.vars[name].value=val
+
+args = getargs.parse()
+verbose = 'wrap' in args.debug

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

@@ +1 @@
+from __future__ import print_function
 import numpy as np
 import cxios
@@ -4 +5 @@
 from dynamico.meshes import radian
 from mpi4py import MPI
+
+import getargs
+args=getargs.parse()
+verbose = 'xios' in args.debug
 
 #-----------------------------------------------------------------------------
@@ -13 +18 @@
         return self
     def __exit__(self, type, value, traceback):
-        print 'xios_finalize()'
+        print('xios_finalize()')
         cxios.finalize()
     def min(self,data): return self.comm.allreduce(data, op=MPI.MIN)
@@ -20 +25 @@
 class Context:
     def __enter__(self):
-        print 'cxios.context_close_definition()'
+        print('cxios.context_close_definition()')
         cxios.context_close_definition()
         return self
     def __exit__(self, type, value, traceback):
-        print 'xios_context_finalize()'
+        print('xios_context_finalize()')
         cxios.context_finalize()
     def init_llm(self, mesh, nqtot):
@@ -62 +67 @@
 def setup_curvilinear(cat,id,nx,ny,own,i_index,j_index,lon,lat):
     mesh = cxios.Handle(cat,id)
-    def log(name,data) : print 'setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max()
+    def log(name,data) : 
+        if verbose : print('setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max() )
     i_index, j_index, lon, lat = [ x[own] for x in i_index, j_index, lon, lat ]
     log('i_index',i_index)

codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py

@@ +1 @@
+from __future__ import print_function
+from dynamico import getargs
+
+getargs.add("--T", type=float, default=3600.,
+                    help="Length of time slice in seconds")
+getargs.add("--perturb", type=float, default=1.,
+                    help="Amplitude of wind perturbation in m/s")
+getargs.add("--N", type=int, default=48,
+                    help="Number of time slices")
+getargs.add("--Davies_N1", type=int, default=3)
+getargs.add("--Davies_N2", type=int, default=3)
+getargs.add("--nx", type=int, default=200)
+getargs.add("--ny", type=int, default=30)
+getargs.add("--betaplane", action='store_true')
+getargs.defaults(dt=360., llm=50)
+
+logging = getargs.Logger('main')
+args = getargs.parse()
+
 from dynamico import unstructured as unst
 from dynamico import dyn
@@ -11 +30 @@
 import numpy as np
 import time
-import argparse
-import logging
 from numpy import pi, log, exp, sin, cos
 
 # Baroclinic instability test based on Ullrich et al. 2015, QJRMS
-
-parser = argparse.ArgumentParser()
-
-parser.add_argument("--mpi_ni", type=int, default=1,
-                    help="number of x processors")
-parser.add_argument("--mpi_nj", type=int, default=1,
-                    help="number of y processors")
-parser.add_argument("--T", type=float, default=3600.,
-                    help="Length of time slice in seconds")
-parser.add_argument("--N", type=int, default=48,
-                    help="Number of time slices")
-parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'),
-                    help="Set verbosity level of logging")
-parser.add_argument("--Davies_N1", type=int, default=3)
-parser.add_argument("--Davies_N2", type=int, default=3)
-parser.add_argument("--llm", type=int, default=50)
-parser.add_argument("--nx", type=int, default=200)
-parser.add_argument("--ny", type=int, default=30)
-parser.add_argument("--dt", type=float, default=360., help='Time step in seconds')
-parser.add_argument("--hydrostatic", action='store_true')
-parser.add_argument("--betaplane", action='store_true')
-
-args = parser.parse_args()
 
 
@@ -52 +46 @@
     Cpd   = 1004.5
     p0    = 1e5
-    up    = 1. # amplitude (m/s)
+    up    = args.perturb # amplitude (m/s)
     xc,yc,lp = 0.,Ly/2.,600000.
 
     omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)    
-    phi0   = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg)
-    f0     = 2*omega*np.sin(phi0) 
-    beta0  = 2*omega*np.cos(phi0)/a if args.betaplane else 0.
+    phi0   = 45.*pi/180.0 # Reference latitude North pi/4 (deg)
+    f0     = 2*omega*sin(phi0) 
+    beta0  = 2*omega*cos(phi0)/a if args.betaplane else 0.
     fb     = f0 - beta0*y0
 
@@ -69 +63 @@
     def ulon(x,y,eta):
         u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 
-#        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))
+        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))
         return  u
 
@@ -103 +97 @@
     eta_ik = eta(alpha_ik)
     eta_ek = eta(alpha_ek)
-    print('min max eta_il', np.min(eta_il),np.max(eta_il))
+#    print('min max eta_il', np.min(eta_il),np.max(eta_il))
 
     Phi_il = Phi_xyeta(y_il, eta_il)
@@ -170 +164 @@
     return div_ik
 
+def KE(mesh, u_ek):
+    llm=u_ek.shape[1]
+    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
+    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
+    div_ik = prec.zeros((primal_num, llm))
+    for i in range(primal_num):
+        div_i = 0.
+        for iedge in range(primal_deg[i]):
+            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
+            div_i = div_i + le_de[e]*(u_ek[e,:]**2)
+        div_ik[i,:] = div_i * (.25/ Ai[i])
+    return div_ik
+
 def curl(mesh, u_ek):
     llm=u_ek.shape[1]
@@ -194 +201 @@
         for i in range(c.size):
             ci=c[i]
-            m[i] = (1.+np.cos((ci-c2)*np.pi/(N2*delta)))/2.0
+            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
             if ci < c2 : m[i]=1.
             if ci > c1 : m[i]=0.
@@ -211 +218 @@
 class myDavies(Davies):
     def mask(self,x,y):
-        logging.debug('davies dy = %f' % dy)
+        logging.debug('davies dy = %f'% dy)
         return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy)
-
-
         
 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
     comm = client.comm
     mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
-
-    #-------------Logging verbosity configuration---------------------------
-    myloglevel = args.log
-    myloglevel = getattr(logging, myloglevel.upper())
-    logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, 
-        format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' )
-
-    logging.debug('%d/%d starting'%(mpi_rank,mpi_size))
+#    getargs.config_log(filename='out.log',filemode='w',
+#                    format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' )
+
+    logging.info('%d/%d starting'%(mpi_rank,mpi_size))
 
     g, Lx, Ly = 9.81, 4e7, 6e6
@@ -237 +238 @@
 
     mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
-    print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk])
     unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
 
@@ -247 +247 @@
     logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt))
     
-    dt=0. # FIXME
-
     caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
 
@@ -264 +262 @@
             caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
         davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e)
-        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )
+#        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )
         logging.debug('mask min/max :')
-        logging.debug(davies.beta_i.min())
-        logging.debug(davies.beta_i.max())
+        logging.debug('%f'% davies.beta_i.min())
+        logging.debug('%f'% davies.beta_i.max())
         def next_flow(m,S,u,Phi,W):
             # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
             caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
             caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
-            kappa = dx**4/43200.
+            kappa = dx**4/432000.
             for i in range(nt):
                 caldyn_step.next()
@@ -282 +280 @@
                 unst.ker.dynamico_scalar_laplacian(s,laps)
                 unst.ker.dynamico_scalar_laplacian(laps,bilaps)
-#                grads = grad(mesh,s)
-#                laps = div(mesh,grads)
-#                gradlaps = grad(mesh,laps) # bilaplacian
-#                bilaps = div(mesh,gradlaps)                
                 caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
@@ -293 +287 @@
     if caldyn_thermo == unst.thermo_theta:
         s=S/m
-        theta = thermo.T0*np.exp(s/thermo.Cpd)
+        theta = thermo.T0*exp(s/thermo.Cpd)
         S=m*theta
         title_format = 'Potential temperature at t=%g s (K)'
@@ -316 +310 @@
             ps = caldyn_step.ps
             curl_vk = curl(mesh, u)
+            KE_ik = KE(mesh,u)
             du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:] 
             div_fast, div_slow = div(mesh,du_fast), div(mesh,du_slow)
+            drhodz, hflux = caldyn_step.drhodz[0,:,:], caldyn_step.hflux[:,:]
+
             # write to disk
             context.send_field_primal('ps', ps)
             context.send_field_primal('temp', gas.T)
+
             context.send_field_primal('p', gas.p)
+#            context.send_field_primal('p', KE_ik)
+#            context.send_field_primal('p', drhodz)
+
             context.send_field_primal('theta', gas.s)
             context.send_field_primal('uz', w)
             context.send_field_primal('z', z)
             context.send_field_primal('div_fast', div_fast)
+
             context.send_field_primal('div_slow', div_slow)
-            print curl_vk.size
+#            context.send_field_primal('div_slow', div(mesh,hflux) )
+
             context.send_field_dual('curl', curl_vk)
 

codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py

@@ +1 @@
+from dynamico import getargs
+getargs.add("--nx", type=int,
+                    default=64, choices=None,
+                    help="number of x points")
+getargs.add("--ny", type=int,
+                    default=64, choices=None,
+                    help="number of y points")
+getargs.add("--Lx", type=float,
+                    default=8., choices=None,
+                    help="Lx")
+getargs.add("--Ly", type=float,
+                    default=8., choices=None,
+                    help="Ly")
+
+args = getargs.parse()
+
 from dynamico.meshes import zeros
 from dynamico import meshes
@@ -4 +20 @@
 import netCDF4 as cdf
 import argparse
-
-parser = argparse.ArgumentParser()
-
-parser.add_argument("-nx", type=int,
-                    default=64, choices=None,
-                    help="number of x points")
-parser.add_argument("-ny", type=int,
-                    default=64, choices=None,
-                    help="number of y points")
-parser.add_argument("-Lx", type=float,
-                    default=8., choices=None,
-                    help="Lx")
-parser.add_argument("-Ly", type=float,
-                    default=8., choices=None,
-                    help="Ly")
-args = parser.parse_args()
 
 nx, ny, Lx, Ly, llm, nqdyn = args.nx, args.ny,args.Lx, args.Ly, 1, 1

codes/icosagcm/devel/Python/test/python.sh

@@ -20 +20 @@
     NB_MPI=$1
     shift
-    mpirun -np $NB_MPI python -u $*
+    rm -f xios_client*.*
+    mpirun -np $NB_MPI python -u $* | tee dynamico.log
 }