- Timestamp:
- 11/14/18 16:54:54 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r771 r774 12 12 import time 13 13 import argparse 14 import logging 14 15 from numpy import pi, log, exp, sin, cos 15 16 … … 22 23 parser.add_argument("--mpi_nj", type=int, default=1, 23 24 help="number of y processors") 24 parser.add_argument("--T", type=float, default= 5.,25 parser.add_argument("--T", type=float, default=3600., 25 26 help="Length of time slice in seconds") 26 parser.add_argument("--Davies_N1", type=int, default=5) 27 parser.add_argument("--Davies_N2", type=int, default=5) 27 parser.add_argument("--N", type=int, default=24, 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) 28 33 29 34 … … 42 47 Cpd = 1004.5 43 48 p0 = 1e5 49 up = 1. # amplitude (m/s) 50 xc,yc,lp = 0.,Ly/2.,600000. 44 51 45 52 omega = 7.292e-5 # Angular velocity of the Earth (s^-1) … … 55 62 56 63 def Phi_xyeta(y,eta): return T0*g/lap*(1-eta**(Rd*lap/g)) + Phi_xy(y)*log(eta)*exp(-((log(eta)/b)**2)) 57 def ulon(x,y,eta): return -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/b/b)) 64 def ulon(x,y,eta): 65 u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 66 # u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) 67 return u 68 58 69 def tmean(eta) : return T0*eta**(Rd*lap/g) 59 70 def T(y,eta) : return tmean(eta)+(Phi_xy(y)/Rd)*(((2/(b*b))*(log(eta))**2)-1)*exp(-((0.5*log(eta))**2)) … … 61 72 62 73 def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels 74 def coriolis(x,y): return f0+beta0*y 63 75 64 76 filename = 'cart_%03d_%03d.nc'%(nx,ny) 65 print 'Reading Cartesian mesh ...' 66 def coriolis(x,y): return f0+beta0*(y+.5*Ly) 77 logging.info('Reading Cartesian mesh ...') 67 78 meshfile = meshes.DYNAMICO_Format(filename) 68 79 pmesh = meshes.Unstructured_PMesh(comm,meshfile) … … 74 85 alpha_k = (np.arange(llm) +.5)/llm 75 86 alpha_l = (np.arange(llm+1)+ 0.)/llm 76 x_ik, alpha_ik = np.meshgrid(mesh.lon_i, alpha_k, indexing='ij') 77 y_ik, alpha_ik = np.meshgrid(mesh.lat_i, alpha_k, indexing='ij') 78 x_il, alpha_il = np.meshgrid(mesh.lon_i, alpha_l, indexing='ij') 79 y_il, alpha_il = np.meshgrid(mesh.lat_i, alpha_l, indexing='ij') 80 x_ek, alpha_ek = np.meshgrid(mesh.lon_e, alpha_k, indexing='ij') 81 y_ek, alpha_ek = np.meshgrid(mesh.lat_e, alpha_k, indexing='ij') 82 83 print('----------------') 84 print 'ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g)) 87 x_ik, alpha_ik = np.meshgrid(mesh.lon_i, alpha_k, indexing='ij') 88 y_ik, alpha_ik = np.meshgrid(mesh.lat_i+.5*Ly, alpha_k, indexing='ij') 89 x_il, alpha_il = np.meshgrid(mesh.lon_i, alpha_l, indexing='ij') 90 y_il, alpha_il = np.meshgrid(mesh.lat_i+.5*Ly, alpha_l, indexing='ij') 91 x_ek, alpha_ek = np.meshgrid(mesh.lon_e, alpha_k, indexing='ij') 92 y_ek, alpha_ek = np.meshgrid(mesh.lat_e+.5*Ly, alpha_k, indexing='ij') 93 94 print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g))) 85 95 print(np.shape(alpha_k),np.shape(alpha_l)) 86 96 thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0) … … 100 110 for l in range(llm): 101 111 mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l]) 102 Sik, ujk, Wil = gas.s*mass_ik, mesh.field_u(), mesh.field_w()112 Sik, Wil = gas.s*mass_ik, mesh.field_w() 103 113 104 114 u_ek = mesh.ucov3D(ulon(x_ek, y_ek, eta_ek), 0.*eta_ek) 105 115 106 print 'ztop (m) = ', Phi_il[0,-1]/g, ztop116 print('ztop (m) = ', Phi_il[0,-1]/g, ztop) 107 117 ptop = p(eta(1.)) 108 print 'ptop (Pa) = ', gas.p[0,-1], ptop118 print( 'ptop (Pa) = ', gas.p[0,-1], ptop) 109 119 110 120 params=dyn.Struct() … … 116 126 # define parameters for lower BC 117 127 pbot = p(eta_il[:,0]) 118 print 'min p, T :', pbot.min(), tmean(pbot/p0)128 print( 'min p, T :', pbot.min(), tmean(pbot/p0)) 119 129 gas_bot = thermo.set_pT(pbot, tmean(pbot/p0)) 120 130 params.pbot = gas_bot.p 121 131 params.rho_bot = 1e6/gas_bot.v 122 132 123 return thermo, mesh, params, prec.asnum([mass_ik,Sik,u jk,Phi_il,Wil]), gas133 return thermo, mesh, params, prec.asnum([mass_ik,Sik,u_ek,Phi_il,Wil]), gas 124 134 125 135 def diagnose(Phi,S,m,W): … … 137 147 self.beta_i = self.mask(x_i,y_i) 138 148 self.beta_e = self.mask(x_e,y_e) 149 self.x_e,self.y_e = x_e,y_e 139 150 def mask0(self,c,c0,delta): # 1D building block for Davies relaxation 140 151 N1, N2 = self.N1, self.N2 141 N3=N1+N2142 152 m = np.zeros(c.size) 143 153 c1,c2 = c0-N1*delta, c0-(N1+N2)*delta 144 154 for i in range(c.size): 145 155 ci=c[i] 146 m[i] = (1.+np.cos((ci-c 0+N3*delta)*np.pi/(N2*delta)))/2.0147 if ci < c 0-N3*delta: m[i]=1.148 if ci > c 0-N1*delta: m[i]=0.156 m[i] = (1.+np.cos((ci-c2)*np.pi/(N2*delta)))/2.0 157 if ci < c2 : m[i]=1. 158 if ci > c1 : m[i]=0. 149 159 return m 150 160 def relax(self, llm, step, flow): … … 159 169 step.W[:,l] = beta_i*step.W[:,l] + (1.-beta_i)*W[:,l] 160 170 161 162 171 class myDavies(Davies): 163 172 def mask(self,x,y): 164 # return self.mask0(y,Ly,dy) 165 return self.mask0(y,-.5*Ly,dy)*self.mask0(-y,-.5*Ly,dy) 173 logging.debug('davies dy = %f' % dy) 174 return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy) 175 176 166 177 167 178 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator 168 179 comm = client.comm 169 180 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 170 print '%d/%d starting'%(mpi_rank,mpi_size) 181 182 #-------------Logging verbosity configuration--------------------------- 183 myloglevel = args.log 184 myloglevel = getattr(logging, myloglevel.upper()) 185 logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, 186 format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' ) 187 188 logging.debug('%d/%d starting'%(mpi_rank,mpi_size)) 171 189 172 190 g, Lx, Ly = 9.81, 4e7, 6e6 173 191 nx, ny, llm = 200, 30, 25 174 dx,dy =Lx/nx,Ly/ny192 dx,dy = Lx/nx,Ly/ny 175 193 176 194 unst.setvar('g',g) … … 179 197 180 198 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 181 print 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk]199 print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk]) 182 200 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 183 201 184 T =3600.185 dt =360.202 T = args.T 203 dt = 360. 186 204 dz = flow0[3].max()/(params.g*llm) 187 205 nt = int(math.ceil(T/dt)) 188 206 dt = T/nt 189 print 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)207 logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) 190 208 191 209 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass … … 201 219 caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 202 220 davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 221 print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) 222 logging.debug('mask min/max :') 223 logging.debug(davies.beta_i.min()) 224 logging.debug(davies.beta_i.max()) 203 225 def next_flow(m,S,u,Phi,W): 204 226 # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) … … 211 233 caldyn_step.geopot.copy(), caldyn_step.W.copy()) 212 234 213 m,S,u,Phi,W =flow0235 m,S,u,Phi,W = flow0 214 236 if caldyn_thermo == unst.thermo_theta: 215 237 s=S/m … … 225 247 226 248 #T, nslice, dt = 3600., 1, 3600. 227 Nslice =24249 Nslice = args.N 228 250 229 251 with xios.Context_Curvilinear(mesh,1, 24*3600) as context: … … 238 260 239 261 # write to disk 262 context.send_field_primal('ps', davies.beta_i) 240 263 context.send_field_primal('temp', gas0.T) 241 264 context.send_field_primal('p', gas.p) … … 243 266 context.send_field_primal('uz', w) 244 267 245 print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g') 246 #if args.mpi_ni*args.mpi_nj==1: plot() 268 print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')) 247 269 248 270 time1, elapsed1 =time.time(), unst.getvar('elapsed') … … 250 272 time2, elapsed2 =time.time(), unst.getvar('elapsed') 251 273 factor = 1000./nt 252 print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)274 print( 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)) 253 275 factor = 1e9/(4*nt*nx*ny*llm) 254 print 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)276 print( 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)) 255 277 256 278 context.update_calendar(Nslice+1) # make sure XIOS writes last iteration 257 print('************DONE************')279 logging.info('************DONE************')
Note: See TracChangeset
for help on using the changeset viewer.