Changeset 779

Timestamp:

11/16/18 14:03:09 (5 years ago)

Author:

dubos

Message:

devel/unstructured : read dual_ne from mesh file when parallel

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/meshes.py (modified) (11 diffs)
• test/py/reorder.py (modified) (5 diffs)

codes/icosagcm/devel/Python/dynamico/meshes.py

@@ -104 +104 @@
         # Python order : nqdyn,ny,nx,llm   / indices start at 0
         # indices below follow Fortran while x,y follow Python/C
-        def periodize(x,nx): return (x+2*nx)%nx
+        def periodize(z,nz): return (z+2*nz)%nz
         def index(x,y):      return 1+periodize(x,nx)+nx*periodize(y,ny)
         def indexu(x,y):     return 2*index(x,y)-1
@@ -183 +183 @@
                          'left','right','down','up',
                          'trisk_deg','trisk','Riv2','wee',
-                         'Aiv','lon_i','lat_i','lon_v','lat_v',
+                         'lon_i','lat_i','lon_v','lat_v',
                          'le_de','le','de','lon_e','lat_e','angle_e',
                          'primal_i','primal_j','edge_i','edge_j')
+        self.Ai, self.Av = Aiv, Aiv
     def ncwrite(self, name):
         """The method writes Cartesian mesh on the disc.
@@ -223 +224 @@
                      ("primal_i","i4",self.primal_i),
                      ("primal_j","i4",self.primal_j),
-                     ("Ai","f8",self.Aiv),
+                     ("Ai","f8",self.Ai),
                      ("lon_i","f8",self.lon_i),
                      ("lat_i","f8",self.lat_i)] )
@@ -229 +230 @@
         create_vars("dual_cell", 
                     [("dual_deg","i4",self.dual_deg), 
-                     ("Av","f8",self.Aiv),
+                     ("Av","f8",self.Av),
                      ("lon_v","f8",self.lon_v),
                      ("lat_v","f8",self.lat_v),
@@ -331 +332 @@
         'lat_i':'latCell','lon_i':'lonCell','lat_v':'latVertex','lon_v':'lonVertex',
         'lat_e':'latEdge','lon_e':'lonEdge','angle_e':'angleEdge',
-        'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex'}
+        'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex',
+        'primal_ne':'signOnCell', 'dual_ne':'signOnVertex'}
     def __init__(self,gridfilename):
         try:
@@ -446 +448 @@
         trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)]
         self.set_members(locals(), 'gridfile','llm','nqdyn','radius','f',
-                         'primal_num', 'dual_num', 'edge_num', 'primal_deg', 'primal_vertex',
+                         'primal_num', 'dual_num', 'edge_num', 'primal_deg', 'primal_vertex', 'primal_ne',
                          'lat_i', 'lon_i', 'lat_v', 'lon_v', 'lat_e', 'lon_e', 'angle_e',
                          'Ai', 'Av', 'fv', 'le', 'de',
-                         'trisk_deg', 'trisk', 'wee', 'Riv2')
+                         'trisk_deg', 'trisk', 'wee', 'Riv2', 'dual_ne')
         gridfile.normalize(self, radius)
-        self.le_de = self.le/self.de
         max_primal_deg, max_dual_deg, max_trisk_deg = [ x.shape[1] for x in primal_edge, dual_edge, trisk]
         unst.init_mesh(llm, nqdyn, edge_num, primal_num,dual_num,
@@ -468 +469 @@
         dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm,
             'primal_cell','edge','dual_cell')
-        primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars(
-            dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' )
+        primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai = get_pvars(
+            dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'primal_ne', 'lon_i', 'lat_i', 'Ai' )
         edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars(
              dim_edge, 'edge_deg', 'trisk_deg', 'edge_left_right', 'edge_down_up', 
              'trisk', 'le', 'de', 'wee', 'lon_e', 'lat_e', 'angle_e')
-        dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars(
-            dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v', 'Av')
+        dual_deg, dual_vertex, dual_edge, dual_ne, Riv2, lon_v, lat_v, Av = get_pvars(
+            dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'dual_ne', 'Riv2', 'lon_v', 'lat_v', 'Av')
         if gridfile.meshtype == 'curvilinear':
             self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j')
@@ -483 +484 @@
             left_right, down_up, trisk ) : x.data = x.data-gridfile.start_index
         self.edge2cell   = Stencil_glob(edge_deg, left_right)
-        self.cell2edge   = Stencil_glob(primal_deg, primal_edge)
+        self.cell2edge   = Stencil_glob(primal_deg, primal_edge, primal_ne)
         self.cell2vertex = Stencil_glob(primal_deg, primal_vertex)
         self.edge2vertex = Stencil_glob(edge_deg, down_up)
-        self.vertex2edge = Stencil_glob(dual_deg, dual_edge)
+        self.vertex2edge = Stencil_glob(dual_deg, dual_edge, dual_ne)
         self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2)
         self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee)
@@ -579 +580 @@
         print 'TRiSK degree min/max : %d %d'%(trisk.degree.min(), trisk.degree.max() )
         print 'Riv2 degree min/max : %d %d'%(Riv2.degree.min(), Riv2.degree.max() )
-        primal_deg, primal_edge     = primal_edge.degree,   primal_edge.neigh_loc
-        dual_deg,   dual_edge       = dual_edge.degree,     dual_edge.neigh_loc
+        primal_deg, primal_edge, primal_ne = primal_edge.degree,   primal_edge.neigh_loc, primal_edge.weight
+        dual_deg,   dual_edge, dual_ne  = dual_edge.degree, dual_edge.neigh_loc, dual_edge.weight
         trisk_deg,  trisk, wee      = trisk.degree,         trisk.neigh_loc, trisk.weight
         dual_deg, dual_vertex, Riv2 = Riv2.degree,          Riv2.neigh_loc,  Riv2.weight
@@ -604 +605 @@
         trisk, dual_vertex, primal_edge, dual_edge =  trisk+1, dual_vertex+1, primal_edge+1, dual_edge+1
         left, right, down, up = left+1, right+1, down+1, up+1
-        primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right)
-        dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
 
         # store what we have computed
+        le_de = self.le/self.de
         self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num',
                          'primal_own_glo', 'primal_own_loc', 'primal_own_deg',
                          'primal_deg', 'primal_vertex', 'displ',
-                        'trisk_deg', 'trisk', 'wee', 'dual_deg', 'dual_vertex', 'Riv2')
+                         'trisk_deg', 'trisk', 'wee', 
+                         'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2',
+                         'left','right','primal_deg','primal_edge','primal_ne','le_de',
+                         'vertices_V1', 'edges_E2')
 
         # normalize and propagate info to Fortran side
@@ -620 +623 @@
                   dual_deg,dual_edge,dual_ne,dual_vertex,
                   left,right,down,up,trisk_deg,trisk,
-                  self.Ai,self.Av,self.fv,self.le/self.de,Riv2,wee)
+                  self.Ai,self.Av,self.fv,le_de,Riv2,wee)
 
         # setup halo transfers - NB : llm must be set before we call set_dynamico_transfer

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

@@ -1 +1 @@
 print 'Loading DYNAMICO modules ...'
-from dynamico.meshes import morton_index, key_to_perm, reorder_indices, reorder_values
+from dynamico.meshes import morton_index, key_to_perm, reorder_indices, reorder_values, MPAS_Format, Unstructured_Mesh
 print '...Done'
 import numpy as np
@@ -7 +7 @@
 import time
 
-grid = 10242 # 2562, 10242, 40962
+#------------ Read command-line arguments ----------#
+
+import argparse
+parser = argparse.ArgumentParser()
+parser.add_argument("-n", type=int, default=2562,
+                    help="number of primal cells in MPAS mesh")
+args = parser.parse_args()
+grid = args.n
 
 #-------------- read MPAS mesh ----------------#
 
-print 'Reading MPAS mesh ...'
-mesh=cdf.Dataset('grids_MPAS/x1.%d.grid.nc'%grid, "r")
+meshfile='grids_MPAS/x1.%d.grid.nc'%grid
+print 'Reading MPAS mesh %s ...'%meshfile
+
+def f(lon,lat): return 0.*lon
+mesh=MPAS_Format(meshfile)
+mesh=Unstructured_Mesh(mesh, 1, 1, 1., f)
+primal_ne, dual_ne = mesh.primal_ne, mesh.dual_ne
+
+mesh=cdf.Dataset(meshfile, "r")
 
 def getdims(*names): return [len(mesh.dimensions[name]) for name in names]
@@ -53 +67 @@
 for idx in primal_edge, left_right, dual_edge, down_up, primal_vertex, dual_vertex, trisk : idx[:]=idx[:]-1
 
-reorder_indices(perm_i, primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai)
-reorder_indices(perm_v, dual_vertex, dual_edge, lon_v, lat_v, Av, Riv2)
+reorder_indices(perm_i, primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai)
+reorder_indices(perm_v, dual_vertex, dual_edge, dual_ne, lon_v, lat_v, Av, Riv2)
 reorder_indices(perm_e, trisk_deg, left_right, down_up, trisk, lon_e, lat_e, angle_e, le, de, wee)
 reorder_values(perm_i, dual_vertex, left_right)
@@ -66 +80 @@
 
 mesh=cdf.Dataset('grids/x1.%d.grid.nc'%grid, "r+")
+
+# create new variables dual_ne, primal_ne
+mesh.createVariable('signOnCell','i4',('nCells','maxEdges'))
+mesh.createVariable('signOnVertex','i4',('nVertices','vertexDegree'))
 
 def putvars(datas,names): 
@@ -83 +101 @@
 putvars( (lat_e,lon_e,angle_e), ('latEdge','lonEdge','angleEdge') )
 putvars( (wee,Riv2), ('weightsOnEdge','kiteAreasOnVertex') )
+putvars( (primal_ne,dual_ne), ('signOnCell','signOnVertex') )
 
 mesh.close()