- Timestamp:
- 10/11/18 18:27:15 (6 years ago)
- Location:
- codes/icosagcm/devel
- Files:
-
- 1 added
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r758 r759 182 182 'trisk_deg','trisk','Riv2','wee', 183 183 'Aiv','lon_i','lat_i','lon_v','lat_v', 184 'le_de','le','de','lon_e','lat_e','angle_e') 184 'le_de','le','de','lon_e','lat_e','angle_e', 185 'primal_i','primal_j') 185 186 def ncwrite(self, name): 186 187 """The method writes Cartesian mesh on the disc. … … 209 210 210 211 f.description = "Cartesian_mesh" 211 f. nx, f.ny =self.nx, self.ny212 f.meshtype, f.nx, f.ny = 'curvilinear', self.nx, self.ny 212 213 213 214 def create_vars(dimname, info): … … 218 219 create_vars("primal_cell", 219 220 [("primal_deg","i4",self.primal_deg), 221 ("primal_i","i4",self.primal_i), 222 ("primal_j","i4",self.primal_j), 220 223 ("Ai","f8",self.Aiv), 221 224 ("lon_i","f8",self.lon_i), 222 225 ("lat_i","f8",self.lat_i)] ) 226 223 227 create_vars("dual_cell", 224 228 [("dual_deg","i4",self.dual_deg), … … 293 297 start_index=1 # indexing starts at 1 as in Fortran 294 298 def __init__(self,gridfilename): 295 self.nc = cdf.Dataset(gridfilename, "r") 299 nc = cdf.Dataset(gridfilename, "r") 300 self.nc, self.meshtype = nc, nc.meshtype 301 if self.meshtype == 'curvilinear': 302 self.nx, self.ny = nc.nx, nc.ny 296 303 print('#NC4: Opening file:',gridfilename) 297 304 def get_pdim(self,comm,name): return parallel.PDim(self.nc.dimensions[name], comm) … … 310 317 """ Helps class Unstructured_Mesh to read a MPAS mesh file.""" 311 318 start_index=1 # indexing starts at 1 as in Fortran 319 meshtype='unstructured' 312 320 translate= { 313 'primal_ num':'nCells', 'edge_num':'nEdges', 'dual_num':'nVertices',321 'primal_cell':'nCells', 'edge':'nEdges', 'dual_cell':'nVertices', 314 322 'primal_deg':'nEdgesOnCell', 'trisk_deg':'nEdgesOnEdge', 315 323 'primal_edge':'edgesOnCell', 'dual_edge':'edgesOnVertex', … … 452 460 453 461 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes 454 def __init__(self,comm, gridfile , meshtype='unstructured'):462 def __init__(self,comm, gridfile): 455 463 self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars 456 464 dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, … … 463 471 dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars( 464 472 dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v', 'Av') 473 if gridfile.meshtype == 'curvilinear': 474 self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 475 465 476 # Let indices start at 0 466 477 for x in (primal_vertex, primal_edge, dual_vertex,dual_edge, … … 478 489 primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 479 490 primal_owner = parallel.LocPArray1D(dim_primal, primal_owner) 480 self.set_members(locals(), ' meshtype', 'dim_primal', 'dim_dual', 'dim_edge',491 self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge', 481 492 'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex', 482 493 'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', … … 530 541 primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2)) 531 542 532 if pmesh.meshtype == 'curvilinear' : 533 self.primal_i, self.primal_j = pmesh.get(get_own_cells, 'primal_i', 'primal_j') 543 self.meshtype = pmesh.gridfile.meshtype 544 if self.meshtype == 'curvilinear' : 545 self.nx, self.ny = pmesh.gridfile.nx, pmesh.gridfile.ny 546 self.primal_i, self.primal_j = pmesh.get(get_all_cells, 'primal_i', 'primal_j') 547 534 548 535 549 # construct local stencils from global stencils … … 566 580 primal_own_all = comm.allgather(primal_own_num) 567 581 displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS 582 568 583 print 'local primal mesh :', primal_vertex.min(), primal_vertex.max(), primal_own_all, displ 569 584 … … 742 757 # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges 743 758 return list(progressive_iter([], cell_lists)) 744 from dynamico.meshes import zeros745 import numpy as np746 import netCDF4 as cdf747 import argparse748 749 750 if False:751 def __init__(self,nx,ny,llm,nqdyn,Lx,Ly):752 dx,dy = Lx/nx, Ly/ny753 self.dx, self.dy = dx,dy754 self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn755 self.field_z = self.field_mass756 # 1D coordinate arrays757 x=(np.arange(nx)-nx/2.)*Lx/nx758 y=(np.arange(ny)-ny/2.)*Ly/ny759 lev=np.arange(llm)760 levp1=np.arange(llm+1)761 self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1762 # 3D coordinate arrays763 self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')764 self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')765 # beware conventions for indexing766 # Fortran order : llm,nx*ny,nqdyn / indices start at 1767 # Python order : nqdyn,ny,nx,llm / indices start at 0768 # indices below follow Fortran while x,y follow Python/C769 index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1770 indexu=lambda x,y : 2*index(x,y)-1771 indexv=lambda x,y : 2*index(x,y)772 indices = lambda shape : np.zeros(shape,np.int32)773 774 for x in range(nx):775 for y in range(ny):776 # NB : Fortran indices start at 1777 # primal cells778 put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),779 ((indexu(x,y),index(x,y),1),780 (indexv(x,y),index(x-1,y),1),781 (indexu(x-1,y),index(x-1,y-1),-1),782 (indexv(x,y-1),index(x,y-1),-1) ))783 # dual cells784 put(index(x,y)-1,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),785 ((indexv(x+1,y),index(x,y),1,.25),786 (indexu(x,y+1),index(x+1,y),-1,.25),787 (indexv(x,y),index(x+1,y+1),-1,.25),788 (indexu(x,y),index(x,y+1),1,.25) ))789 # edges :790 # left and right are adjacent primal cells791 # flux is positive when going from left to right792 # up and down are adjacent dual cells (vertices)793 # circulation is positive when going from down to up794 795 Aiv=np.zeros(nx*ny)+dx*dy796 Ai=Av=np.zeros(nx*ny)+dx*dy797 798 self.llm=llm799 self.nqdyn=nqdyn800 self.nx=nx801 self.ny=ny802 self.primal_deg=primal_deg803 self.primal_edge=primal_edge804 self.primal_ne=primal_ne805 self.dual_deg=dual_deg806 self.dual_edge=dual_edge807 self.dual_ne=dual_ne808 self.dual_vertex=dual_vertex809 self.primal_vertex=primal_vertex810 self.left=left811 self.right=right812 self.down=down813 self.up=up814 self.trisk_deg=trisk_deg815 self.trisk=trisk816 self.Aiv=Aiv817 self.Ai=Ai818 self.Av=Av819 self.le=le820 self.de=de821 self.angle_e=angle_e822 self.Riv2=Riv2823 self.wee=wee824 self.lon_i = lon_i825 self.lon_v = lon_v826 self.lon_e = lon_e827 #self.lon_ev = indices(2*nx*ny)828 self.lat_i = lat_i829 self.lat_v = lat_v830 self.lat_e = lat_e831 #self.lat_ev = indices(2*nx*ny)832 833 834 def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))835 def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))836 def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))837 def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))838 def field_u(self,n=1):839 if n==1:840 return np.zeros((self.ny,2*self.nx,self.llm))841 else:842 return np.zeros((n,self.ny,2*self.nx,self.llm))843 def field_ps(self,n=1): return zeros((n,self.ny,self.nx))844 def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]845 def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u846 def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]847 848 -
codes/icosagcm/devel/Python/dynamico/xios.py
r696 r759 4 4 from dynamico.meshes import radian 5 5 from cxios import finalize 6 7 #----------------------------------------------------------------------------- 8 9 def setup_curvilinear(cat,id,nx,ny,displ,index_own_i,index_own_j,lon,lat): 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 13 mesh.set_attr(type='curvilinear', 14 ni_glo=nx, i_index=index_own_i, 15 nj_glo=ny, j_index=index_own_j) 16 mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten()) 6 17 7 18 def cf_boundaries(degree,points,lon,lat): … … 19 30 return bnd_lon, bnd_lat 20 31 21 def setup_domain(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_glo, displ): # cat in 'domain','domaingroup' 32 def setup_unstructured(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_own_glo, displ): # cat in 'domain','domaingroup' 33 # ncell_glo (int) : global number of cells 34 # index_own_glo : global indices of cells own by this MPI process 35 # displ : min of index_own across all MPI processes 22 36 bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v) 23 37 ncell, nvertex = bnd_lon.shape 24 38 mesh = cxios.Handle(cat,id) 25 39 mesh.set_attr(type='unstructured') 26 mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_ glo)40 mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_own_glo) 27 41 mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat) 28 42 mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat) 29 43 44 #----------------------------------------------------------------------------- 45 30 46 class XIOS_Context: 31 def __init__(self, pmesh, mesh,nqtot, step):47 def init_llm(self, mesh, nqtot): 32 48 self.mesh=mesh 33 49 ker.dynamico_setup_xios() … … 38 54 levp1.set_attr(n_glo=llm+1, value=np.arange(llm+1)) 39 55 nq.set_attr(n_glo=nqtot, value=np.arange(nqtot)) 40 # primal mesh 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] 42 mesh_i=setup_domain('domaingroup','i', mesh.primal_own_deg, mesh.primal_vertex, 43 lon_i[mesh.primal_own_loc], lat_i[mesh.primal_own_loc], 44 lon_v, lat_v, 45 pmesh.dim_primal.n, mesh.primal_own_glo, mesh.displ) 46 # calendar 56 def init_calendar(self, step): 47 57 calendar=cxios.Handle('calendar_wrapper') 48 58 dtime = cxios.Duration(second=step) … … 53 63 freq_offset = cxios.Duration(second=0) 54 64 std.set_attr(freq_offset=freq_offset) 65 def update_calendar(self, i): 66 ker.dynamico_xios_update_calendar(i) 67 def finalize(self): 68 cxios.context_finalize() 69 70 class XIOS_Context_Curvilinear(XIOS_Context): 71 def __init__(self, mesh, nqtot, step): 72 self.init_llm(mesh, nqtot) 73 # primal mesh 74 # lat_i,lon_i = np.meshgrid(mesh.y,mesh.x, indexing='ij') 75 # ny,nx=lat_i.shape 76 own = mesh.primal_own_loc 77 lon_i, lat_i = [ x[own]*radian for x in mesh.lon_i, mesh.lat_i ] 78 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.size 80 setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.displ, 81 primal_i, primal_j, lon_i, lat_i) 82 # calendar 83 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 slices 89 data = np.ascontiguousarray(data) 90 cxios.send_field(name, data) 91 92 class XIOS_Context_Unstructured(XIOS_Context): 93 def __init__(self, pmesh, mesh,nqtot, step): 94 self.init_llm(mesh, nqtot) 95 # primal mesh 96 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] 97 setup_unstructured('domaingroup','i', mesh.primal_own_deg, mesh.primal_vertex, 98 lon_i[mesh.primal_own_loc], lat_i[mesh.primal_own_loc], 99 lon_v, lat_v, 100 pmesh.dim_primal.n, mesh.primal_own_glo, mesh.displ) 101 # calendar 102 self.init_calendar(step) 55 103 print 'cxios.context_close_definition()' 56 104 cxios.context_close_definition() … … 64 112 data = np.ascontiguousarray(data) 65 113 cxios.send_field(name, data) 66 def update_calendar(self, i):67 ker.dynamico_xios_update_calendar(i)68 def finalize(self):69 cxios.context_finalize() -
codes/icosagcm/devel/Python/src
- Property svn:ignore
-
old new 1 1 build 2 xios.c2 cxios.c 3 3 unstructured.c
-
- Property svn:ignore
-
codes/icosagcm/devel/Python/src/cxios.pyx
r694 r759 67 67 fun = cxios.vardict[prefix+key] 68 68 extra=value 69 if type(value) in (int,c_int ):69 if type(value) in (int,c_int,np.int64,np.int32): 70 70 fun(handle,c_int(value)) 71 71 elif type(value) is np.ndarray: … … 102 102 ['update_calendar',c_int] ]) 103 103 104 import_set_attr(['domain','domaingroup'], [c_int,'ni_glo','ibegin','ni','nvertex','data_dim'], [str,'type'], 105 [np.ndarray, 'i_index','lonvalue_1d','latvalue_1d','bounds_lat_1d','bounds_lon_1d']) 104 import_set_attr(['domain','domaingroup'], [c_int,'ni_glo','nj_glo','ibegin','jbegin','ni','nvertex','data_dim'], [str,'type'], 105 [np.ndarray, 'i_index','j_index', 106 'lonvalue_1d','latvalue_1d','bounds_lat_1d','bounds_lon_1d']) 106 107 import_set_attr(['field','fieldgroup'], [Duration,'freq_op','freq_offset']) 107 108 import_set_attr(['axis'], [c_int,'n_glo'], [np.ndarray,'value']) -
codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py
r756 r759 15 15 dx,dy=Lx/nx,Ly/ny 16 16 17 filename, llm, nqdyn, g, f, radius = 'cart_ 128.nc', 1, 1, 1., 1., None17 filename, llm, nqdyn, g, f, radius = 'cart_%03d_%03d.nc'%(nx,ny), 1, 1, 1., 1., None 18 18 unst.setvar('g',g) 19 19 … … 27 27 caldyn = unst.Caldyn_RSW(mesh) 28 28 29 xx = (mesh.lat_i-nx/2.)*dx 30 yy = (mesh.lon_i-ny/2.)*dy 29 #xx = (mesh.lat_i-nx/2.)*dx 30 #yy = (mesh.lon_i-ny/2.)*dy 31 xx,yy = mesh.lat_i, mesh.lon_i 31 32 32 33 x1,x2,yy = xx-1., xx+1., yy -
codes/icosagcm/devel/Python/test/py/test_xios.py
r696 r759 28 28 #--------------------------------- write some data ---------------------------------------- 29 29 30 context=xios.XIOS_Context (pmesh,mesh,nqtot, 3600)30 context=xios.XIOS_Context_Unstructured(pmesh,mesh,nqtot, 3600) 31 31 32 32 lat_i = radian*mesh.lat_i -
codes/icosagcm/devel/Python/test/python.sh
r717 r759 8 8 { 9 9 LD_PRELOAD=$PYTHON_PRELOAD python -u $* 10 } 11 12 function cmd_gdb() 13 { 14 set -x 15 LD_PRELOAD=$PYTHON_PRELOAD gdb --args python -u $* 10 16 } 11 17 -
codes/icosagcm/devel/Python/test/xml/filedef_simple.xml
r755 r759 44 44 <variable name="time_frequency" type="string" > 24h </variable> 45 45 46 </file> 47 --> 46 </file> --> 48 47 49 48 </file_definition> -
codes/icosagcm/devel/Python/test/xml/link.sh
r755 r759 7 7 ln -s xml/${name}_$1.xml ${name}.xml 8 8 done 9 ln -s iodef.xml9 ln -s xml/iodef.xml 10 10 } 11 11 -
codes/icosagcm/devel/src/unstructured/data_unstructured.F90
r749 r759 164 164 dual_ne(:,:) = dual_ne_(:,:) 165 165 dual_vertex(:,:) = dual_vertex_(:,:) 166 IF(MINVAL(dual_deg)< 3) THEN167 STOP 'At least one dual cell has less than 3vertices'168 END IF 169 IF(MINVAL(primal_deg)< 3) THEN170 STOP 'At least one primal cell has less than 3vertices'166 IF(MINVAL(dual_deg)<2) THEN 167 STOP 'At least one dual cell has less than 2 vertices' 168 END IF 169 IF(MINVAL(primal_deg)<2) THEN 170 STOP 'At least one primal cell has less than 2 vertices' 171 171 END IF 172 172 left(:)=left_(:)
Note: See TracChangeset
for help on using the changeset viewer.