Changeset 631 for codes/icosagcm/devel
- Timestamp:
- 12/12/17 16:04:23 (7 years ago)
- Location:
- codes/icosagcm/devel
- Files:
-
- 1 added
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/src/unstructured.pyx
r622 r631 2 2 import math 3 3 import numpy as np 4 import netCDF4 as cdf5 import matplotlib.tri as tri6 import matplotlib.pyplot as plt7 8 4 import dynamico.wrap as wrap 9 5 from ctypes import c_void_p, c_int, c_double, c_bool 10 radian=180/math.pi11 6 12 7 #------------- direct Cython interface to DYNAMICO routines -------------# … … 39 34 ['dynamico_xios_set_timestep',c_double], 40 35 ['dynamico_xios_update_calendar',c_int], 41 # ['init_params',c_double],42 36 ['dynamico_init_mesh',c_void_p,13], 43 37 ['dynamico_init_metric', c_void_p,6], … … 146 140 (dm,dS,du_slow,dPhi_slow,dW_slow)) 147 141 148 #------------------------ -------- Hybrid mass-based coordinate-------------142 #------------------------ Copy mesh info to Fortran side ------------------- 149 143 150 # compute hybrid coefs from distribution of mass 151 def compute_hybrid_coefs(mass): 152 nx,llm=mass.shape 153 mass_dak = np.zeros((nx,llm)) 154 mass_dbk = np.zeros((nx,llm)) 155 mass_bl = np.zeros((nx,llm+1))+1. 156 for i in range(nx): 157 m_i = mass[i,:] 158 dbk_i = m_i/sum(m_i) 159 mass_dbk[i,:] = dbk_i 160 mass_bl[i,1:]= 1.-dbk_i.cumsum() 161 return mass_bl, mass_dak, mass_dbk 144 def init_mesh(llm, nqdyn, edge_num, primal_num, dual_num, 145 max_trisk_deg, max_primal_deg, max_dual_deg, 146 primal_nb, primal_edge, primal_ne, 147 dual_nb,dual_edge,dual_ne,dual_vertex, 148 left,right,down,up,trisk_deg,trisk, 149 Ai, Av, fv, le_de, Riv2, wee): 150 setvars( ('llm','nqdyn','edge_num','primal_num','dual_num', 151 'max_trisk_deg','max_primal_deg','max_dual_deg'), 152 (llm, nqdyn, edge_num, primal_num, dual_num, 153 max_trisk_deg, max_primal_deg, max_dual_deg) ) 154 print('init_mesh ...') 155 ker.dynamico_init_mesh(primal_nb,primal_edge,primal_ne, 156 dual_nb,dual_edge,dual_ne,dual_vertex, 157 left,right,down,up,trisk_deg,trisk) 158 print ('...done') 159 print('init_metric ...') 160 ker.dynamico_init_metric(Ai,Av,fv,le_de,Riv2,wee) 161 print ('...done') 162 162 163 #----------------------- Cartesian mesh ----------------------- 164 165 def squeeze(dims): 166 # return np.zeros(dims) 167 return np.zeros([n for n in dims if n>1]) 168 169 # arrays is a list of arrays 170 # vals is a list of tuples 171 # each tuple is stored in each array 172 def put(ij, deg, arrays, vals): 173 k=0 174 for vv in vals: # vv is a tuple of values to be stored in arrays 175 for array,v in zip(arrays,vv): 176 array[ij,k]=v 177 k=k+1 178 deg[ij]=k 179 180 class Cartesian_mesh: 181 def __init__(self,nx,ny,llm,nqdyn,Lx,Ly,f): 182 dx,dy = Lx/nx, Ly/ny 183 self.dx, self.dy, self.f = dx,dy,f 184 self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn 185 self.field_z = self.field_mass 186 # 1D coordinate arrays 187 x=(np.arange(nx)-nx/2.)*Lx/nx 188 y=(np.arange(ny)-ny/2.)*Ly/ny 189 lev=np.arange(llm) 190 levp1=np.arange(llm+1) 191 self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1 192 # 3D coordinate arrays 193 self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij') 194 self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij') 195 # beware conventions for indexing 196 # Fortran order : llm,nx*ny,nqdyn / indices start at 1 197 # Python order : nqdyn,ny,nx,llm / indices start at 0 198 # indices below follow Fortran while x,y follow Python/C 199 index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1 200 indexu=lambda x,y : 2*index(x,y)-1 201 indexv=lambda x,y : 2*index(x,y) 202 indices = lambda shape : np.zeros(shape,np.int32) 203 204 primal_nb = indices(nx*ny) 205 primal_edge = indices((nx*ny,4)) 206 primal_ne = indices((nx*ny,4)) 207 208 dual_nb = indices(nx*ny) 209 dual_edge = indices((nx*ny,4)) 210 dual_ne = indices((nx*ny,4)) 211 dual_vertex = indices((nx*ny,4)) 212 Riv2 = np.zeros((nx*ny,4)) 213 214 left = indices(2*nx*ny) 215 right = indices(2*nx*ny) 216 up = indices(2*nx*ny) 217 down = indices(2*nx*ny) 218 le_de = np.zeros(2*nx*ny) 219 220 trisk_deg = indices(2*nx*ny) 221 trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge 222 wee = np.zeros((2*nx*ny,4)) 223 224 for x in range(nx): 225 for y in range(ny): 226 # NB : Fortran indices start at 1 227 # primal cells 228 put(index(x,y)-1,primal_nb,(primal_edge,primal_ne), 229 ((indexu(x,y),1), 230 (indexv(x,y),1), 231 (indexu(x-1,y),-1), 232 (indexv(x,y-1),-1) )) 233 # dual cells 234 put(index(x,y)-1,dual_nb,(dual_edge,dual_vertex,dual_ne,Riv2), 235 ((indexv(x+1,y),index(x,y),1,.25), 236 (indexu(x,y+1),index(x+1,y),-1,.25), 237 (indexv(x,y),index(x+1,y+1),-1,.25), 238 (indexu(x,y),index(x,y+1),1,.25) )) 239 # edges : 240 # left and right are adjacent primal cells 241 # flux is positive when going from left to right 242 # up and down are adjacent dual cells (vertices) 243 # circulation is positive when going from down to up 244 # u-points 245 ij =indexu(x,y)-1 246 left[ij]=index(x,y) 247 right[ij]=index(x+1,y) 248 down[ij]=index(x,y-1) 249 up[ij]=index(x,y) 250 le_de[ij]=dy/dx 251 put(ij,trisk_deg,(trisk,wee),( 252 (indexv(x,y),-.25), 253 (indexv(x+1,y),-.25), 254 (indexv(x,y-1),-.25), 255 (indexv(x+1,y-1),-.25))) 256 # v-points 257 ij = indexv(x,y)-1 258 left[ij]=index(x,y) 259 right[ij]=index(x,y+1) 260 down[ij]=index(x,y) 261 up[ij]=index(x-1,y) 262 le_de[ij]=dx/dy 263 put(ij,trisk_deg,(trisk,wee),( 264 (indexu(x,y),.25), 265 (indexu(x-1,y),.25), 266 (indexu(x,y+1),.25), 267 (indexu(x-1,y+1),.25))) 268 setvars(('llm','nqdyn','edge_num','primal_num','dual_num', 269 'max_trisk_deg','max_primal_deg','max_dual_deg'), 270 (llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4) ) 271 print('init_mesh ...') 272 ker.dynamico_init_mesh(primal_nb,primal_edge,primal_ne, 273 dual_nb,dual_edge,dual_ne,dual_vertex, 274 left,right,down,up,trisk_deg,trisk) 275 print ('...done') 276 277 Aiv=np.zeros(nx*ny)+dx*dy # Ai=Av=dx*dy 278 279 print('init_metric ...') 280 ker.dynamico_init_metric(Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee) 281 print ('...done') 282 def field_theta(self): return squeeze((self.nqdyn,self.ny,self.nx,self.llm)) 283 def field_mass(self): return squeeze((self.ny,self.nx,self.llm)) 284 def field_z(self): return squeeze((self.ny,self.nx,self.llm)) 285 def field_w(self): return squeeze((self.ny,self.nx,self.llm+1)) 286 def field_u(self): return np.zeros((self.ny,2*self.nx,self.llm)) 287 def field_ps(self): return squeeze((self.ny,self.nx)) 288 def ucomp(self,u): return u[:,range(0,2*self.nx,2),:] 289 def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u 290 def vcomp(self,u): return u[:,range(1,2*self.nx,2),:] 291 292 #---------------------- MPAS fully unstructured mesh ------------------------ 293 294 def compute_ne(num,deg,edges,left,right): 295 ne = np.zeros_like(edges) 296 for cell in range(num): 297 for iedge in range(deg[cell]): 298 edge = edges[cell,iedge]-1 299 if left[edge]==cell+1: ne[cell,iedge]+=1 300 if right[edge]==cell+1: ne[cell,iedge]-=1 301 if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge 302 return ne 303 304 def plot(tri,data): 305 plt.figure(figsize=(12,4)) 306 plt.gca().set_aspect('equal') 307 plt.tricontourf(tri, data, 20) 308 plt.colorbar() 309 plt.ylim((-90,90)) 310 plt.xlim((0,360)) 311 312 class MPAS_Mesh: 313 def __init__(self, gridfile, llm, nqdyn, radius, f): 314 self.gridfile, self.llm, self.nqdyn = gridfile,llm,nqdyn 315 self.radius, self.f = radius, f 316 # open mesh file, get main dimensions 317 try: 318 nc = cdf.Dataset(gridfile, "r") 319 except RuntimeError: 320 print """ 321 Unable to open grid file %s, maybe you forgot to download it ? 322 To do so, go to the 'Python/' dir and execute './get_MPAS_grids.sh'. 323 """ % gridfile 324 raise 325 326 def getdims(*names): return [len(nc.dimensions[name]) for name in names] 327 def getvars(*names): 328 for name in names : print "getvar %s ..."%name 329 time1=time.time() 330 ret=[nc.variables[name][:] for name in names] 331 print "... Done in %f seconds"%(time.time()-time1) 332 return ret 333 primal_num, edge_num, dual_num = getdims('nCells','nEdges','nVertices') 334 print 'Number of primal cells, dual cells and edges : %d, %d, %d '%(primal_num,dual_num,edge_num) 335 primal_deg, trisk_deg = getvars('nEdgesOnCell','nEdgesOnEdge') 336 dual_deg = [3 for i in range(dual_num) ] 337 dual_deg = np.ascontiguousarray(dual_deg,dtype=np.int32) # NB : Fortran code expects 32-bit ints 338 339 # get indices for stencils 340 # primal -> vertices (unused) 341 primal_vertex, dual_vertex = getvars('verticesOnCell','cellsOnVertex') 342 # primal <-> edges 343 primal_edge, left_right = getvars('edgesOnCell','cellsOnEdge') 344 left,right = left_right[:,0], left_right[:,1] 345 primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right) 346 # dual <-> edges 347 dual_edge, down_up = getvars('edgesOnVertex','verticesOnEdge') 348 down,up = down_up[:,0], down_up[:,1] 349 dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down) 350 # primal <-> dual, edges <-> edges (TRiSK) 351 dual_vertex, trisk = getvars('cellsOnVertex','edgesOnEdge') 352 # get positions, lengths, surfaces and weights 353 le,de,Ai,Av = getvars('dvEdge','dcEdge','areaCell','areaTriangle') 354 lat_i,lon_i = getvars('latCell','lonCell') 355 lat_v,lon_v = getvars('latVertex','lonVertex') 356 lat_e,lon_e,angle_e = getvars('latEdge','lonEdge','angleEdge') 357 wee,Riv2 = getvars('weightsOnEdge','kiteAreasOnVertex') 358 359 # fix normalization of wee and Riv2 weights 360 for edge1 in range(edge_num): 361 for i in range(trisk_deg[edge1]): 362 edge2=trisk[edge1,i]-1 # NB Fortran vs C/Python indexing 363 wee[edge1,i] = de[edge1]*wee[edge1,i]/le[edge2] 364 for ivertex in range(dual_num): 365 for j in range(3): 366 Riv2[ivertex,j]=Riv2[ivertex,j]/Av[ivertex] 367 r=Riv2[ivertex,:] 368 r=sum(r) 369 if abs(r-1.)>1e-6 : print 'error with Riv2 at vertex ', ivertex 370 371 max_primal_deg, max_dual_deg, max_trisk_deg = getdims('maxEdges','vertexDegree','maxEdges2') 372 373 # CRITICAL : force arrays left, etc. to be contiguous in memory: 374 left,right,up,down = [np.ascontiguousarray(x) for x in (left,right,up,down)] 375 trisk,wee,primal_edge = [np.ascontiguousarray(x) for x in (trisk,wee,primal_edge)] 376 377 print ('Max stencil sizes (div,curl,trisk) : %d, %d, %d' 378 % (max_primal_deg, max_dual_deg, max_trisk_deg) ) 379 380 r2 = radius**2 381 Av = r2*Av 382 fv = f(lon_v,lat_v) # Coriolis parameter 383 self.Ai, self.Av, self.fv = r2*Ai,Av,fv 384 self.le, self.de, self.le_de = radius*le, radius*de, le/de 385 self.trisk_deg, self.trisk, self.wee = trisk_deg, trisk, wee 386 setvars(('llm','nqdyn','edge_num','primal_num','dual_num', 387 'max_trisk_deg','max_primal_deg','max_dual_deg'), 388 (llm, nqdyn, edge_num, primal_num,dual_num, 389 max_trisk_deg, max_primal_deg, max_dual_deg) ) 390 391 ker.dynamico_init_mesh(primal_deg,primal_edge,primal_ne, 392 dual_deg,dual_edge,dual_ne,dual_vertex, 393 left,right,down,up,trisk_deg,trisk) 394 ker.dynamico_init_metric(self.Ai,self.Av,self.fv,le/de,Riv2,wee) 395 396 for edge in range(edge_num): 397 iedge = trisk[edge,0:trisk_deg[edge]] 398 if iedge.min()<1 : 399 print 'error' 400 401 self.primal_num, self.edge_num, self.dual_num = primal_num, edge_num, dual_num 402 def period(x) : return (x+2*math.pi)%(2*math.pi) 403 lon_i, lon_v, lon_e = map(period, (lon_i,lon_v,lon_e)) 404 self.lon_i, self.lat_i = lon_i, lat_i 405 self.lon_v, self.lat_v = lon_v, lat_v 406 self.lon_e, self.lat_e, self.angle_e = lon_e, lat_e, angle_e 407 self.primal_deg, self.primal_vertex = primal_deg, primal_vertex 408 self.primal = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi) 409 self.dual_deg, self.dual_vertex = dual_deg, dual_vertex 410 self.dual = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi) 411 self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi) 412 self.dx = de.min() 413 self.lon3D_i, self.ll3D = np.meshgrid(lon_i, range(llm)) 414 self.lat3D_i, self.ll3D = np.meshgrid(lat_i, range(llm)) 415 def field_theta(self): return squeeze((self.nqdyn,self.primal_num,self.llm)) 416 def field_mass(self): return squeeze((self.primal_num,self.llm)) 417 def field_z(self): return squeeze((self.dual_num,self.llm)) 418 def field_w(self): return squeeze((self.primal_num,self.llm+1)) 419 def field_u(self): return squeeze((self.edge_num,self.llm)) 420 def field_ps(self): return squeeze((self.primal_num,)) 421 def ucov2D(self, ulon, ulat): 422 return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e)) 423 def ucov3D(self, ulon, ulat): 424 ucov = np.zeros((self.edge_num,ulon.shape[1])) 425 for edge in range(self.edge_num): 426 angle=self.angle_e[edge] 427 ucov[edge,:] = self.de[edge]*(ulon[edge,:]*math.cos(angle)+ulat[edge,:]*math.sin(angle)) 428 return ucov 429 def plot_i(self,data): 430 plot(self.primal,data) 431 def plot_v(self,data): 432 plot(self.dual,data) 433 def plot_e(self,data): 434 plot(self.triedge,data) 435 436 #-------------------------------------- Mesh partitioning ------------------------------------------# 163 #------------------------ Mesh partitioning ------------------------ 437 164 438 165 # Helper functions and interface to ParMETIS … … 464 191 ker.dynamico_partition_graph(mpi_rank, mpi_size, vtxdist, xadj_loc, adjncy_loc, nparts, owner) 465 192 return owner 466 467 def partition_from_stencil(owner2, degree, stencil):468 # given a stencil dim1->dim2 and owner2 on dim2, define owner[i] on dim1 as min(stencil[i,:] if i is even, max if odd469 dim1, dim2= degree.dim, owner2.dim470 degree, stencil, n = degree.data, stencil.data, dim1.end-dim1.start471 cells2 = list_stencil(degree, stencil)472 cells2 = sorted(list(set(list(cells2)))) # list of cells for which we need to know owner2473 get2 = dim2.getter(cells2)474 owner1, owner2, glob2loc = np.zeros(n, dtype=np.int32), get2(owner2), get2.dict475 for i in range(n):476 owners = [ owner2[glob2loc[stencil[i,j]]] for j in range(degree[i]) ]477 if i%2 == 0 : owner1[i] = min(owners)478 else : owner1[i] = max(owners)479 return owner1480 481 def find_my_cells(owner): # a PArray1D containing the data returned by partition_mesh482 dim, comm, owner = owner.dim, owner.dim.comm, owner.data483 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()484 cells=[set() for i in range(mpi_size)]485 for i in range(len(owner)):486 cells[owner[i]].add(dim.start+i)487 cells = [sorted(list(cells[i])) for i in range(mpi_size)]488 mycells = comm.alltoall(cells)489 mycells = sorted(sum(mycells, [])) # concatenate into a single list490 return mycells491 492 #---------------------------------- Stencil management ---------------------------------------#493 494 # Class Stencil represents an adjacency relationship (e.g. cell=>edges)495 # using adjacency information read from PArrays496 # It computes a list of "edges" adjacent to a given list of "cells"497 # This is used to form the sets E0 -> C0 -> E1 -> V1 -> E2 -> C1498 # which are then used to form lists of global indices for V1,C1,E2 such that499 # C0, E0, E1 form contiguous subsets of C1, E2 starting from 0500 501 def reindex(vertex_dict, degree, bounds):502 for i in range(degree.size):503 for j in range(degree[i]):504 bounds[i,j] = vertex_dict[bounds[i,j]]505 506 class Stencil_glob:507 def __init__(self, degree, neigh, weight=None):508 self.degree, self.neigh, self.weight = degree, neigh, weight509 def __call__(self, cells, neigh_dict=None):510 return Stencil(cells, self.degree, self.neigh, neigh_dict, self.weight)511 512 class Stencil:513 def __init__(self, cells, degree, neigh, neigh_dict, weight=None):514 get = degree.dim.getter(cells)515 mydegree, myneigh = [get(x) for x in (degree, neigh)]516 if not weight is None : myweight = get(weight)517 if neigh_dict is None :518 keep = lambda n : True519 else : # if neigh_dict is present, only neighbours in dict are retained520 keep = lambda n : n in neigh_dict521 neigh_set = list_stencil(mydegree, myneigh, keep)522 self.neigh_set = list(set(list(neigh_set) ))523 rej=0524 for i in range(len(mydegree)): # keep only elements in neigh_dict, in-place525 k=0526 for j in range(mydegree[i]):527 n=myneigh[i,j]528 if keep(n):529 myneigh[i,k]=n530 if not weight is None : myweight[i,k]=myweight[i,j]531 k=k+1532 else:533 rej=rej+1534 mydegree[i]=k535 if neigh_dict is None:536 neigh_dict = {j:i for i,j in enumerate(self.neigh_set)}537 myneigh_loc = reindex(neigh_dict, mydegree, myneigh)538 self.degree, self.neigh_glob, self.neigh_loc = mydegree, myneigh, myneigh_loc539 540 def progressive_iter(mylist, cell_lists):541 for thelist in cell_lists:542 mylist = mylist + list(set(thelist)-set(mylist))543 yield mylist544 def progressive_list(*cell_lists):545 # cell_lists : a tuple of lists of global indices, with each list a subset of the next546 # returns : a list 'mylist' such that for each list 'thelist' in cell_lists, thelist = mylist[0:len(thelist)]547 # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges548 return list(progressive_iter([], cell_lists))549 -
codes/icosagcm/devel/Python/test/py/RSW_2D.py
r617 r631 1 from dynamico.meshes import Cartesian_mesh as Mesh 1 2 from dynamico import unstructured as unst 2 3 from dynamico import dyn … … 15 16 16 17 unst.setvar('g',g) 17 mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,f)18 mesh = Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) 18 19 caldyn = unst.Caldyn_RSW(mesh) 19 20 -
codes/icosagcm/devel/Python/test/py/RSW_MPAS_W02.py
r617 r631 1 import sys2 1 print 'Loading modules ...' 3 sys.stdout.flush()4 5 2 import math as math 6 # select non-interactive backend, cf http://stackoverflow.com/questions/4931376/generating-matplotlib-graphs-without-a-running-x-server7 import matplotlib8 matplotlib.use('Agg')9 3 import matplotlib.pyplot as plt 10 4 import numpy as np 5 print '...Done' 11 6 12 7 print 'Loading DYNAMICO modules ...' 13 sys.stdout.flush()14 8 from dynamico import unstructured as unst 9 from dynamico.meshes import MPAS_Mesh as Mesh 15 10 from dynamico import time_step 16 11 print '...Done' 17 sys.stdout.flush()18 12 19 13 grid, llm, nqdyn = 10242, 1,1 # 2562, 10242, 40962 … … 22 16 23 17 print 'Omega, planetary PV', Omega, 2*Omega/gh0 24 sys.stdout.flush()25 18 26 19 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 27 20 print 'Reading MPAS mesh ...' 28 sys.stdout.flush() 29 mesh = unst.MPAS_Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 21 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 30 22 print '...Done' 31 sys.stdout.flush()32 23 lon, lat = mesh.lon_i, mesh.lat_i 33 24 x,y,z = np.cos(lat)*np.cos(lon), np.cos(lat)*np.sin(lon), np.sin(lat) -
codes/icosagcm/devel/Python/test/py/bubble.py
r617 r631 4 4 import matplotlib.animation as manimation 5 5 6 from dynamico.meshes import Cartesian_mesh as Mesh 6 7 import dynamico.dyn as dyn 7 8 import dynamico.time_step as time_step … … 94 95 95 96 unst.setvar('nb_threads', 1) 96 mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,ny*dx,0.)97 mesh = Mesh(nx,ny,llm,nqdyn,Lx,ny*dx,0.) 97 98 xx_ik, xx_il, ll = mesh.xx[:,0,:]/1000, mesh.xxp1[:,0,:]/1000, mesh.ll[:,0,:] 98 99 -
codes/icosagcm/devel/Python/test/py/slice_GW_NH.py
r617 r631 1 from dynamico.meshes import Cartesian_mesh as Mesh 1 2 from dynamico import unstructured as unst 2 3 from dynamico import dyn … … 91 92 Lx, nx, ztop, llm = 2e5, 400, 3e4, 60 92 93 nqdyn, ny, dx = 1, 1, Lx/nx 93 mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)94 mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.) 94 95 xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:] 95 96 -
codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py
r617 r631 1 from dynamico.meshes import Cartesian_mesh as Mesh 1 2 from dynamico import unstructured as unst 2 3 from dynamico import dyn … … 58 59 Lx, nx, llm = 3e5, 300, 20 59 60 nqdyn, ny, dx = 1, 1, Lx/nx 60 mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)61 mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.) 61 62 xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:] 62 63 metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP31(Lx, nx, llm) -
codes/icosagcm/devel/Python/test/py/test_xios.py
r622 r631 4 4 print '%d/%d starting'%(mpi_rank,mpi_size) 5 5 6 import sys7 6 import numpy as np 8 #import dynamico.wrap as wrap9 7 print 'import dynamico.unstructured' 10 8 import dynamico.unstructured as unst … … 12 10 import dynamico.xios as xios 13 11 print 'Done.' 12 from dynamico.meshes import MPAS_Mesh as Mesh, radian 14 13 15 14 def boundaries(degree,points,lon,lat): … … 29 28 #----------------------------- read MPAS mesh -------------------------------- 30 29 31 radian=unst.radian32 30 grid, llm, nqdyn, nqtot = 10242, 1,1,1 # 2562, 10242, 40962 33 31 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 34 32 N, T, courant = 40, 10800., 1.2 # simulation length = N*T 35 33 print 'Omega, planetary PV', Omega, 2*Omega/gh0 36 sys.stdout.flush()37 34 38 35 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 39 36 print 'Reading MPAS mesh ...' 40 sys.stdout.flush() 41 mesh = unst.MPAS_Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 37 mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) 42 38 print '...Done' 43 sys.stdout.flush()44 39 45 40 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] -
codes/icosagcm/devel/make_python
r626 r631 94 94 cd $DYNAMICO_ROOT 95 95 96 for module in DCMIP time_step dyn xios unstructured; do96 for module in xios meshes dyn time_step DCMIP; do 97 97 echo "from dynamico import $module" 98 98 python -c "from dynamico import $module"
Note: See TracChangeset
for help on using the changeset viewer.