- Timestamp:
- 10/15/18 01:04:44 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r759 r760 167 167 168 168 primal_i, primal_j = [ x.astype(np.int32) for x in lon_i, lat_i] 169 edge_i, edge_j = [ x.astype(np.int32) for x in lon_e, lat_e] 169 170 lon_i, lon_v, lon_e = [x*dx-Lx/2 for x in lon_i, lon_v, lon_e] 170 171 lat_i, lat_v, lat_e = [y*dy-Ly/2 for y in lat_i, lat_v, lat_e] … … 183 184 'Aiv','lon_i','lat_i','lon_v','lat_v', 184 185 'le_de','le','de','lon_e','lat_e','angle_e', 185 'primal_i','primal_j' )186 'primal_i','primal_j','edge_i','edge_j') 186 187 def ncwrite(self, name): 187 188 """The method writes Cartesian mesh on the disc. … … 249 250 ("lon_e","f8", self.lon_e), 250 251 ("lat_e","f8", self.lat_e), 251 ("angle_e","f8", self.angle_e)] ) 252 ("angle_e","f8", self.angle_e), 253 ("edge_i","i4",self.edge_i), 254 ("edge_j","i4",self.edge_j)] ) 252 255 253 256 create_vars( ("edge","TWO"), … … 473 476 if gridfile.meshtype == 'curvilinear': 474 477 self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 478 self.edge_i, self.edge_j = get_pvars(dim_edge, 'edge_i','edge_j') 475 479 476 480 # Let indices start at 0 … … 484 488 self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2) 485 489 self.edge2edge = Stencil_glob(trisk_deg, trisk, wee) 486 print 'Partitioning ...'487 edge_owner = unst.partition_mesh(trisk_deg, trisk, comm.Get_size())488 edge_owner = parallel.LocPArray1D(dim_edge, edge_owner)489 primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge)490 primal_owner = parallel.LocPArray1D(dim_primal, primal_owner)491 490 self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge', 492 ' edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex',491 'primal_deg', 'primal_vertex', 'primal_edge', 'trisk_deg', 'trisk', 493 492 'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 494 493 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 495 # if meshtype == 'curvilinear' : # read extra information from mesh file 496 # self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 497 494 def partition(self, edge_owner): 495 self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner) 496 primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge) 497 self.primal_owner = parallel.LocPArray1D(self.dim_primal, primal_owner) 498 499 def partition_metis(self): 500 print 'Partitioning with ParMetis...' 501 self.partition(unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())) 502 503 def partition_curvilinear(self, ni,nj): 504 def owner(dim,i,j): 505 owner_i, owner_j = (ni*i.data)/nx, (nj*j.data)/ny 506 return parallel.LocPArray1D(dim, owner_i + ni*owner_j) 507 nx, ny = self.gridfile.nx, self.gridfile.ny 508 print 'Partitioning %d x %d cells in %d x %d blocks ...'%(nx,ny,ni,nj) 509 n = self.comm.Get_size() 510 assert n == ni*nj, 'Mismatch between ni x nj = %d and MPI_SIZE = %d.'%(ni*nj, n) 511 self.edge_owner = owner(self.dim_edge, self.edge_i, self.edge_j) 512 self.primal_owner = owner(self.dim_primal, self.primal_i, self.primal_j) 513 print 'partition_curvilinear %d :'%(self.comm.Get_rank()), self.primal_owner.data 514 498 515 def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ] 499 516 -
codes/icosagcm/devel/Python/dynamico/xios.py
r759 r760 3 3 from unstructured import ker 4 4 from dynamico.meshes import radian 5 from cxios import finalize5 from mpi4py import MPI 6 6 7 7 #----------------------------------------------------------------------------- … … 9 9 def setup_curvilinear(cat,id,nx,ny,displ,index_own_i,index_own_j,lon,lat): 10 10 mesh = cxios.Handle(cat,id) 11 print 'setup_curvilinear : index_own_i.shape', index_own_i.shape 12 print 'setup_curvilinear : index_own_j.shape', index_own_j.shape 11 def log(name,data) : print 'setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max() 12 log('index_own_i',index_own_i) 13 log('index_own_j',index_own_j) 14 log('lon',lon) 15 log('lat',lat) 13 16 mesh.set_attr(type='curvilinear', 14 17 ni_glo=nx, i_index=index_own_i, … … 44 47 #----------------------------------------------------------------------------- 45 48 46 class XIOS_Context: 49 class Client: 50 def __enter__(self): 51 ker.dynamico_setup_xios() 52 self.comm=MPI.COMM_WORLD # FIXME : we should use the communicator returned by XIOS 53 return self 54 def __exit__(self, type, value, traceback): 55 print 'xios_finalize()' 56 cxios.finalize() 57 58 class Context: 59 def __enter__(self): 60 print 'cxios.context_close_definition()' 61 cxios.context_close_definition() 62 return self 63 def __exit__(self, type, value, traceback): 64 print 'xios_context_finalize()' 65 cxios.context_finalize() 47 66 def init_llm(self, mesh, nqtot): 48 67 self.mesh=mesh 49 ker.dynamico_setup_xios()50 68 # vertical axis, nq 51 69 llm = mesh.llm … … 65 83 def update_calendar(self, i): 66 84 ker.dynamico_xios_update_calendar(i) 67 def finalize(self): 68 cxios.context_finalize() 85 def send_field_primal(self,name, data): 86 own_loc = self.mesh.primal_own_loc 87 if data.ndim==1: 88 data = data[own_loc] 89 if data.ndim==2: 90 data = data[own_loc,:] 91 data = data.transpose() # XIOS expects contiguous horizontal slices 92 data = np.ascontiguousarray(data) 93 cxios.send_field(name, data) 69 94 70 class XIOS_Context_Curvilinear(XIOS_Context):95 class Context_Curvilinear(Context): 71 96 def __init__(self, mesh, nqtot, step): 72 97 self.init_llm(mesh, nqtot) 73 98 # primal mesh 74 # lat_i,lon_i = np.meshgrid(mesh.y,mesh.x, indexing='ij')75 # ny,nx=lat_i.shape76 99 own = mesh.primal_own_loc 77 100 lon_i, lat_i = [ x[own]*radian for x in mesh.lon_i, mesh.lat_i ] 78 101 primal_i, primal_j = [ x[own] for x in mesh.primal_i, mesh.primal_j ] 79 print 'setup_curvilinear : ', lon_i.size, lat_i.size, primal_i.size, primal_j.size80 102 setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.displ, 81 103 primal_i, primal_j, lon_i, lat_i) 82 104 # calendar 83 105 self.init_calendar(step) 84 print 'cxios.context_close_definition()'85 cxios.context_close_definition()86 def send_field_primal(self,name, data):87 if data.ndim==2:88 data = data.transpose() # XIOS expects contiguous horizontal slices89 data = np.ascontiguousarray(data)90 cxios.send_field(name, data)91 106 92 class XIOS_Context_Unstructured(XIOS_Context):107 class Context_Unstructured(Context): 93 108 def __init__(self, pmesh, mesh,nqtot, step): 94 109 self.init_llm(mesh, nqtot) … … 101 116 # calendar 102 117 self.init_calendar(step) 103 print 'cxios.context_close_definition()'104 cxios.context_close_definition()105 def send_field_primal(self,name, data):106 own_loc = self.mesh.primal_own_loc107 if data.ndim==1:108 data = data[own_loc]109 if data.ndim==2:110 data = data[own_loc,:]111 data = data.transpose() # XIOS expects contiguous horizontal slices112 data = np.ascontiguousarray(data)113 cxios.send_field(name, data) -
codes/icosagcm/devel/Python/src/cxios.pyx
r759 r760 28 28 cdef char* idptr = id 29 29 cdef int idlen = len(id) 30 print 'cxios.write_data', id, data.shape30 # print 'cxios.write_data', id, data.shape 31 31 ndim=data.ndim 32 32 if ndim==1: send_field1(idptr,idlen,data) … … 61 61 create_fun(byref(self.handle), id, c_int(len(id))) 62 62 self.cat, self.prefix_set, self.id = cat, 'set_'+cat+'_', id 63 print 'Handle', cat, self.id, self.handle63 # print 'Handle', cat, self.id, self.handle 64 64 def set_attr(self, **kwargs): 65 65 prefix, handle = self.prefix_set, self.handle -
codes/icosagcm/devel/Python/test/py/DCMIP2008c5_MPI.py
r701 r760 58 58 # mesh = meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 59 59 pmesh = meshes.Unstructured_PMesh(comm,meshfile) 60 pmesh.partition_metis() 60 61 mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, f) 61 62 mesh.pmesh=pmesh -
codes/icosagcm/devel/Python/test/py/RSW2_MPAS_W02.py
r749 r760 29 29 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 30 30 pmesh = PMesh(comm,meshfile) 31 pmesh.partition_metis() 31 32 mesh = Mesh(pmesh, llm, nqdyn, radius, f) 32 33 print '...Done' -
codes/icosagcm/devel/Python/test/py/partition.py
r692 r760 80 80 meshfile = meshes.MPAS_Format('grids/%s.grid.nc'%grid) 81 81 pmesh = meshes.Unstructured_PMesh(comm, meshfile) 82 pmesh.partition_metis() 82 83 83 84 def coriolis(lon,lat): return 0.*lat -
codes/icosagcm/devel/Python/test/py/test_xios.py
r759 r760 1 import dynamico.unstructured as unst2 import dynamico.xios asxios1 from dynamico import unstructured as unst 2 from dynamico import xios 3 3 from dynamico.meshes import radian, MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh 4 import numpy as np 4 5 5 from mpi4py import MPI 6 comm = MPI.COMM_WORLD 7 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 8 print '%d/%d starting'%(mpi_rank,mpi_size) 9 10 import numpy as np 11 import cProfile 12 import re 6 with xios.Client() as client: 7 comm = client.comm 8 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 9 print '%d/%d starting'%(mpi_rank,mpi_size) 13 10 14 11 #----------------------------- read MPAS mesh -------------------------------- 15 12 16 grid, llm, nqdyn, nqtot = 2562, 4,1,1 # 2562, 10242, 4096217 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e418 N, T, courant = 40, 10800., 1.2 # simulation length = N*T19 print 'Omega, planetary PV', Omega, 2*Omega/gh013 grid, llm, nqdyn, nqtot = 2562, 4,1,1 # 2562, 10242, 40962 14 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 15 N, T, courant = 40, 10800., 1.2 # simulation length = N*T 16 print 'Omega, planetary PV', Omega, 2*Omega/gh0 20 17 21 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 22 print 'Reading MPAS mesh ...' 23 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 24 pmesh = PMesh(comm,meshfile) 25 mesh = Mesh(pmesh, llm, nqdyn, radius, f) 26 print '...Done' 18 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter 19 print 'Reading MPAS mesh ...' 20 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 21 pmesh = PMesh(comm,meshfile) 22 pmesh.partition_metis() 23 mesh = Mesh(pmesh, llm, nqdyn, radius, f) 24 print '...Done' 27 25 28 26 #--------------------------------- write some data ---------------------------------------- 29 27 30 context=xios.XIOS_Context_Unstructured(pmesh,mesh,nqtot, 3600) 28 with xios.Context_Unstructured(pmesh,mesh,nqtot, 3600) as context: 29 lat_i = radian*mesh.lat_i 30 lat_ik, junk = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 31 31 32 lat_i = radian*mesh.lat_i 33 lat_ik, junk = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 34 35 no_error=True 36 try: 37 for i in range(100): 38 context.update_calendar(i) 39 print 'send_field', i, lat_ik.shape 40 context.send_field_primal('ps', lat_i) 41 context.send_field_primal('theta', lat_ik) 42 except: 43 print 'An error has occured, trying to close XIOS properly' 44 no_error=False 45 46 print 'xios.context_finalize()' 47 context.finalize() 48 print 'xios.finalize()' 49 xios.finalize() 50 print 'Done' 51 52 if not no_error : raise 32 for i in range(100): 33 context.update_calendar(i) 34 print 'send_field', i, lat_ik.shape 35 context.send_field_primal('ps', lat_i) 36 context.send_field_primal('theta', lat_ik) -
codes/icosagcm/devel/Python/test/py/test_xios_cartesian.py
r759 r760 3 3 from dynamico import unstructured as unst 4 4 from dynamico import meshes 5 ker=unst.ker6 7 from mpi4py import MPI8 comm = MPI.COMM_WORLD9 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()10 print '%d/%d starting'%(mpi_rank,mpi_size)11 12 5 import numpy as np 13 6 14 nqdyn,nqtot,llm = 1,1,4 7 # client.__exit__() and context.__exit__() are called at the end of the 'with' blocks 8 # even if an exception is raised in the body of the block 15 9 16 if False: 17 nx,ny,llm = 64,64 18 Lx,Ly,f = 1.,1.,1. 19 mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) 20 else: 21 filename, g, f, radius = 'cart_128_128.nc', 1., 1., None 22 filename = 'cart_004_004.nc' 23 unst.setvar('g',g) 24 print 'Reading Cartesian mesh ...' 25 def coriolis(lon,lat): 26 return f+0.*lon 27 meshfile = meshes.DYNAMICO_Format(filename) 28 nx,ny=meshfile.nx,meshfile.ny 10 with xios.Client() as client: 11 comm = client.comm 12 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 13 print '%d/%d starting'%(mpi_rank,mpi_size) 29 14 30 pmesh = meshes.Unstructured_PMesh(comm,meshfile) 31 mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 15 nqdyn,nqtot,llm = 1,1,4 16 17 if False: 18 nx,ny,llm = 64,64 19 Lx,Ly,f = 1.,1.,1. 20 mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) 21 else: 22 filename, g, f, radius = 'cart_128_128.nc', 1., 1., None 23 # filename = 'cart_008_008.nc' 24 unst.setvar('g',g) 25 print 'Reading Cartesian mesh ...' 26 def coriolis(lon,lat): 27 return f+0.*lon 28 meshfile = meshes.DYNAMICO_Format(filename) 29 nx,ny=meshfile.nx,meshfile.ny 30 31 pmesh = meshes.Unstructured_PMesh(comm,meshfile) 32 # pmesh.partition_metis() 33 pmesh.partition_curvilinear(2,2) 34 mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 32 35 33 36 #--------------------------------- write some data ---------------------------------------- 34 37 35 context=xios.XIOS_Context_Curvilinear(mesh,nqtot, 3600) 36 #lon_i,lat_i = mesh.xx[:,:,0].flatten(), mesh.yy[:,:,0].flatten() 37 lon_i,lat_i = mesh.lon_i, mesh.lat_i 38 lat_ik, junk = np.meshgrid(lat_i, np.arange(llm), indexing='ij')38 with xios.Context_Curvilinear(mesh,nqtot, 3600) as context: 39 lon_i,lat_i = mesh.lon_i, mesh.lat_i 40 lat_ik, junk = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 41 lon_ik, junk = np.meshgrid(lon_i, np.arange(llm), indexing='ij') 39 42 40 no_error=True 41 try: 42 for i in range(3): 43 context.update_calendar(i) 44 print 'send_field', i, lat_ik.shape 45 context.send_field_primal('ps', lat_i) 46 context.send_field_primal('theta', lat_ik) 47 except: 48 print 'An error has occured, trying to close XIOS properly' 49 no_error=False 50 51 print 'xios.context_finalize()' 52 context.finalize() 53 print 'xios.finalize()' 54 xios.finalize() 55 print 'Done' 56 57 if not no_error : raise 43 for i in range(48): 44 context.update_calendar(i) 45 # print 'send_field', i, lat_ik.shape 46 context.send_field_primal('ps', lat_i+lon_i) 47 context.send_field_primal('theta', lat_ik+lon_ik) -
codes/icosagcm/devel/Python/test/xml/iodef.xml
r755 r760 3 3 <context id="xios"> 4 4 <variable_definition> 5 <variable id="print_file" type="bool"> false </variable>5 <variable id="print_file" type="bool"> true </variable> 6 6 <variable_group id="buffer"> 7 7 </variable_group>
Note: See TracChangeset
for help on using the changeset viewer.