Changeset 680 for codes/icosagcm/devel
- Timestamp:
- 02/09/18 16:24:33 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 1 added
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r674 r680 5 5 import matplotlib.tri as tri 6 6 import matplotlib.pyplot as plt 7 from unstructured import init_mesh, list_stencil, ker 7 from matplotlib.patches import Polygon 8 from matplotlib.collections import PatchCollection 9 10 import unstructured as unst 11 import parallel 12 from util import list_stencil, Base_class 8 13 9 14 radian=180/math.pi … … 35 40 mass_dbk[i,j,:] = dbk_ij 36 41 mass_bl[i,j,1:]= 1.-dbk_ij.cumsum() 37 42 38 43 return mass_bl, mass_dak, mass_dbk 44 45 # from theta.k90 : rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij) 46 # from caldyn_wflux.k90: wflux(l,ij) = mass_bl(l,ij) * dmass_col(ij) - convm(l,ij) 47 def mass_coefs_from_pressure_coefs(g, ap_il, bp_il): 48 # p = Mg + ptop = a + b.ps 49 # ps = Mcol.g+ptop 50 # M = mass_a + mass_b.Mcol 51 # M = (a+b.ps-ptop)/g 52 # = (a+b.ptop-ptop)/g + b.Mcol 53 # mass_a = (a+b.ptop)/g 54 # mass_da = (da+db.ptop)/g 55 # mass_b = b 56 # mass_db = db 57 n, llm, ptop = ap_il.shape[0], ap_il.shape[1]-1, ap_il[0,-1] 58 print 'hybrid ptop : ', ptop 59 mass_al = (ap_il+ptop*bp_il)/g 60 mass_dak, mass_dbk = np.zeros((n,llm)), np.zeros((n,llm)) 61 for k in range(llm): 62 mass_dak[:,k] = mass_al[:,k]-mass_al[:,k+1] 63 mass_dbk[:,k] = bp_il[:,k]-bp_il[:,k+1] 64 print mass_al[0,:] 65 print mass_dak[0,:] 66 print mass_dbk[0,:] 67 return [ np.asarray(x) for x in bp_il, mass_dak, mass_dbk ] 39 68 40 69 #----------------------- Cartesian mesh ----------------------- … … 142 171 (indexu(x-1,y+1),.25))) 143 172 Aiv=np.zeros(nx*ny)+dx*dy # Ai=Av=dx*dy 144 init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4,173 unst.init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4, 145 174 primal_nb,primal_edge,primal_ne, 146 175 dual_nb,dual_edge,dual_ne,dual_vertex, … … 164 193 #---------------------- MPAS fully unstructured mesh ------------------------ 165 194 166 def compute_ne(num,deg,edges,left,right): 167 ne = np.zeros_like(edges) 168 for cell in range(num): 169 for iedge in range(deg[cell]): 170 edge = edges[cell,iedge]-1 171 if left[edge]==cell+1: ne[cell,iedge]+=1 172 if right[edge]==cell+1: ne[cell,iedge]-=1 173 if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge 174 return ne 175 176 def plot(tri,data, **kwargs): 177 plt.figure(figsize=(12,4)) 178 plt.gca().set_aspect('equal') 179 plt.tricontourf(tri, data, 20, **kwargs) 180 plt.colorbar() 181 plt.ylim((-90,90)) 182 plt.xlim((0,360)) 183 184 class MPAS_Mesh: 185 def __init__(self, gridfile, llm, nqdyn, radius, f): 186 self.gridfile, self.llm, self.nqdyn = gridfile,llm,nqdyn 187 self.radius, self.f = radius, f 188 # open mesh file, get main dimensions 195 class Mesh_Format(Base_class): 196 def getdims(self,*names): return [len(self.nc.dimensions[name]) for name in names] 197 def get_pdims(self,comm,*names): return [self.get_pdim(comm,name) for name in names] 198 def get_pvars(self,pdim,*names): return [self.get_pvar(pdim,name) for name in names] 199 def getvars(self,*names): 200 for name in names : print "getvar %s ..."%name 201 time1=time.time() 202 ret=[self.getvar(name) for name in names] 203 print "... Done in %f seconds"%(time.time()-time1) 204 return ret 205 206 class MPAS_Format(Mesh_Format): 207 start_index=1 # indexing starts at 1 as in Fortran 208 translate= { 209 'primal_num':'nCells', 210 'edge_num':'nEdges', 211 'dual_num':'nVertices', 212 'primal_deg':'nEdgesOnCell', 'trisk_deg':'nEdgesOnEdge', 213 'primal_edge':'edgesOnCell', 'left_right':'cellsOnEdge', 214 'dual_edge':'edgesOnVertex', 'down_up':'verticesOnEdge', 215 'primal_vertex':'verticesOnCell', 'dual_vertex':'cellsOnVertex', 216 'trisk':'edgesOnEdge', 'down_up':'verticesOnEdge', 'left_right':'cellsOnEdge', 217 'le':'dvEdge', 'de':'dcEdge', 'Ai':'areaCell', 'Av':'areaTriangle', 218 'lat_i':'latCell','lon_i':'lonCell','lat_v':'latVertex','lon_v':'lonVertex', 219 'lat_e':'latEdge','lon_e':'lonEdge','angle_e':'angleEdge', 220 'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex'} 221 def __init__(self,gridfilename): 189 222 try: 190 nc = cdf.Dataset(gridfile, "r")223 self.nc = cdf.Dataset(gridfilename, "r") 191 224 except RuntimeError: 192 225 print """ … … 195 228 """ % gridfile 196 229 raise 197 198 def getdims(*names): return [len(nc.dimensions[name]) for name in names] 199 def getvars(*names): 200 for name in names : print "getvar %s ..."%name 201 time1=time.time() 202 ret=[nc.variables[name][:] for name in names] 203 print "... Done in %f seconds"%(time.time()-time1) 204 return ret 205 primal_num, edge_num, dual_num = getdims('nCells','nEdges','nVertices') 206 print 'Number of primal cells, dual cells and edges : %d, %d, %d '%(primal_num,dual_num,edge_num) 207 primal_deg, trisk_deg = getvars('nEdgesOnCell','nEdgesOnEdge') 208 dual_deg = [3 for i in range(dual_num) ] 209 dual_deg = np.ascontiguousarray(dual_deg,dtype=np.int32) # NB : Fortran code expects 32-bit ints 210 211 # get indices for stencils 212 # primal <-> edges 213 primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge') 214 left,right = left_right[:,0], left_right[:,1] 215 primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right) 216 # dual <-> edges 217 dual_edge, down_up = getvars('edgesOnVertex','verticesOnEdge') 218 down,up = down_up[:,0], down_up[:,1] 219 dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down) 220 # primal <-> dual, edges <-> edges (TRiSK) 221 primal_vertex, dual_vertex, trisk = getvars('verticesOnCell','cellsOnVertex','edgesOnEdge') 222 # get positions, lengths, surfaces and weights 223 le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle') 224 lat_i,lon_i = getvars('latCell','lonCell') 225 lat_v,lon_v = getvars('latVertex','lonVertex') 226 lat_e,lon_e,angle_e = getvars('latEdge','lonEdge','angleEdge') 227 wee,Riv2 = getvars('weightsOnEdge','kiteAreasOnVertex') 228 230 def get_pdim(self,comm,name): 231 name = self.translate[name] 232 return parallel.PDim(self.nc.dimensions[name], comm) 233 def getvar(self,name): 234 # first deal with special cases 235 if name == 'dual_deg': 236 dual_deg = [3 for i in range(self.dual_num) ] 237 dual_deg = np.ascontiguousarray(dual_deg,dtype=np.int32) 238 # NB : Fortran code expects 32-bit ints 239 return dual_deg 240 # general case 241 name=self.translate[name] 242 return self.nc.variables[name][:] 243 def get_pvar(self,pdim,name): 244 # provides parallel access to a NetCDF variable, with name translation 245 # first deal with special cases 246 if name == 'dual_deg': return parallel.CstPArray1D(pdim, np.int32, 3) 247 if name == 'edge_deg': return parallel.CstPArray1D(pdim, np.int32, 2) 248 # general case 249 name=self.translate[name] 250 return parallel.PArray(pdim, self.nc.variables[name]) 251 def normalize(self, mesh, radius): 252 (trisk_deg, trisk, Ai, 253 de, le, wee, Av, Riv2) = (mesh.trisk_deg, mesh.trisk, mesh.Ai, 254 mesh.de, mesh.le, mesh.wee, mesh.Av, mesh.Riv2) 255 edge_num, dual_num, r2 = de.size, Av.size, radius**2 229 256 # fix normalization of wee and Riv2 weights 230 257 for edge1 in range(edge_num): … … 238 265 r=sum(r) 239 266 if abs(r-1.)>1e-6 : print 'error with Riv2 at vertex ', ivertex 240 241 max_primal_deg, max_dual_deg, max_trisk_deg = getdims('maxEdges','vertexDegree','maxEdges2') 242 243 # CRITICAL : force arrays left, etc. to be contiguous in memory: 244 left,right,up,down = [np.ascontiguousarray(x) for x in (left,right,up,down)] 245 trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)] 267 # scale lengths and areas 268 mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai 246 269 247 print ('Max stencil sizes (div,curl,trisk) : %d, %d, %d' 248 % (max_primal_deg, max_dual_deg, max_trisk_deg) ) 249 250 r2 = radius**2 251 Av = r2*Av 252 fv = f(lon_v,lat_v) # Coriolis parameter 253 self.Ai, self.Av, self.fv = r2*Ai,Av,fv 254 self.le, self.de, self.le_de = radius*le, radius*de, le/de 255 self.trisk_deg, self.trisk, self.wee = trisk_deg, trisk, wee 256 init_mesh(llm, nqdyn, edge_num, primal_num,dual_num, 257 max_trisk_deg, max_primal_deg, max_dual_deg, 258 primal_deg,primal_edge,primal_ne, 259 dual_deg,dual_edge,dual_ne,dual_vertex, 260 left,right,down,up,trisk_deg,trisk, 261 self.Ai,self.Av,self.fv,le/de,Riv2,wee) 262 263 for edge in range(edge_num): 264 iedge = trisk[edge,0:trisk_deg[edge]] 265 if iedge.min()<1 : 266 print 'error' 267 268 self.primal_num, self.edge_num, self.dual_num = primal_num, edge_num, dual_num 269 def period(x) : return (x+2*math.pi)%(2*math.pi) 270 lon_i, lon_v, lon_e = map(period, (lon_i,lon_v,lon_e)) 271 self.lon_i, self.lat_i = lon_i, lat_i 272 self.lon_v, self.lat_v = lon_v, lat_v 273 self.lon_e, self.lat_e, self.angle_e = lon_e, lat_e, angle_e 274 self.primal_deg, self.primal_vertex = primal_deg, primal_vertex 275 self.primal = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi) 276 self.dual_deg, self.dual_vertex = dual_deg, dual_vertex 277 self.dual = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi) 278 self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi) 279 self.dx = de.min() 280 self.lon3D_i, self.ll3D = np.meshgrid(lon_i, range(llm)) 281 self.lat3D_i, self.ll3D = np.meshgrid(lat_i, range(llm)) 270 def compute_ne(num,deg,edges,left,right): 271 ne = np.zeros_like(edges) 272 for cell in range(num): 273 for iedge in range(deg[cell]): 274 edge = edges[cell,iedge]-1 275 if left[edge]==cell+1: ne[cell,iedge]+=1 276 if right[edge]==cell+1: ne[cell,iedge]-=1 277 if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge 278 return ne 279 280 class Abstract_Mesh(Base_class): 282 281 def field_theta(self,n=1): return squeeze((n,self.nqdyn,self.primal_num,self.llm)) 283 282 def field_mass(self,n=1): return squeeze((n,self.primal_num,self.llm)) … … 294 293 ucov[edge,:] = self.de[edge]*(ulon[edge,:]*math.cos(angle)+ulat[edge,:]*math.sin(angle)) 295 294 return ucov 295 def plot(self, tri,data, **kwargs): 296 plt.figure(figsize=(12,4)) 297 plt.gca().set_aspect('equal') 298 plt.tricontourf(tri, data, 20, **kwargs) 299 plt.colorbar() 300 plt.ylim((-90,90)) 301 plt.xlim((0,360)) 296 302 def plot_i(self,data, **kwargs): 297 plot(self.primal,data, **kwargs)303 self.plot(self.primal,data, **kwargs) 298 304 def plot_v(self,data, **kwargs): 299 plot(self.dual,data, **kwargs)305 self.plot(self.dual,data, **kwargs) 300 306 def plot_e(self,data, **kwargs): 301 plot(self.triedge,data, **kwargs) 307 self.plot(self.triedge,data, **kwargs) 308 309 class Unstructured_Mesh(Abstract_Mesh): 310 def __init__(self, gridfile, llm, nqdyn, radius, f): 311 # open mesh file, get main dimensions 312 getvars = gridfile.getvars 313 # get positions, lengths, surfaces and weights 314 le,de,Ai,Av,wee,Riv2 = getvars('le','de','Ai','Av','wee','Riv2') 315 lat_i,lon_i,lat_v,lon_v = getvars('lat_i','lon_i','lat_v','lon_v') 316 lat_e,lon_e,angle_e = getvars('lat_e','lon_e','angle_e') 317 fv = f(lon_v,lat_v) # Coriolis parameter 318 primal_num, dual_num, edge_num = [ x.size for x in Ai, Av, le ] 319 gridfile.primal_num, gridfile.dual_num = primal_num, dual_num 320 321 primal_deg, dual_deg, trisk_deg = getvars('primal_deg','dual_deg','trisk_deg') 322 323 # get indices for stencils 324 # primal <-> edges 325 primal_edge, left_right = getvars('primal_edge','left_right') 326 left,right = left_right[:,0], left_right[:,1] 327 primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right) 328 # dual <-> edges 329 dual_edge, down_up = getvars('dual_edge','down_up') 330 down,up = down_up[:,0], down_up[:,1] 331 dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down) 332 # primal <-> dual, edges <-> edges (TRiSK) 333 primal_vertex, dual_vertex, trisk = getvars('primal_vertex','dual_vertex','trisk') 334 # CRITICAL : force arrays left, etc. to be contiguous in memory: 335 left,right,up,down = [np.ascontiguousarray(x) for x in (left,right,up,down)] 336 trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)] 337 self.set_members(locals(), 'gridfile','llm','nqdyn','radius','f', 338 'primal_num', 'dual_num', 'edge_num', 'primal_deg', 'primal_vertex', 339 'lat_i', 'lon_i', 'lat_v', 'lon_v', 'lat_e', 'lon_e', 'angle_e', 340 'Ai', 'Av', 'fv', 'le', 'de', 341 'trisk_deg', 'trisk', 'wee', 'Riv2') 342 gridfile.normalize(self, radius) 343 self.le_de = self.le/self.de 344 max_primal_deg, max_dual_deg, max_trisk_deg = [ x.shape[1] for x in primal_edge, dual_edge, trisk] 345 unst.init_mesh(llm, nqdyn, edge_num, primal_num,dual_num, 346 max_trisk_deg, max_primal_deg, max_dual_deg, 347 primal_deg,primal_edge,primal_ne, 348 dual_deg,dual_edge,dual_ne,dual_vertex, 349 left,right,down,up,trisk_deg,trisk, 350 self.Ai,self.Av,self.fv,le/de,Riv2,wee) 351 self.primal = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi) 352 self.dual = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi) 353 self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi) 354 355 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes 356 def __init__(self,comm, gridfile): 357 self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars 358 dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, 359 'primal_num','edge_num','dual_num') 360 (primal_deg, primal_vertex, primal_edge, lon_i, 361 lat_i) = get_pvars(dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i' ) 362 (edge_deg, trisk_deg, left_right, down_up, 363 trisk) = get_pvars(dim_edge, 'edge_deg', 'trisk_deg', 'left_right', 'down_up', 'trisk') 364 dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v = get_pvars(dim_dual, 365 'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v') 366 # Let indices start at 0 367 for x in (primal_vertex, primal_edge, dual_vertex,dual_edge, 368 left_right, down_up, trisk ) : x.data = x.data-gridfile.start_index 369 self.edge2cell = Stencil_glob(edge_deg, left_right) 370 self.cell2edge = Stencil_glob(primal_deg, primal_edge) 371 self.edge2vertex = Stencil_glob(edge_deg, down_up) 372 self.vertex2edge = Stencil_glob(dual_deg, dual_edge) 373 self.edge2edge = Stencil_glob(trisk_deg, trisk) 374 print 'Partitioning ...' 375 edge_owner = unst.partition_mesh(trisk_deg, trisk, comm.Get_size()) 376 edge_owner = parallel.LocPArray1D(dim_edge, edge_owner) 377 primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 378 primal_owner = parallel.LocPArray1D(dim_primal, primal_owner) 379 self.set_members(locals(), 'dim_primal', 'dim_dual', 'edge_owner', 380 'primal_owner', 'primal_deg', 'primal_vertex', 381 'lon_v', 'lat_v', 'lon_i', 'lat_i') 382 # self.dim_primal, self.edge_owner, self.primal_owner = dim_primal, edge_owner, primal_owner 383 384 def chain(start, links): 385 for link in links: 386 start = link(start).neigh_set 387 yield start 388 389 def patches(degree, bounds, lon, lat): 390 for i in range(degree.size): 391 nb_edge=degree[i] 392 bounds_cell = bounds[i,0:nb_edge] 393 lat_cell = lat[bounds_cell] 394 lon_cell = lon[bounds_cell] 395 orig=lon_cell[0] 396 lon_cell = lon_cell-orig+180. 397 lon_cell = np.mod(lon_cell,360.) 398 lon_cell = lon_cell+orig-180. 399 # if np.abs(lon_cell-orig).max()>10. : 400 # print '%d patches :'%mpi_rank, lon_cell 401 lonlat_cell = np.zeros((nb_edge,2)) 402 lonlat_cell[:,0],lonlat_cell[:,1] = lon_cell,lat_cell 403 polygon = Polygon(lonlat_cell, True) 404 yield polygon 405 406 class Local_Mesh(Base_class): 407 def __init__(self, pmesh): 408 dim_primal, edge_owner, primal_owner = pmesh.dim_primal, pmesh.edge_owner, pmesh.primal_owner 409 print 'Constructing halos ...' 410 edges_E0 = find_my_cells(edge_owner) 411 cells_C0, edges_E1, vertices_V1, edges_E2, cells_C1 = chain( 412 edges_E0, ( pmesh.edge2cell, pmesh.cell2edge, pmesh.edge2vertex, pmesh.vertex2edge, pmesh.edge2cell) ) 413 edges_E0, edges_E1, edges_E2 = progressive_list(edges_E0, edges_E1, edges_E2) 414 cells_C0, cells_C1 = progressive_list(cells_C0, cells_C1) 415 print 'E2,E1,E0 ; C1,C0 : ', map(len, (edges_E2, edges_E1, edges_E0, cells_C1, cells_C0)) 416 mycells, halo_cells = cells_C0, cells_C1 417 get_mycells, get_halo_cells = dim_primal.getter(mycells), dim_primal.getter(halo_cells) 418 self.com_primal = parallel.Halo_Xchange(42, dim_primal, halo_cells, get_halo_cells(primal_owner)) 419 def plot(self, ax, clim, degree, bounds, lon, lat, data): 420 nb_vertex = lon.size # global 421 p = list(patches(degree, bounds, lon, lat)) 422 p = PatchCollection(p, linewidth=0.01) 423 p.set_array(data) # set values at each polygon (cell) 424 p.set_clim(clim) 425 ax.add_collection(p) 302 426 303 427 #-------------------------------------- Mesh reordering ------------------------------------------# … … 341 465 nn, xx, yy, zz = x.size, toint(x), toint(y), toint(z) 342 466 idx=np.zeros(nn, dtype=np.int64) 343 ker.dynamico_morton_encode(nn, xx,yy,zz, idx)467 unst.ker.dynamico_morton_encode(nn, xx,yy,zz, idx) 344 468 return idx 345 469 -
codes/icosagcm/devel/Python/src/unstructured.pyx
r674 r680 1 1 import time 2 2 import math 3 import dynamico.wrap as wrap 4 from dynamico.libs import libicosa 3 import wrap 4 from libs import libicosa 5 from util import list_stencil 5 6 6 7 from ctypes import c_void_p, c_int, c_double, c_bool … … 289 290 290 291 # Helper functions and interface to ParMETIS 291 # list_stencil converts an adjacency graph from array format index[num_cells, MAX_EDGES] to compressed format292 292 # loc_stencil returns the start/end indices (vtxdist) expected by ParMETIS 293 293 # i.e. index[start:end] with start=vtxdist[cell], end=vtxdist[cell+1] lists the edges of cell 'cell' … … 295 295 import numpy as np 296 296 297 def list_stencil(degree, stencil, cond=lambda x:True):298 for i in range(degree.size):299 for j in range(degree[i]):300 s=stencil[i,j]301 if cond(s): yield stencil[i,j]302 303 297 def loc_stencil(degree, stencil): 304 298 loc=0 -
codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py
r679 r680 4 4 from dynamico import DCMIP 5 5 from dynamico import meshes 6 from dynamico.meshes import MPAS_Mesh as Mesh7 6 import math as math 8 7 import matplotlib.pyplot as plt … … 45 44 nqdyn, preff, Treff = 1, 1e5, 300. 46 45 thermo = dyn.Ideal_perfect(Cpd, Rd, preff, Treff) 47 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 46 meshfile = meshes.MPAS_Format('grids/x1.%d.grid.nc'%grid) 47 mesh = meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 48 48 lev,levp1 = np.arange(llm),np.arange(llm+1) 49 49 lon_i, lat_i, lon_e, lat_e = mesh.lon_i, mesh.lat_i, mesh.lon_e, mesh.lat_e -
codes/icosagcm/devel/Python/test/py/NH_3D_DCMIP31.py
r678 r680 4 4 from dynamico import DCMIP 5 5 from dynamico import meshes 6 from dynamico.meshes import MPAS_Mesh as Mesh7 6 import math as math 8 7 import matplotlib.pyplot as plt … … 49 48 nqdyn, preff, Treff = 1, 1e5, 300. 50 49 thermo = dyn.Ideal_perfect(Cpd, Rd, preff, Treff) 51 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 50 meshfile = meshes.MPAS_Format('grids/x1.%d.grid.nc'%grid) 51 mesh = meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 52 52 lev,levp1 = np.arange(llm),np.arange(llm+1) 53 53 lat_ik,k_i = np.meshgrid(mesh.lat_i,lev, indexing='ij') -
codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py
r666 r680 1 1 print 'Loading DYNAMICO modules ...' 2 2 from dynamico import unstructured as unst 3 from dynamico.meshes import MPAS_ Mesh as Mesh3 from dynamico.meshes import MPAS_Format, Unstructured_Mesh as Mesh 4 4 from dynamico import time_step 5 5 print '...Done' … … 19 19 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 20 20 print 'Reading MPAS mesh ...' 21 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 21 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 22 mesh = Mesh(meshfile, llm, nqdyn, radius, f) 22 23 print '...Done' 23 24 lon, lat = mesh.lon_i, mesh.lat_i -
codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py
r672 r680 1 1 print 'Loading DYNAMICO modules ...' 2 2 from dynamico import unstructured as unst 3 from dynamico.meshes import MPAS_ Mesh asMesh3 from dynamico.meshes import MPAS_Format, Unstructured_Mesh 4 4 from dynamico import time_step 5 5 print '...Done' … … 19 19 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 20 20 print 'Reading MPAS mesh ...' 21 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 21 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 22 mesh=Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 22 23 print '...Done' 23 24 lon, lat = mesh.lon_i, mesh.lat_i -
codes/icosagcm/devel/Python/test/py/partition.py
r672 r680 57 57 # Helper functions to plot unstructured graph 58 58 59 def patches(degree, bounds, lon, lat):60 for i in range(degree.size):61 nb_edge=degree[i]62 bounds_cell = bounds[i,0:nb_edge]63 lat_cell = lat[bounds_cell]64 lon_cell = lon[bounds_cell]65 orig=lon_cell[0]66 lon_cell = lon_cell-orig+180.67 lon_cell = np.mod(lon_cell,360.)68 lon_cell = lon_cell+orig-180.69 # if np.abs(lon_cell-orig).max()>10. :70 # print '%d patches :'%mpi_rank, lon_cell71 lonlat_cell = np.zeros((nb_edge,2))72 lonlat_cell[:,0],lonlat_cell[:,1] = lon_cell,lat_cell73 polygon = Polygon(lonlat_cell, True)74 yield polygon75 76 def plot_mesh(ax, clim, degree, bounds, lon, lat, data):77 nb_vertex = lon.size # global78 p = list(patches(degree, bounds, lon, lat))79 print '%d : plot_mesh %d %d %d'%( mpi_rank, degree.size, len(p), len(data) )80 p = PatchCollection(p, linewidth=0.01)81 p.set_array(data) # set values at each polygon (cell)82 p.set_clim(clim)83 ax.add_collection(p)84 85 59 def local_mesh(get_mycells): 86 mydegree, mybounds = [get_mycells(x) for x in nEdgesOnCell, verticesOnCell] 60 # mydegree, mybounds = [get_mycells(x) for x in nEdgesOnCell, verticesOnCell] 61 mydegree, mybounds = [get_mycells(x) for x in primal_deg, primal_vertex] 87 62 print '%d : len(mydegree)=%d'%(mpi_rank, len(mydegree)) 88 63 vertex_list = sorted(set(unst.list_stencil(mydegree,mybounds))) … … 94 69 return vertex_list, mydegree, mybounds, mylon, mylat 95 70 71 def members(struct, *names): return [struct.__dict__ [name] for name in names] 72 96 73 #--------------- read MPAS grid file ---------------# 97 74 98 #grid = 'x1.2562'99 grid = 'x1.10242'75 grid = 'x1.2562' 76 #grid = 'x1.10242' 100 77 #grid = 'x4.163842' 101 78 print 'Reading MPAS file %s ...'%grid 102 79 103 nc = cdf.Dataset('grids/%s.grid.nc'%grid, "r") 104 dim_cell, dim_edge, dim_vertex = [ 105 parallel.PDim(nc.dimensions[name], comm) 106 for name in 'nCells','nEdges','nVertices'] 107 edge_degree = parallel.CstPArray1D(dim_edge, np.int32, 2) 108 vertex_degree = parallel.CstPArray1D(dim_vertex, np.int32, 3) 109 nEdgesOnCell, verticesOnCell, edgesOnCell, cellsOnCell, latCell = [ 110 parallel.PArray(dim_cell, nc.variables[var]) 111 for var in 'nEdgesOnCell', 'verticesOnCell', 'edgesOnCell', 'cellsOnCell', 'latCell' ] 112 cellsOnVertex, edgesOnVertex, kiteAreasOnVertex, lonVertex, latVertex = [ 113 parallel.PArray(dim_vertex, nc.variables[var]) 114 for var in 'cellsOnVertex', 'edgesOnVertex', 'kiteAreasOnVertex', 'lonVertex', 'latVertex'] 115 nEdgesOnEdge, cellsOnEdge, edgesOnEdge, verticesOnEdge, weightsOnEdge = [ 116 parallel.PArray(dim_edge, nc.variables[var]) 117 for var in 'nEdgesOnEdge', 'cellsOnEdge', 'edgesOnEdge', 'verticesOnEdge', 'weightsOnEdge'] 80 meshfile = meshes.MPAS_Format('grids/%s.grid.nc'%grid) 81 pmesh = meshes.Unstructured_PMesh(comm, meshfile) 82 lmesh = meshes.Local_Mesh(pmesh) 118 83 119 # Indices start at 0 on the C/Python side and at 1 on the Fortran/MPAS side 120 # hence an offset of 1 is added/substracted where needed. 121 for x in (verticesOnCell, edgesOnCell, cellsOnCell, cellsOnVertex, edgesOnVertex, 122 cellsOnEdge, edgesOnEdge, verticesOnEdge) : x.data = x.data-1 123 edge2cell, cell2edge, edge2vertex, vertex2edge, cell2cell, edge2edge = [ 124 meshes.Stencil_glob(a,b) for a,b in 125 (edge_degree, cellsOnEdge), (nEdgesOnCell, edgesOnCell), 126 (edge_degree, verticesOnEdge), (vertex_degree, edgesOnVertex), 127 (nEdgesOnCell, cellsOnCell), (nEdgesOnEdge, edgesOnEdge) ] 84 (primal_deg, primal_vertex, dim_vertex, dim_cell, cell_owner, 85 lonVertex, latVertex, lonCell, latCell) = members( 86 pmesh, 'primal_deg', 'primal_vertex', 'dim_dual', 'dim_primal', 'primal_owner', 87 'lon_v', 'lat_v', 'lon_i', 'lat_i') 128 88 129 #---------------- partition edges and cells ------------------# 130 131 print 'Partitioning ...' 132 133 edge_owner = unst.partition_mesh(nEdgesOnEdge, edgesOnEdge, mpi_size) 134 edge_owner = parallel.LocPArray1D(dim_edge, edge_owner) 135 cell_owner = meshes.partition_from_stencil(edge_owner, nEdgesOnCell, edgesOnCell) 136 cell_owner = parallel.LocPArray1D(dim_cell, cell_owner) 137 138 #--------------------- construct halos -----------------------# 139 140 print 'Constructing halos ...' 141 142 def chain(start, links): 143 for link in links: 144 start = link(start).neigh_set 145 yield start 146 147 edges_E0 = meshes.find_my_cells(edge_owner) 148 cells_C0, edges_E1, vertices_V1, edges_E2, cells_C1 = chain( 149 edges_E0, ( edge2cell, cell2edge, edge2vertex, vertex2edge, edge2cell) ) 150 151 edges_E0, edges_E1, edges_E2 = meshes.progressive_list(edges_E0, edges_E1, edges_E2) 152 cells_C0, cells_C1 = meshes.progressive_list(cells_C0, cells_C1) 153 154 print 'E2,E1,E0 ; C1,C0 : ', map(len, (edges_E2, edges_E1, edges_E0, cells_C1, cells_C0)) 155 156 #com_edges = parallel.Halo_Xchange(24, dim_edge, edges_E2, dim_edge.get(edges_E2, edge_owner)) 157 158 mycells, halo_cells = cells_C0, cells_C1 159 get_mycells, get_halo_cells = dim_cell.getter(mycells), dim_cell.getter(halo_cells) 160 com_cells = parallel.Halo_Xchange(42, dim_cell, halo_cells, get_halo_cells(cell_owner)) 161 162 local_num, total_num = np.zeros(1), np.zeros(1) 89 local_num, total_num, com_cells = np.zeros(1), np.zeros(1), lmesh.com_primal 163 90 local_num[0]=com_cells.own_len 164 91 comm.Reduce(local_num, total_num, op=MPI.SUM, root=0) … … 167 94 #---------------------------- plot -----------------------------# 168 95 169 if True: 170 print 'Plotting ...' 96 print 'Plotting ...' 171 97 172 173 98 halo_vertex_list, mydegree, mybounds, mylon, mylat = local_mesh(com_cells.get_all) 99 buf = parallel.LocalArray1(com_cells) 174 100 175 176 177 178 179 plot_mesh(ax,[-math.pi/2,math.pi/2], mydegree, mybounds, mylon, mylat, buf.data)180 181 182 101 fig, ax = plt.subplots() 102 buf.read_own(latCell) # reads only own values 103 buf.data = np.cos(10.*buf.data) 104 buf.update() # updates halo 105 lmesh.plot(ax,[-math.pi/2,math.pi/2], mydegree, mybounds, mylon, mylat, buf.data) 106 plt.xlim(-190.,190.) 107 plt.ylim(-90.,90.) 108 plt.savefig('fig_partition/A%03d.pdf'%mpi_rank, dpi=1600) 183 109 184 185 186 187 plot_mesh(ax,[0,mpi_rank+1], mydegree, mybounds, mylon, mylat, buf.data)188 189 190 110 fig, ax = plt.subplots() 111 buf.read_own(cell_owner) 112 buf.update() 113 lmesh.plot(ax,[0,mpi_rank+1], mydegree, mybounds, mylon, mylat, buf.data) 114 plt.xlim(-190.,190.) 115 plt.ylim(-90.,90.) 116 plt.savefig('fig_partition/B%03d.pdf'%mpi_rank, dpi=1600) -
codes/icosagcm/devel/Python/test/py/test_xios.py
r663 r680 1 1 #print 'import dynamico' 2 2 #import dynamico 3 print 'import dynamico.unstructured'4 3 import dynamico.unstructured as unst 5 print 'import dynamico.xios'6 4 import dynamico.xios as xios 7 print 'Done.' 8 from dynamico.meshes import MPAS_Mesh as Mesh, radian 5 from dynamico.meshes import MPAS_Format, Unstructured_Mesh, radian 9 6 10 7 from mpi4py import MPI … … 38 35 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 39 36 print 'Reading MPAS mesh ...' 40 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 37 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 38 mesh = Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 41 39 print '...Done' 42 40 43 41 lon_i, lat_i, lon_v, lat_v = [x*radian for x in mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v] 44 #lon_i, lat_i, lon_v, lat_v = mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v45 42 46 43 #--------------------------------- initialize XIOS --------------------------
Note: See TracChangeset
for help on using the changeset viewer.