Changeset 701 for codes/icosagcm/devel
- Timestamp:
- 05/24/18 15:00:36 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/DCMIP2008c5_MPI.py
r697 r701 1 # apparently we should initialize MPI before doing anything else 2 from mpi4py import MPI 3 comm = MPI.COMM_WORLD 4 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 5 print '%d/%d starting'%(mpi_rank,mpi_size) 6 7 # now start doing something useful 1 8 from dynamico import unstructured as unst 2 9 from dynamico import dyn … … 4 11 from dynamico import DCMIP 5 12 from dynamico import meshes 6 import dynamico.xios as xios13 #import dynamico.xios as xios 7 14 8 15 import math as math … … 10 17 import numpy as np 11 18 import time 12 13 from mpi4py import MPI 14 comm = MPI.COMM_WORLD 15 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 16 print '%d/%d starting'%(mpi_rank,mpi_size) 19 import argparse 17 20 18 21 #------------------------ initial condition ------------------------- … … 85 88 0.07810627, 0.1069411, 0.1408637, 0.1807720, 0.2277220, 0.2829562, 0.3479364, 0.4243822, 86 89 0.5143168, 0.6201202, 0.7235355, 0.8176768, 0.8962153, 0.9534761, 0.9851122, 1.0 ] 90 if llm==49: 91 ap_l=[0.002251865, 0.003983890, 0.006704364, 0.01073231, 0.01634233, 0.02367119, 92 0.03261456, 0.04274527, 0.05382610, 0.06512175, 0.07569850, 0.08454283, 93 0.08396310, 0.08334103, 0.08267352, 0.08195725, 0.08118866, 0.08036393, 94 0.07947895, 0.07852934, 0.07751036, 0.07641695, 0.07524368, 0.07398470, 95 0.07263375, 0.07118414, 0.06962863, 0.06795950, 0.06616846, 0.06424658, 96 0.06218433, 0.05997144, 0.05759690, 0.05504892, 0.05231483, 0.04938102, 97 0.04623292, 0.04285487, 0.03923006, 0.03534049, 0.03116681, 0.02668825, 98 0.02188257, 0.01676371, 0.01208171, 0.007959612, 0.004510297, 0.001831215, 99 0.0, 0.0 ] 100 bp_l=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 101 0.006755112, 0.01400364, 0.02178164, 0.03012778, 0.03908356, 0.04869352, 102 0.05900542, 0.07007056, 0.08194394, 0.09468459, 0.1083559, 0.1230258, 103 0.1387673, 0.1556586, 0.1737837, 0.1932327, 0.2141024, 0.2364965, 104 0.2605264, 0.2863115, 0.3139801, 0.3436697, 0.3755280, 0.4097133, 105 0.4463958, 0.4857576, 0.5279946, 0.5733168, 0.6219495, 0.6741346, 106 0.7301315, 0.7897776, 0.8443334, 0.8923650, 0.9325572, 0.9637744, 107 0.9851122, 1.0] 87 108 88 109 ap_l, bp_l = ref_press*np.asarray(ap_l[-1::-1]), bp_l[-1::-1] # reverse indices 89 110 ptop = ap_l[-1] # pressure BC 90 111 91 print ptop, ap_l, bp_l112 if mpi_rank==0: print ptop, ap_l, bp_l 92 113 ps_il,ap_il = np.meshgrid(ps_i,ap_l, indexing='ij') 93 114 ps_il,bp_il = np.meshgrid(ps_i,bp_l, indexing='ij') … … 109 130 ujk, Wil = mesh.ucov3D(ulon(lat_ek),ulat(lat_ek)), mesh.field_w() 110 131 111 print 'ztop (m) = ', Phi_il[0,-1]/g 112 print 'ptop (Pa) = ', gas.p[0,-1], ptop 132 if mpi_rank==0: 133 print 'ztop (m) = ', Phi_il[0,-1]/g 134 print 'ptop (Pa) = ', gas.p[0,-1], ptop 113 135 dx=mesh.de.min() 114 136 params=dyn.Struct() … … 118 140 return thermo, mesh, hybrid_coefs, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas 119 141 120 121 142 #------------------------ main program ------------------------- 122 143 123 #grid, llm = 40962, 26 124 nqtot, llm, grid = 1, 26, 2562 144 unst.setvar('dynamico_mpi_rank', mpi_rank) 145 146 parser = argparse.ArgumentParser() 147 parser.add_argument("-r", "--refinement", type=int, 148 default=5, choices=[4, 5, 6, 7], 149 help="grid refinement level") 150 parser.add_argument("-l", "--llm", type=int, 151 default=49, choices=[18, 26, 49], 152 help="number of vertical levels") 153 args = parser.parse_args() 154 nqtot, llm, grid = 1, args.llm, 2+10*(4**args.refinement) 155 #nqtot, llm, grid = 1,26,40962 156 125 157 T, Nslice, courant = 14400, 24, 3.0 126 158 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag 127 159 #caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 160 128 161 thermo, mesh, hybrid_coefs, params, flow0, gas0 = DCMIP2008c5(grid,llm) 129 162 llm, dx = mesh.llm, params.dx 130 print 'llm, nb_hex, dx =', llm, mesh.Ai.size, dx 131 if caldyn_eta == unst.eta_lag: 132 print 'Lagrangian coordinate.' 163 if mpi_rank==0: 164 print 'grid, llm, local_gridsize, dx =', grid, llm, mesh.Ai.size, dx 165 if caldyn_eta == unst.eta_lag: 166 print 'Lagrangian coordinate.' 167 else: 168 print 'Mass-based coordinate.' 169 170 unst.ker.dynamico_init_hybrid(*hybrid_coefs) 171 172 dt = courant*.5*dx/np.sqrt(gas0.c2.max()) 173 174 if False: 175 nt = int(math.ceil(T/dt)) 176 dt = T/nt 133 177 else: 134 print 'Mass-based coordinate.' 135 136 unst.ker.dynamico_init_hybrid(*hybrid_coefs) 137 138 dt = courant*.5*dx/np.sqrt(gas0.c2.max()) 139 140 nt = int(math.ceil(T/dt)) 141 dt = T/nt 142 print 'Time step : %g s' % dt 178 nt=100 179 if mpi_rank==0: print 'Time step : %g s' % dt 143 180 144 181 #mesh.plot_e(mesh.le/mesh.de) ; plt.title('le/de') … … 150 187 scheme = time_step.ARK2(None, dt, a32=0.7) 151 188 caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta, thermo,params,params.g) 189 152 190 def next_flow(m,S,u): 153 191 caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 154 caldyn_step.remap()192 # caldyn_step.remap() 155 193 caldyn_step.next() 156 194 return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), … … 186 224 #plots(0) 187 225 188 context=xios.XIOS_Context(mesh.pmesh,mesh,nqtot, T)226 #context=xios.XIOS_Context(mesh.pmesh,mesh,nqtot, T) 189 227 190 228 for it in range(Nslice): 191 context.update_calendar(it)229 # context.update_calendar(it) 192 230 unst.setvar('debug_hevi_solver',False) 193 231 time1, elapsed1 =time.time(), unst.getvar('elapsed') 194 232 m,S,u,Phi = next_flow(m,S,u) 195 233 time2, elapsed2 = time.time(), unst.getvar('elapsed') 196 factor = 1000./nt 197 print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 198 factor = 1e9/(4*nt*m.size) 199 print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 200 201 s=S/m 202 s=.5*(s+abs(s)) 203 for l in range(llm): 204 z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g 205 vol[:,l]=(Phi[:,l+1]-Phi[:,l])/params.g/m[:,l] # specific volume 206 gas = thermo.set_vs(vol, s) 207 ss = np.asarray(gas.T, dtype=np.double) 208 context.send_field_primal('theta', ss) 209 #plots(it+1) 210 211 print 'xios.context_finalize()' 212 context.finalize() 213 print 'xios.finalize()' 214 xios.finalize() 215 print 'Done' 234 if mpi_rank==0: 235 factor = 1000./nt 236 print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 237 factor = 1e9/(4*nt*m.size) 238 print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1) 239 240 if False: 241 s=S/m 242 s=.5*(s+abs(s)) 243 for l in range(llm): 244 z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g 245 vol[:,l]=(Phi[:,l+1]-Phi[:,l])/params.g/m[:,l] # specific volume 246 gas = thermo.set_vs(vol, s) 247 ss = np.asarray(gas.T, dtype=np.double) 248 # context.send_field_primal('theta', ss) 249 #plots(it+1) 250 251 unst.ker.dynamico_print_trace() 252 253 #print 'xios.context_finalize()' 254 #context.finalize() 255 #print 'xios.finalize()' 256 #xios.finalize() 257 print 'MPI Rank %d Done'%mpi_rank
Note: See TracChangeset
for help on using the changeset viewer.