- Timestamp:
- 11/05/18 11:51:07 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r765 r769 18 18 parser = argparse.ArgumentParser() 19 19 20 parser.add_argument("- mpi_ni", type=int,20 parser.add_argument("--mpi_ni", type=int, 21 21 default=64, choices=None, 22 22 help="number of x processors") 23 parser.add_argument("- mpi_nj", type=int,23 parser.add_argument("--mpi_nj", type=int, 24 24 default=64, choices=None, 25 25 help="number of y processors") 26 parser.add_argument("--T", type=float, default=5., 27 help="Length of time slice in seconds") 28 26 29 args = parser.parse_args() 27 30 … … 40 43 41 44 omega = 7.292e-5 # Angular velocity of the Earth (s^-1) 42 phi0 = 45. 45 phi0 = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg) 43 46 f0 = 2*omega*np.sin(phi0) 44 47 beta0 = 2*omega*np.cos(phi0)/a 45 48 fb = 2*omega*np.sin(phi0) - y0*2*omega*np.cos(phi0)/a 46 49 47 def Phi_xy( x,y):50 def Phi_xy(y): 48 51 fc = y*y - (Ly*y/pi)*sin(2*pi*y/Ly) 49 52 fd = Ly*Ly/(2*pi*pi)*cos(2*pi*y/Ly) 50 53 return .5*u0*( fb*(y-y0-Ly/(2*pi)*sin(2*pi*y/Ly)) + .5*beta0*(fc-fd-(Ly*Ly/3.)- Ly*Ly/(2*pi*pi)) ) 51 54 52 def Phi_xyeta( x,y,eta): return T0*g/lap*(1-eta**(Rd*lap/g)) + Phi_xy(x,y)*log(eta)*exp(-((log(eta)/b)**2))55 def Phi_xyeta(y,eta): return T0*g/lap*(1-eta**(Rd*lap/g)) + Phi_xy(y)*log(eta)*exp(-((log(eta)/b)**2)) 53 56 def ulon(x,y,eta): return -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/b/b)) 54 57 def tmean(eta) : return T0*eta**(Rd*lap/g) 55 def T( x,y,eta) : return tmean(eta)+(Phi_xy(x,y)/Rd)*(((2/(b*b))*(log(eta))**2)-1)*exp(-((0.5*log(eta))**2))58 def T(y,eta) : return tmean(eta)+(Phi_xy(y)/Rd)*(((2/(b*b))*(log(eta))**2)-1)*exp(-((0.5*log(eta))**2)) 56 59 def p(eta): return p0*eta # eta = p/p_s 57 60 … … 62 65 def coriolis(lon,lat): return f0+0.*lon 63 66 meshfile = meshes.DYNAMICO_Format(filename) 64 nqdyn, radius = 1, None65 67 pmesh = meshes.Unstructured_PMesh(comm,meshfile) 66 68 pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj) 69 nqdyn, radius = 1, None 67 70 mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 71 68 72 69 73 alpha_k = (np.arange(llm) +.5)/llm … … 86 90 print('min max eta_il', np.min(eta_il),np.max(eta_il)) 87 91 88 Phi_il = Phi_xyeta( x_il,y_il, eta_il)89 Phi_ik = Phi_xyeta( x_ik,y_ik, eta_ik)92 Phi_il = Phi_xyeta(y_il, eta_il) 93 Phi_ik = Phi_xyeta(y_ik, eta_ik) 90 94 p_ik = p(eta_ik) 91 T_ik = T( x_ik,y_ik, eta_ik) #ik full level(40), il(41)95 T_ik = T(y_ik, eta_ik) #ik full level(40), il(41) 92 96 93 97 gas = thermo.set_pT(p_ik,T_ik) … … 118 122 return thermo, mesh, params, prec.asnum([mass_ik,Sik,ujk,Phi_il,Wil]), gas 119 123 124 def diagnose(Phi,S,m,W): 125 s=S/m ; s=.5*(s+abs(s)) 126 for l in range(llm): 127 v[:,l]=(Phi[:,l+1]-Phi[:,l])/(g*m[:,l]) 128 w[:,l]=.5*params.g*(W[:,l+1]+W[:,l])/m[:,l] 129 z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g 130 gas = thermo.set_vs(v,s) 131 return gas, w, z 132 120 133 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator 121 134 comm = client.comm … … 124 137 125 138 g, Lx, Ly = 9.81, 4e7, 6e6 126 nx, ny, llm = 200, 30, 2 2139 nx, ny, llm = 200, 30, 25 127 140 dx,dy=Lx/nx,Ly/ny 128 141 … … 131 144 thermo, mesh, params, flow0, gas0 = baroclinic_3D(Lx,nx,Ly,ny,llm) 132 145 133 T, nslice, dt = 3600., 1, 3600. 134 135 with xios.Context_Curvilinear(mesh,1, dt) as context: 146 mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) 147 print 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk] 148 unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) 149 150 T=3600. 151 dt=360. 152 dz = flow0[3].max()/(params.g*llm) 153 nt = int(math.ceil(T/dt)) 154 dt = T/nt 155 print 'Time step : %d x %g = %g s' % (nt,dt,nt*dt) 156 157 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 158 159 if False: # time stepping in Python 160 caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g) 161 scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) 162 def next_flow(m,S,u,Phi,W): 163 return scheme.advance((m,S,u,Phi,W),nt) 164 165 else: # time stepping in Fortran 166 scheme = time_step.ARK2(None, dt) 167 caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, caldyn_thermo,caldyn_eta, 168 thermo,params,params.g) 169 def next_flow(m,S,u,Phi,W): 170 # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) 171 caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 172 caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W 173 caldyn_step.next() 174 return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), 175 caldyn_step.geopot.copy(), caldyn_step.W.copy()) 176 177 m,S,u,Phi,W=flow0 178 if caldyn_thermo == unst.thermo_theta: 179 s=S/m 180 theta = thermo.T0*np.exp(s/thermo.Cpd) 181 S=m*theta 182 title_format = 'Potential temperature at t=%g s (K)' 183 else: 184 title_format = 'Specific entropy at t=%g s (J/K/kg)' 185 186 w=mesh.field_mass() 187 z=mesh.field_mass() 188 xx,yy = mesh.lat_i, mesh.lon_i 189 190 191 #T, nslice, dt = 3600., 1, 3600. 192 Nslice=24 193 194 with xios.Context_Curvilinear(mesh,1, 24*3600) as context: 136 195 # now XIOS knows about the mesh and we can write to disk 137 for i in range(48): 196 v = mesh.field_mass() # specific volume (diagnosed) 197 198 for i in range(Nslice): 138 199 context.update_calendar(i) 139 print 'send_field', i, gas0.T.shape 140 # context.send_field_primal('ps', lat_i) 200 201 # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W 202 gas, w, z = diagnose(Phi,S,m,W) 203 204 # write to disk 141 205 context.send_field_primal('temp', gas0.T) 142 context.send_field_primal('p', gas0.p) 206 context.send_field_primal('p', gas.p) 207 context.send_field_primal('theta', gas.s) 208 context.send_field_primal('uz', w) 209 210 print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g') 211 #if args.mpi_ni*args.mpi_nj==1: plot() 212 213 time1, elapsed1 =time.time(), unst.getvar('elapsed') 214 m,S,u,Phi,W = next_flow(m,S,u,Phi,W) 215 time2, elapsed2 =time.time(), unst.getvar('elapsed') 216 factor = 1000./nt 217 print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 218 factor = 1e9/(4*nt*nx*ny*llm) 219 print 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 220 221 context.update_calendar(Nslice+1) # make sure XIOS writes last iteration 222 print('************DONE************')
Note: See TracChangeset
for help on using the changeset viewer.