Changeset 701

Timestamp:

05/24/18 15:00:36 (6 years ago)

Author:

dubos

Message:

devel/unstructured : use command line parameters for test case DCMIP2008c5

File:

• codes/icosagcm/devel/Python/test/py/DCMIP2008c5_MPI.py (modified) (8 diffs)

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

@@ +1 @@
+# apparently we should initialize MPI before doing anything else
+from mpi4py import MPI
+comm = MPI.COMM_WORLD
+mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
+print '%d/%d starting'%(mpi_rank,mpi_size)
+
+# now start doing something useful
 from dynamico import unstructured as unst
 from dynamico import dyn
@@ -4 +11 @@
 from dynamico import DCMIP
 from dynamico import meshes
-import dynamico.xios as xios
+#import dynamico.xios as xios
 
 import math as math
@@ -10 +17 @@
 import numpy as np
 import time
-
-from mpi4py import MPI
-comm = MPI.COMM_WORLD
-mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
-print '%d/%d starting'%(mpi_rank,mpi_size)
+import argparse
 
 #------------------------ initial condition -------------------------
@@ -85 +88 @@
                 0.07810627, 0.1069411, 0.1408637, 0.1807720, 0.2277220, 0.2829562, 0.3479364, 0.4243822, 
                 0.5143168, 0.6201202, 0.7235355, 0.8176768, 0.8962153, 0.9534761, 0.9851122, 1.0 ]
+    if llm==49:
+        ap_l=[0.002251865, 0.003983890, 0.006704364, 0.01073231, 0.01634233, 0.02367119, 
+              0.03261456, 0.04274527, 0.05382610, 0.06512175, 0.07569850, 0.08454283, 
+              0.08396310, 0.08334103, 0.08267352, 0.08195725, 0.08118866, 0.08036393, 
+              0.07947895, 0.07852934, 0.07751036, 0.07641695, 0.07524368, 0.07398470, 
+              0.07263375, 0.07118414, 0.06962863, 0.06795950, 0.06616846, 0.06424658, 
+              0.06218433, 0.05997144, 0.05759690, 0.05504892, 0.05231483, 0.04938102, 
+              0.04623292, 0.04285487, 0.03923006, 0.03534049, 0.03116681, 0.02668825, 
+              0.02188257, 0.01676371, 0.01208171, 0.007959612, 0.004510297, 0.001831215, 
+              0.0, 0.0 ]
+        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, 
+              0.006755112, 0.01400364, 0.02178164, 0.03012778, 0.03908356, 0.04869352, 
+              0.05900542, 0.07007056, 0.08194394, 0.09468459, 0.1083559, 0.1230258, 
+              0.1387673, 0.1556586, 0.1737837, 0.1932327, 0.2141024, 0.2364965, 
+              0.2605264, 0.2863115, 0.3139801, 0.3436697, 0.3755280, 0.4097133, 
+              0.4463958, 0.4857576, 0.5279946, 0.5733168, 0.6219495, 0.6741346, 
+              0.7301315, 0.7897776, 0.8443334, 0.8923650, 0.9325572, 0.9637744, 
+              0.9851122, 1.0]
 
     ap_l, bp_l = ref_press*np.asarray(ap_l[-1::-1]), bp_l[-1::-1] # reverse indices
     ptop = ap_l[-1] # pressure BC
     
-    print ptop, ap_l, bp_l
+    if mpi_rank==0: print ptop, ap_l, bp_l
     ps_il,ap_il = np.meshgrid(ps_i,ap_l, indexing='ij')
     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()
     
-    print 'ztop (m) = ', Phi_il[0,-1]/g
-    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    if mpi_rank==0: 
+        print 'ztop (m) = ', Phi_il[0,-1]/g
+        print 'ptop (Pa) = ', gas.p[0,-1], ptop
     dx=mesh.de.min()
     params=dyn.Struct()
@@ -118 +140 @@
     return thermo, mesh, hybrid_coefs, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
 
-
 #------------------------ main program -------------------------
 
