Changeset 790
- Timestamp:
- 12/03/18 17:58:11 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 2 added
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/wrap.py
r747 r790 1 1 #------------------- wrap C/Fortran routines -------------------# 2 2 import getargs 3 3 import numpy as np 4 4 from ctypes import * … … 55 55 name = prefix+name 56 56 soname = self.prefix_so+name 57 print 'wrap.Sharedlib : importing', soname57 if verbose : print 'wrap.Sharedlib : importing', soname 58 58 self.vardict[self.prefix_py+name] = Wrap(self.lib[soname], self.check_args, proto) 59 59 def addvars(self,*types_and_names): … … 72 72 for name,val in zip(names,vals): 73 73 self.vars[name].value=val 74 75 args = getargs.parse() 76 verbose = 'wrap' in args.debug -
codes/icosagcm/devel/Python/dynamico/xios.py
r787 r790 1 from __future__ import print_function 1 2 import numpy as np 2 3 import cxios … … 4 5 from dynamico.meshes import radian 5 6 from mpi4py import MPI 7 8 import getargs 9 args=getargs.parse() 10 verbose = 'xios' in args.debug 6 11 7 12 #----------------------------------------------------------------------------- … … 13 18 return self 14 19 def __exit__(self, type, value, traceback): 15 print 'xios_finalize()'20 print('xios_finalize()') 16 21 cxios.finalize() 17 22 def min(self,data): return self.comm.allreduce(data, op=MPI.MIN) … … 20 25 class Context: 21 26 def __enter__(self): 22 print 'cxios.context_close_definition()'27 print('cxios.context_close_definition()') 23 28 cxios.context_close_definition() 24 29 return self 25 30 def __exit__(self, type, value, traceback): 26 print 'xios_context_finalize()'31 print('xios_context_finalize()') 27 32 cxios.context_finalize() 28 33 def init_llm(self, mesh, nqtot): … … 62 67 def setup_curvilinear(cat,id,nx,ny,own,i_index,j_index,lon,lat): 63 68 mesh = cxios.Handle(cat,id) 64 def log(name,data) : print 'setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max() 69 def log(name,data) : 70 if verbose : print('setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max() ) 65 71 i_index, j_index, lon, lat = [ x[own] for x in i_index, j_index, lon, lat ] 66 72 log('i_index',i_index) -
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r787 r790 1 from __future__ import print_function 2 from dynamico import getargs 3 4 getargs.add("--T", type=float, default=3600., 5 help="Length of time slice in seconds") 6 getargs.add("--perturb", type=float, default=1., 7 help="Amplitude of wind perturbation in m/s") 8 getargs.add("--N", type=int, default=48, 9 help="Number of time slices") 10 getargs.add("--Davies_N1", type=int, default=3) 11 getargs.add("--Davies_N2", type=int, default=3) 12 getargs.add("--nx", type=int, default=200) 13 getargs.add("--ny", type=int, default=30) 14 getargs.add("--betaplane", action='store_true') 15 getargs.defaults(dt=360., llm=50) 16 17 logging = getargs.Logger('main') 18 args = getargs.parse() 19 1 20 from dynamico import unstructured as unst 2 21 from dynamico import dyn … … 11 30 import numpy as np 12 31 import time 13 import argparse14 import logging15 32 from numpy import pi, log, exp, sin, cos 16 33 17 34 # Baroclinic instability test based on Ullrich et al. 2015, QJRMS 18 19 parser = argparse.ArgumentParser()20 21 parser.add_argument("--mpi_ni", type=int, default=1,22 help="number of x processors")23 parser.add_argument("--mpi_nj", type=int, default=1,24 help="number of y processors")25 parser.add_argument("--T", type=float, default=3600.,26 help="Length of time slice in seconds")27 parser.add_argument("--N", type=int, default=48,28 help="Number of time slices")29 parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'),30 help="Set verbosity level of logging")31 parser.add_argument("--Davies_N1", type=int, default=3)32 parser.add_argument("--Davies_N2", type=int, default=3)33 parser.add_argument("--llm", type=int, default=50)34 parser.add_argument("--nx", type=int, default=200)35 parser.add_argument("--ny", type=int, default=30)36 parser.add_argument("--dt", type=float, default=360., help='Time step in seconds')37 parser.add_argument("--hydrostatic", action='store_true')38 parser.add_argument("--betaplane", action='store_true')39 40 args = parser.parse_args()41 35 42 36 … … 52 46 Cpd = 1004.5 53 47 p0 = 1e5 54 up = 1.# amplitude (m/s)48 up = args.perturb # amplitude (m/s) 55 49 xc,yc,lp = 0.,Ly/2.,600000. 56 50 57 51 omega = 7.292e-5 # Angular velocity of the Earth (s^-1) 58 phi0 = 45.* np.pi/180.0 # Reference latitude North pi/4 (deg)59 f0 = 2*omega* np.sin(phi0)60 beta0 = 2*omega* np.cos(phi0)/a if args.betaplane else 0.52 phi0 = 45.*pi/180.0 # Reference latitude North pi/4 (deg) 53 f0 = 2*omega*sin(phi0) 54 beta0 = 2*omega*cos(phi0)/a if args.betaplane else 0. 61 55 fb = f0 - beta0*y0 62 56 … … 69 63 def ulon(x,y,eta): 70 64 u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 71 #u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))65 u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) 72 66 return u 73 67 … … 103 97 eta_ik = eta(alpha_ik) 104 98 eta_ek = eta(alpha_ek) 105 print('min max eta_il', np.min(eta_il),np.max(eta_il))99 # print('min max eta_il', np.min(eta_il),np.max(eta_il)) 106 100 107 101 Phi_il = Phi_xyeta(y_il, eta_il) … … 170 164 return div_ik 171 165 166 def KE(mesh, u_ek): 167 llm=u_ek.shape[1] 168 primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge 169 primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai 170 div_ik = prec.zeros((primal_num, llm)) 171 for i in range(primal_num): 172 div_i = 0. 173 for iedge in range(primal_deg[i]): 174 e = primal_edge[i,iedge]-1 # Fortran => Python indexing 175 div_i = div_i + le_de[e]*(u_ek[e,:]**2) 176 div_ik[i,:] = div_i * (.25/ Ai[i]) 177 return div_ik 178 172 179 def curl(mesh, u_ek): 173 180 llm=u_ek.shape[1] … … 194 201 for i in range(c.size): 195 202 ci=c[i] 196 m[i] = (1.+ np.cos((ci-c2)*np.pi/(N2*delta)))/2.0203 m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0 197 204 if ci < c2 : m[i]=1. 198 205 if ci > c1 : m[i]=0. … … 211 218 class myDavies(Davies): 212 219 def mask(self,x,y): 213 logging.debug('davies dy = %f' 220 logging.debug('davies dy = %f'% dy) 214 221 return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy) 215 216 217 222 218 223 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator 219 224 comm = client.comm 220 225 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 221 222 #-------------Logging verbosity configuration--------------------------- 223 myloglevel = args.log 224 myloglevel = getattr(logging, myloglevel.upper()) 225 logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, 226 format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' ) 227 228 logging.debug('%d/%d starting'%(mpi_rank,mpi_size)) 226 # getargs.config_log(filename='out.log',filemode='w', 227 # format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' ) 228 229 logging.info('%d/%d starting'%(mpi_rank,mpi_size)) 229 230 230 231 g, Lx, Ly = 9.81, 4e7, 6e6 … … 237 238 238 239 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 239 print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk])240 240 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 241 241 … … 247 247 logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) 248 248 249 dt=0. # FIXME250 251 249 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 252 250 … … 264 262 caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 265 263 davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 266 print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )264 # print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) 267 265 logging.debug('mask min/max :') 268 logging.debug( davies.beta_i.min())269 logging.debug( davies.beta_i.max())266 logging.debug('%f'% davies.beta_i.min()) 267 logging.debug('%f'% davies.beta_i.max()) 270 268 def next_flow(m,S,u,Phi,W): 271 269 # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 272 270 caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 273 271 caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W 274 kappa = dx**4/43200 .272 kappa = dx**4/432000. 275 273 for i in range(nt): 276 274 caldyn_step.next() … … 282 280 unst.ker.dynamico_scalar_laplacian(s,laps) 283 281 unst.ker.dynamico_scalar_laplacian(laps,bilaps) 284 # grads = grad(mesh,s)285 # laps = div(mesh,grads)286 # gradlaps = grad(mesh,laps) # bilaplacian287 # bilaps = div(mesh,gradlaps)288 282 caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step 289 283 return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), … … 293 287 if caldyn_thermo == unst.thermo_theta: 294 288 s=S/m 295 theta = thermo.T0* np.exp(s/thermo.Cpd)289 theta = thermo.T0*exp(s/thermo.Cpd) 296 290 S=m*theta 297 291 title_format = 'Potential temperature at t=%g s (K)' … … 316 310 ps = caldyn_step.ps 317 311 curl_vk = curl(mesh, u) 312 KE_ik = KE(mesh,u) 318 313 du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:] 319 314 div_fast, div_slow = div(mesh,du_fast), div(mesh,du_slow) 315 drhodz, hflux = caldyn_step.drhodz[0,:,:], caldyn_step.hflux[:,:] 316 320 317 # write to disk 321 318 context.send_field_primal('ps', ps) 322 319 context.send_field_primal('temp', gas.T) 320 323 321 context.send_field_primal('p', gas.p) 322 # context.send_field_primal('p', KE_ik) 323 # context.send_field_primal('p', drhodz) 324 324 325 context.send_field_primal('theta', gas.s) 325 326 context.send_field_primal('uz', w) 326 327 context.send_field_primal('z', z) 327 328 context.send_field_primal('div_fast', div_fast) 329 328 330 context.send_field_primal('div_slow', div_slow) 329 print curl_vk.size 331 # context.send_field_primal('div_slow', div(mesh,hflux) ) 332 330 333 context.send_field_dual('curl', curl_vk) 331 334 -
codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py
r773 r790 1 from dynamico import getargs 2 getargs.add("--nx", type=int, 3 default=64, choices=None, 4 help="number of x points") 5 getargs.add("--ny", type=int, 6 default=64, choices=None, 7 help="number of y points") 8 getargs.add("--Lx", type=float, 9 default=8., choices=None, 10 help="Lx") 11 getargs.add("--Ly", type=float, 12 default=8., choices=None, 13 help="Ly") 14 15 args = getargs.parse() 16 1 17 from dynamico.meshes import zeros 2 18 from dynamico import meshes … … 4 20 import netCDF4 as cdf 5 21 import argparse 6 7 parser = argparse.ArgumentParser()8 9 parser.add_argument("-nx", type=int,10 default=64, choices=None,11 help="number of x points")12 parser.add_argument("-ny", type=int,13 default=64, choices=None,14 help="number of y points")15 parser.add_argument("-Lx", type=float,16 default=8., choices=None,17 help="Lx")18 parser.add_argument("-Ly", type=float,19 default=8., choices=None,20 help="Ly")21 args = parser.parse_args()22 22 23 23 nx, ny, Lx, Ly, llm, nqdyn = args.nx, args.ny,args.Lx, args.Ly, 1, 1 -
codes/icosagcm/devel/Python/test/python.sh
r777 r790 20 20 NB_MPI=$1 21 21 shift 22 mpirun -np $NB_MPI python -u $* 22 rm -f xios_client*.* 23 mpirun -np $NB_MPI python -u $* | tee dynamico.log 23 24 } 24 25
Note: See TracChangeset
for help on using the changeset viewer.