Changeset 807 for codes/icosagcm/devel/Python/test

Timestamp:

02/08/19 18:53:34 (5 years ago)

Author:

jisesh

Message:

devel/Python : polar projection + Williamson test case 2 on LAM

Location:

codes/icosagcm/devel/Python/test/py

Files:

• RSW2_MPAS_W02.py (modified) (5 diffs)
• write_Cartesian_mesh.py (modified) (1 diff)

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

@@ -1 +1 @@
 from dynamico import getargs
+getargs.add("--LAM", action='store_true')
+# Args for global mesh
+getargs.add("--grid", type=int, help='Number of hexagons', default=2562)
+# Args for LAM
+getargs.add("--nx", type=int, help='Zonal dimension of LAM mesh', default=100)
+getargs.add("--ny", type=int, help='Meridional dimension of LAM mesh', default=100)
+getargs.add("--dx", type=float, help='Resolution at center of LAM domain', default=1e5)
+getargs.add("--center_lat", type=float, help='Latitude in degrees of LAM center', default=0.)
+getargs.add("--Davies_N1", type=int, default=3)
+getargs.add("--Davies_N2", type=int, default=3)
+
+args = getargs.parse()
+
 log_master, log_world = getargs.getLogger()
 INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
@@ -14 +27 @@
 from dynamico import unstructured as unst
 from dynamico.dev.meshes import MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh
+from dynamico.dev import meshes
 from dynamico import time_step
 from dynamico import maps
+from dynamico.LAM import Davies
+
 print '...Done'
 
@@ -24 +40 @@
 print '...Done'
 
-grid, llm, nqdyn = 2562, 1,1 # 2562, 10242, 40962
+#--------------------------- functions and classes -----------------------------
+
+class myDavies(Davies):
+    def mask(self,X,Y):
+        # X and Y are coordinates in the reference domain (cell = unit square)
+        # numerical domain extends from -nx/2 ... nx/2 and -ny/2 ... ny/2
+        # useful domain extends from -X0 ... X0 and -Y0 ... Y0
+        N3 = args.Davies_N1+args.Davies_N2
+        X0 = args.nx/2. - N3
+        Y0 = args.ny/2. - N3
+        mask  = self.mask0( X,X0,1.)  # Western boundary
+        mask *= self.mask0(-X,X0,1.) # Eastern boundary
+        mask *= self.mask0( Y,Y0,1.) # Northern boundary
+        mask *= self.mask0(-Y,Y0,1.) # Southern boundary
+        return mask
+
+#----------------------------- main program --------------------------------
+
+llm, nqdyn = 1,1 # 2562, 10242, 40962
 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4
 N, T, courant = 40, 10400., 1.2 # simulation length = N*T
 
 print 'Omega, planetary PV', Omega, 2*Omega/gh0
-planet = maps.SphereMap(radius, Omega)
 
-print 'Reading MPAS mesh ...'
-meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid)
-pmesh = PMesh(comm,meshfile)
-pmesh.partition_metis()
+if args.LAM:
+    nx, ny = args.nx, args.ny
+    filename = 'cart_%03d_%03d.nc'%(nx,ny)
+    INFO('Reading Cartesian mesh ...')
+    meshfile = meshes.DYNAMICO_Format(filename)
+    pmesh = meshes.Unstructured_PMesh(comm,meshfile)
+    pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj)
+    planet = maps.PolarStereoMap(radius,Omega, args.dx, args.center_lat*np.pi/180.)
+else:
+    planet = maps.SphereMap(radius, Omega)
+    print 'Reading MPAS mesh ...'
+    meshfile = MPAS_Format('grids/x1.%d.grid.nc'%args.grid)
+    pmesh = PMesh(comm,meshfile)
+    pmesh.partition_metis()
+
 mesh = Mesh(pmesh, llm, nqdyn, planet)
+
 print '...Done'
 lon, lat = mesh.lon_i, mesh.lat_i
@@ -53 +98 @@
 scheme = time_step.RK4(None, dt)
 print dt, scheme.csjl, scheme.cfjl
-step = unst.caldyn_step_TRSW(mesh,scheme,nt)
+step = unst.caldyn_step_TRSW(mesh,scheme,1)
 
 u0 = Omega*radius/12. # cf Williamson (1991), p.13
@@ -59 +104 @@
 print 'Williamson (1991) test 2, u0=', u0
 ulon = u0*np.cos(mesh.lat_e)
-Phi0 = gh0 - gh1*(np.sin(mesh.lat_i)**2)
-zeta0 = (2*u0/radius+2*Omega)*np.sin(mesh.lat_v)
-Phi0v = gh0 - (radius*Omega*u0+.5*u0**2)*(np.sin(mesh.lat_v)**2)
-q0 = zeta0/Phi0v
-fu_perp = mesh.ucov2D(0.*ulon,2*Omega*np.sin(mesh.lat_e)*ulon)
 
-# initial condition
-step.mass[:]=Phi0
-step.theta_rhodz[:]=Phi0
-step.u[:]=mesh.ucov2D(ulon,0.*ulon)
+gh_ini = gh0 - gh1*(np.sin(mesh.lat_i)**2)
+u_ini = mesh.ucov2D(ulon,0.*ulon)
 
-gh,u = step.mass.copy(), step.u.copy()
+# unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, u.shape[0], u) 
 
+if args.LAM:
+    davies = myDavies(args.Davies_N1, args.Davies_N2, 
+                      mesh.ref_lon_i, mesh.ref_lat_i, mesh.ref_lon_e,mesh.ref_lat_e)
+    def relax():
+        davies.relax_RSW(llm, step, (gh_ini, u_ini) )
+else:
+    def relax(): pass
+
+def next_flow(Phi, u):
+    step.mass[:], step.theta_rhodz[:], step.u[:] = Phi, Phi, u
+    for i in range(nt):
+        step.next()
+        relax()
+    return step.mass.copy(), step.u.copy()
+
+gh, u = gh_ini, u_ini
 print gh.shape, gh.min(), gh.max(), u.min(), u.max()
-flow=gh,u
 
-unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, u.shape[0], u) 
-                                  
 for i in range(20):
     if True:
-        step.next()
-        gh,u = step.mass, step.u
+        gh, u = next_flow(gh,u)
+        # step.next()
+        # gh,u = step.mass, step.u
         print i, gh.shape, gh.min(), gh.max(), u.min(), u.max()
-        mesh.plot_i(gh-Phi0) ; plt.title('err(gh)');
+        mesh.plot_i(gh-gh_ini) ; plt.title('err(gh)');
         plt.savefig('%s_err_gh_%02d.png'%(prefix,i))
         plt.close()

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

@@ -24 +24 @@
 
 dx,dy=Lx/nx,Ly/ny
-mesh = meshes.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,0.)
+mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,0.)
 print('Successfully initialized Cartesian mesh')
 mesh.ncwrite('cart_%03d_%03d.nc'%(nx,ny))