-#grid, llm = 40962, 26
-nqtot, llm, grid = 1, 26, 2562
+unst.setvar('dynamico_mpi_rank', mpi_rank)
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-r", "--refinement", type=int, 
+                    default=5, choices=[4, 5, 6, 7],
+                    help="grid refinement level")
+parser.add_argument("-l", "--llm", type=int, 
+                    default=49, choices=[18, 26, 49],
+                    help="number of vertical levels")
+args = parser.parse_args()
+nqtot, llm, grid = 1, args.llm, 2+10*(4**args.refinement)
+#nqtot, llm, grid = 1,26,40962
+
 T, Nslice, courant = 14400, 24, 3.0
 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag
 #caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
+
 thermo, mesh, hybrid_coefs, params, flow0, gas0 = DCMIP2008c5(grid,llm)
 llm, dx = mesh.llm, params.dx
-print 'llm, nb_hex, dx =', llm, mesh.Ai.size, dx
-if caldyn_eta == unst.eta_lag:
-    print 'Lagrangian coordinate.'
+if mpi_rank==0: 
+    print 'grid, llm, local_gridsize, dx =', grid, llm, mesh.Ai.size, dx
+    if caldyn_eta == unst.eta_lag:
+        print 'Lagrangian coordinate.'
+    else:
+        print 'Mass-based coordinate.'
+
+unst.ker.dynamico_init_hybrid(*hybrid_coefs)
+    
+dt = courant*.5*dx/np.sqrt(gas0.c2.max())
+
+if False:
+    nt = int(math.ceil(T/dt))
+    dt = T/nt
 else:
-    print 'Mass-based coordinate.'
-
-unst.ker.dynamico_init_hybrid(*hybrid_coefs)
-    
-dt = courant*.5*dx/np.sqrt(gas0.c2.max())
-
-nt = int(math.ceil(T/dt))
-dt = T/nt
-print 'Time step : %g s' % dt
+    nt=100
+if mpi_rank==0: print 'Time step : %g s' % dt
 
 #mesh.plot_e(mesh.le/mesh.de) ; plt.title('le/de')
@@ -150 +187 @@
 scheme = time_step.ARK2(None, dt, a32=0.7)
 caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta, thermo,params,params.g)
+
 def next_flow(m,S,u):
     caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
-    caldyn_step.remap()
+#    caldyn_step.remap()
     caldyn_step.next()
     return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), 
@@ -186 +224 @@
 #plots(0)
 
-context=xios.XIOS_Context(mesh.pmesh,mesh,nqtot, T)
+#context=xios.XIOS_Context(mesh.pmesh,mesh,nqtot, T)
 
 for it in range(Nslice):
-    context.update_calendar(it)
+#    context.update_calendar(it)
     unst.setvar('debug_hevi_solver',False)
     time1, elapsed1 =time.time(), unst.getvar('elapsed')
     m,S,u,Phi = next_flow(m,S,u)
     time2, elapsed2 = time.time(), unst.getvar('elapsed')
-    factor = 1000./nt
-    print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
-    factor = 1e9/(4*nt*m.size)
-    print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
-    
-    s=S/m
-    s=.5*(s+abs(s))
-    for l in range(llm):
-        z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g
-        vol[:,l]=(Phi[:,l+1]-Phi[:,l])/params.g/m[:,l] # specific volume
-    gas = thermo.set_vs(vol, s)
-    ss = np.asarray(gas.T, dtype=np.double)
-    context.send_field_primal('theta', ss)
-    #plots(it+1)
-
-print 'xios.context_finalize()'
-context.finalize()
-print 'xios.finalize()'
-xios.finalize()
-print 'Done'
+    if mpi_rank==0: 
+        factor = 1000./nt
+        print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+        factor = 1e9/(4*nt*m.size)
+        print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+    
+    if False:
+        s=S/m
+        s=.5*(s+abs(s))
+        for l in range(llm):
+            z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g
+            vol[:,l]=(Phi[:,l+1]-Phi[:,l])/params.g/m[:,l] # specific volume
+        gas = thermo.set_vs(vol, s)
+        ss = np.asarray(gas.T, dtype=np.double)
+        #    context.send_field_primal('theta', ss)
+        #plots(it+1)
+
+unst.ker.dynamico_print_trace()
+ 
+#print 'xios.context_finalize()'
+#context.finalize()
+#print 'xios.finalize()'
+#xios.finalize()
+print 'MPI Rank %d Done'%mpi_rank