Changeset 779
- Timestamp:
- 11/16/18 14:03:09 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r773 r779 104 104 # Python order : nqdyn,ny,nx,llm / indices start at 0 105 105 # indices below follow Fortran while x,y follow Python/C 106 def periodize( x,nx): return (x+2*nx)%nx106 def periodize(z,nz): return (z+2*nz)%nz 107 107 def index(x,y): return 1+periodize(x,nx)+nx*periodize(y,ny) 108 108 def indexu(x,y): return 2*index(x,y)-1 … … 183 183 'left','right','down','up', 184 184 'trisk_deg','trisk','Riv2','wee', 185 ' Aiv','lon_i','lat_i','lon_v','lat_v',185 'lon_i','lat_i','lon_v','lat_v', 186 186 'le_de','le','de','lon_e','lat_e','angle_e', 187 187 'primal_i','primal_j','edge_i','edge_j') 188 self.Ai, self.Av = Aiv, Aiv 188 189 def ncwrite(self, name): 189 190 """The method writes Cartesian mesh on the disc. … … 223 224 ("primal_i","i4",self.primal_i), 224 225 ("primal_j","i4",self.primal_j), 225 ("Ai","f8",self.Ai v),226 ("Ai","f8",self.Ai), 226 227 ("lon_i","f8",self.lon_i), 227 228 ("lat_i","f8",self.lat_i)] ) … … 229 230 create_vars("dual_cell", 230 231 [("dual_deg","i4",self.dual_deg), 231 ("Av","f8",self.A iv),232 ("Av","f8",self.Av), 232 233 ("lon_v","f8",self.lon_v), 233 234 ("lat_v","f8",self.lat_v), … … 331 332 'lat_i':'latCell','lon_i':'lonCell','lat_v':'latVertex','lon_v':'lonVertex', 332 333 'lat_e':'latEdge','lon_e':'lonEdge','angle_e':'angleEdge', 333 'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex'} 334 'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex', 335 'primal_ne':'signOnCell', 'dual_ne':'signOnVertex'} 334 336 def __init__(self,gridfilename): 335 337 try: … … 446 448 trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)] 447 449 self.set_members(locals(), 'gridfile','llm','nqdyn','radius','f', 448 'primal_num', 'dual_num', 'edge_num', 'primal_deg', 'primal_vertex', 450 'primal_num', 'dual_num', 'edge_num', 'primal_deg', 'primal_vertex', 'primal_ne', 449 451 'lat_i', 'lon_i', 'lat_v', 'lon_v', 'lat_e', 'lon_e', 'angle_e', 450 452 'Ai', 'Av', 'fv', 'le', 'de', 451 'trisk_deg', 'trisk', 'wee', 'Riv2' )453 'trisk_deg', 'trisk', 'wee', 'Riv2', 'dual_ne') 452 454 gridfile.normalize(self, radius) 453 self.le_de = self.le/self.de454 455 max_primal_deg, max_dual_deg, max_trisk_deg = [ x.shape[1] for x in primal_edge, dual_edge, trisk] 455 456 unst.init_mesh(llm, nqdyn, edge_num, primal_num,dual_num, … … 468 469 dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, 469 470 'primal_cell','edge','dual_cell') 470 primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars(471 dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', ' lon_i', 'lat_i', 'Ai' )471 primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai = get_pvars( 472 dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'primal_ne', 'lon_i', 'lat_i', 'Ai' ) 472 473 edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars( 473 474 dim_edge, 'edge_deg', 'trisk_deg', 'edge_left_right', 'edge_down_up', 474 475 'trisk', 'le', 'de', 'wee', 'lon_e', 'lat_e', 'angle_e') 475 dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars(476 dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', ' Riv2', 'lon_v', 'lat_v', 'Av')476 dual_deg, dual_vertex, dual_edge, dual_ne, Riv2, lon_v, lat_v, Av = get_pvars( 477 dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'dual_ne', 'Riv2', 'lon_v', 'lat_v', 'Av') 477 478 if gridfile.meshtype == 'curvilinear': 478 479 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 484 485 self.edge2cell = Stencil_glob(edge_deg, left_right) 485 self.cell2edge = Stencil_glob(primal_deg, primal_edge )486 self.cell2edge = Stencil_glob(primal_deg, primal_edge, primal_ne) 486 487 self.cell2vertex = Stencil_glob(primal_deg, primal_vertex) 487 488 self.edge2vertex = Stencil_glob(edge_deg, down_up) 488 self.vertex2edge = Stencil_glob(dual_deg, dual_edge )489 self.vertex2edge = Stencil_glob(dual_deg, dual_edge, dual_ne) 489 490 self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2) 490 491 self.edge2edge = Stencil_glob(trisk_deg, trisk, wee) … … 579 580 print 'TRiSK degree min/max : %d %d'%(trisk.degree.min(), trisk.degree.max() ) 580 581 print 'Riv2 degree min/max : %d %d'%(Riv2.degree.min(), Riv2.degree.max() ) 581 primal_deg, primal_edge = primal_edge.degree, primal_edge.neigh_loc582 dual_deg, dual_edge = dual_edge.degree, dual_edge.neigh_loc582 primal_deg, primal_edge, primal_ne = primal_edge.degree, primal_edge.neigh_loc, primal_edge.weight 583 dual_deg, dual_edge, dual_ne = dual_edge.degree, dual_edge.neigh_loc, dual_edge.weight 583 584 trisk_deg, trisk, wee = trisk.degree, trisk.neigh_loc, trisk.weight 584 585 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 605 606 left, right, down, up = left+1, right+1, down+1, up+1 606 primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right)607 dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)608 607 609 608 # store what we have computed 609 le_de = self.le/self.de 610 610 self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num', 611 611 'primal_own_glo', 'primal_own_loc', 'primal_own_deg', 612 612 'primal_deg', 'primal_vertex', 'displ', 613 'trisk_deg', 'trisk', 'wee', 'dual_deg', 'dual_vertex', 'Riv2') 613 'trisk_deg', 'trisk', 'wee', 614 'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2', 615 'left','right','primal_deg','primal_edge','primal_ne','le_de', 616 'vertices_V1', 'edges_E2') 614 617 615 618 # normalize and propagate info to Fortran side … … 620 623 dual_deg,dual_edge,dual_ne,dual_vertex, 621 624 left,right,down,up,trisk_deg,trisk, 622 self.Ai,self.Av,self.fv, self.le/self.de,Riv2,wee)625 self.Ai,self.Av,self.fv,le_de,Riv2,wee) 623 626 624 627 # setup halo transfers - NB : llm must be set before we call set_dynamico_transfer -
codes/icosagcm/devel/Python/test/py/reorder.py
r674 r779 1 1 print 'Loading DYNAMICO modules ...' 2 from dynamico.meshes import morton_index, key_to_perm, reorder_indices, reorder_values 2 from dynamico.meshes import morton_index, key_to_perm, reorder_indices, reorder_values, MPAS_Format, Unstructured_Mesh 3 3 print '...Done' 4 4 import numpy as np … … 7 7 import time 8 8 9 grid = 10242 # 2562, 10242, 40962 9 #------------ Read command-line arguments ----------# 10 11 import argparse 12 parser = argparse.ArgumentParser() 13 parser.add_argument("-n", type=int, default=2562, 14 help="number of primal cells in MPAS mesh") 15 args = parser.parse_args() 16 grid = args.n 10 17 11 18 #-------------- read MPAS mesh ----------------# 12 19 13 print 'Reading MPAS mesh ...' 14 mesh=cdf.Dataset('grids_MPAS/x1.%d.grid.nc'%grid, "r") 20 meshfile='grids_MPAS/x1.%d.grid.nc'%grid 21 print 'Reading MPAS mesh %s ...'%meshfile 22 23 def f(lon,lat): return 0.*lon 24 mesh=MPAS_Format(meshfile) 25 mesh=Unstructured_Mesh(mesh, 1, 1, 1., f) 26 primal_ne, dual_ne = mesh.primal_ne, mesh.dual_ne 27 28 mesh=cdf.Dataset(meshfile, "r") 15 29 16 30 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 54 68 55 reorder_indices(perm_i, primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai)56 reorder_indices(perm_v, dual_vertex, dual_edge, lon_v, lat_v, Av, Riv2)69 reorder_indices(perm_i, primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai) 70 reorder_indices(perm_v, dual_vertex, dual_edge, dual_ne, lon_v, lat_v, Av, Riv2) 57 71 reorder_indices(perm_e, trisk_deg, left_right, down_up, trisk, lon_e, lat_e, angle_e, le, de, wee) 58 72 reorder_values(perm_i, dual_vertex, left_right) … … 66 80 67 81 mesh=cdf.Dataset('grids/x1.%d.grid.nc'%grid, "r+") 82 83 # create new variables dual_ne, primal_ne 84 mesh.createVariable('signOnCell','i4',('nCells','maxEdges')) 85 mesh.createVariable('signOnVertex','i4',('nVertices','vertexDegree')) 68 86 69 87 def putvars(datas,names): … … 83 101 putvars( (lat_e,lon_e,angle_e), ('latEdge','lonEdge','angleEdge') ) 84 102 putvars( (wee,Riv2), ('weightsOnEdge','kiteAreasOnVertex') ) 103 putvars( (primal_ne,dual_ne), ('signOnCell','signOnVertex') ) 85 104 86 105 mesh.close()
Note: See TracChangeset
for help on using the changeset viewer.