- Timestamp:
- 10/10/18 18:12:01 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r757 r758 81 81 deg[ij]=k 82 82 83 class Cartesian_mesh: 83 def concat(x,y): 84 z = np.asarray([x,y]) 85 return z.transpose() 86 87 class Cartesian_mesh(Base_class): 84 88 def __init__(self,nx,ny,llm,nqdyn,Lx,Ly,f): 85 89 dx,dy = Lx/nx, Ly/ny … … 103 107 indexu=lambda x,y : 2*index(x,y)-1 104 108 indexv=lambda x,y : 2*index(x,y) 105 indices = lambda shape : np.zeros(shape,np.int32) 106 107 primal_nb = indices(nx*ny) 108 primal_edge = indices((nx*ny,4)) 109 primal_ne = indices((nx*ny,4)) 109 def zeros(shape,n=1): return [np.zeros(shape) for i in range(n)] if n>1 else np.zeros(shape) 110 def indices(shape,n=1): return [np.zeros(shape,np.int32) for i in range(n)] if n>1 else np.zeros(shape,np.int32) 111 112 Aiv, lon_i,lon_v,lat_i,lat_v = zeros( nx*ny, 5) 113 lon_e,lat_e,le,de,angle_e = zeros( 2*nx*ny, 5) 114 Riv2 = zeros((nx*ny,4) ) 115 wee = zeros((2*nx*ny,4)) 116 117 primal_deg,dual_deg = indices( nx*ny, 2) 118 primal_edge,primal_ne,primal_vertex = indices((nx*ny,4),3) 119 dual_edge,dual_ne,dual_vertex = indices((nx*ny,4),3) 120 trisk_deg,left,right,up,down = indices( 2*nx*ny, 5) 121 trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge 110 122 111 dual_nb = indices(nx*ny) 112 dual_edge = indices((nx*ny,4)) 113 dual_ne = indices((nx*ny,4)) 114 dual_vertex = indices((nx*ny,4)) 115 Riv2 = np.zeros((nx*ny,4)) 116 117 left = indices(2*nx*ny) 118 right = indices(2*nx*ny) 119 up = indices(2*nx*ny) 120 down = indices(2*nx*ny) 121 le_de = np.zeros(2*nx*ny) 122 123 trisk_deg = indices(2*nx*ny) 124 trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge 125 wee = np.zeros((2*nx*ny,4)) 126 123 def putval(ij, arrays, vals): 124 for array,v in zip(arrays,vals): array[ij]=v 127 125 for x in range(nx): 128 126 for y in range(ny): 129 127 # NB : Fortran indices start at 1 130 128 # primal cells 131 put(index(x,y)-1,primal_nb,(primal_edge,primal_ne), 132 ((indexu(x,y),1), 133 (indexv(x,y),1), 134 (indexu(x-1,y),-1), 135 (indexv(x,y-1),-1) )) 129 ij=index(x,y)-1 130 lon_i[ij], lat_i[ij] = x, y 131 put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne), 132 ((indexu(x,y),index(x,y),1), 133 (indexv(x,y),index(x-1,y),1), 134 (indexu(x-1,y),index(x-1,y-1),-1), 135 (indexv(x,y-1),index(x,y-1),-1) )) 136 136 # dual cells 137 put(index(x,y)-1,dual_nb,(dual_edge,dual_vertex,dual_ne,Riv2), 137 ij=index(x,y)-1 138 lon_v[ij], lat_v[ij] = x+.5, y+.5 139 put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2), 138 140 ((indexv(x+1,y),index(x,y),1,.25), 139 141 (indexu(x,y+1),index(x+1,y),-1,.25), … … 147 149 # u-points 148 150 ij =indexu(x,y)-1 149 left[ij]=index(x,y) 150 right[ij]=index(x+1,y) 151 down[ij]=index(x,y-1) 152 up[ij]=index(x,y) 153 le_de[ij]=dy/dx 151 putval(ij, (left, right, down, up, le, de, angle_e, lon_e, lat_e), 152 (index(x,y), index(x+1,y), index(x,y-1), index(x,y), dy, dx, 0., x+.5, y ) ) 154 153 put(ij,trisk_deg,(trisk,wee),( 155 (indexv(x,y), -.25),156 (indexv(x+1,y), -.25),157 (indexv(x,y-1), -.25),158 (indexv(x+1,y-1), -.25)))154 (indexv(x,y),.25), 155 (indexv(x+1,y),.25), 156 (indexv(x,y-1),.25), 157 (indexv(x+1,y-1),.25))) 159 158 # v-points 160 159 ij = indexv(x,y)-1 161 left[ij]=index(x,y) 162 right[ij]=index(x,y+1) 163 down[ij]=index(x,y) 164 up[ij]=index(x-1,y) 165 le_de[ij]=dx/dy 160 putval(ij, (left, right, down, up, le, de, angle_e, lon_e, lat_e), 161 (index(x,y), index(x,y+1), index(x,y), index(x-1,y), dx, dy, 5.*math.pi, x, y+.5 ) ) 166 162 put(ij,trisk_deg,(trisk,wee),( 167 (indexu(x,y),.25), 168 (indexu(x-1,y),.25), 169 (indexu(x,y+1),.25), 170 (indexu(x-1,y+1),.25))) 171 Aiv=np.zeros(nx*ny) 163 (indexu(x,y),-.25), 164 (indexu(x-1,y),-.25), 165 (indexu(x,y+1),-.25), 166 (indexu(x-1,y+1),-.25))) 167 168 primal_i, primal_j = [ x.astype(np.int32) for x in lon_i, lat_i] 169 lon_i, lon_v, lon_e = [x*dx-Lx/2 for x in lon_i, lon_v, lon_e] 170 lat_i, lat_v, lat_e = [y*dy-Ly/2 for y in lat_i, lat_v, lat_e] 171 172 172 Aiv[:]=dx*dy # Ai=Av=dx*dy 173 le_de=le/de 173 174 unst.init_mesh(llm,nqdyn,2*nx*ny,nx*ny,nx*ny,4,4,4, 174 primal_ nb,primal_edge,primal_ne,175 dual_ nb,dual_edge,dual_ne,dual_vertex,175 primal_deg,primal_edge,primal_ne, 176 dual_deg,dual_edge,dual_ne,dual_vertex, 176 177 left,right,down,up,trisk_deg,trisk, 177 178 Aiv,Aiv,f+0.*Aiv,le_de,Riv2,-wee) 179 self.set_members(locals(), 'primal_deg', 'primal_edge', 'primal_ne', 'primal_vertex', 180 'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 181 'left','right','down','up', 182 'trisk_deg','trisk','Riv2','wee', 183 'Aiv','lon_i','lat_i','lon_v','lat_v', 184 'le_de','le','de','lon_e','lat_e','angle_e') 185 def ncwrite(self, name): 186 """The method writes Cartesian mesh on the disc. 187 Args: Mesh parameters""" 188 189 #----------------OPENING NETCDF OUTPUT FILE------------ 190 191 try: 192 f = cdf.Dataset(name, 'w', format='NETCDF4') 193 except: 194 print("Cartesian_mesh.ncwrite : Error occurred while opening new netCDF mesh file") 195 raise 196 197 #----------------DEFINING DIMENSIONS-------------------- 198 199 for dimname, dimsize in [("primal_cell", self.primal_deg.size), 200 ("dual_cell", self.dual_deg.size), 201 ("edge", self.left.size), 202 ("primal_edge_or_vertex", self.primal_edge.shape[1]), 203 ("dual_edge_or_vertex", self.dual_edge.shape[1]), 204 ("TWO", 2), 205 ("trisk_edge", 4)]: 206 f.createDimension(dimname,dimsize) 207 208 #----------------DEFINING VARIABLES--------------------- 209 210 f.description = "Cartesian_mesh" 211 f.nx, f.ny = self.nx, self.ny 212 213 def create_vars(dimname, info): 214 for vname, vtype, vdata in info: 215 var = f.createVariable(vname,vtype,dimname) 216 var[:] = vdata 217 218 create_vars("primal_cell", 219 [("primal_deg","i4",self.primal_deg), 220 ("Ai","f8",self.Aiv), 221 ("lon_i","f8",self.lon_i), 222 ("lat_i","f8",self.lat_i)] ) 223 create_vars("dual_cell", 224 [("dual_deg","i4",self.dual_deg), 225 ("Av","f8",self.Aiv), 226 ("lon_v","f8",self.lon_v), 227 ("lat_v","f8",self.lat_v), 228 ] ) 229 230 create_vars( ("primal_cell","primal_edge_or_vertex"), 231 [("primal_edge", "i4", self.primal_edge), 232 ("primal_ne", "i4", self.primal_ne), 233 ("primal_vertex", "i4", self.primal_vertex)] ) 234 235 create_vars( ("dual_cell","dual_edge_or_vertex"), 236 [("dual_edge", "i4", self.dual_edge), 237 ("dual_vertex","i4",self.dual_vertex), 238 ("dual_ne","i4",self.dual_ne), 239 ("Riv2","f8",self.Riv2)] ) 240 241 create_vars("edge", 242 [("trisk_deg","i4",self.trisk_deg), 243 ("le","f8",self.le), 244 ("de","f8",self.de), 245 ("lon_e","f8", self.lon_e), 246 ("lat_e","f8", self.lat_e), 247 ("angle_e","f8", self.angle_e)] ) 248 249 create_vars( ("edge","TWO"), 250 [("edge_left_right","i4", concat(self.left,self.right)), 251 ("edge_down_up","i4", concat(self.down,self.up))] ) 252 253 254 create_vars( ("edge","trisk_edge"), 255 [("trisk","i4",self.trisk), 256 ("wee","f8",self.wee)] ) 257 258 f.close() 178 259 179 260 def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm)) … … 371 452 372 453 class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes 373 def __init__(self,comm, gridfile ):454 def __init__(self,comm, gridfile, meshtype='unstructured'): 374 455 self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars 375 456 dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, 376 'primal_ num','edge_num','dual_num')457 'primal_cell','edge','dual_cell') 377 458 primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars( 378 459 dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' ) … … 397 478 primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 398 479 primal_owner = parallel.LocPArray1D(dim_primal, primal_owner) 399 self.set_members(locals(), ' dim_primal', 'dim_dual', 'dim_edge',480 self.set_members(locals(), 'meshtype', 'dim_primal', 'dim_dual', 'dim_edge', 400 481 'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex', 401 482 'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 402 483 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 484 # if meshtype == 'curvilinear' : # read extra information from mesh file 485 # self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 486 403 487 def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ] 404 488 … … 445 529 get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 446 530 primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2)) 531 532 if pmesh.meshtype == 'curvilinear' : 533 self.primal_i, self.primal_j = pmesh.get(get_own_cells, 'primal_i', 'primal_j') 447 534 448 535 # construct local stencils from global stencils … … 655 742 # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges 656 743 return list(progressive_iter([], cell_lists)) 744 from dynamico.meshes import zeros 745 import numpy as np 746 import netCDF4 as cdf 747 import argparse 748 749 750 if False: 751 def __init__(self,nx,ny,llm,nqdyn,Lx,Ly): 752 dx,dy = Lx/nx, Ly/ny 753 self.dx, self.dy = dx,dy 754 self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn 755 self.field_z = self.field_mass 756 # 1D coordinate arrays 757 x=(np.arange(nx)-nx/2.)*Lx/nx 758 y=(np.arange(ny)-ny/2.)*Ly/ny 759 lev=np.arange(llm) 760 levp1=np.arange(llm+1) 761 self.x, self.y, self.lev, self.levp1 = x,y,lev,levp1 762 # 3D coordinate arrays 763 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 indexing 766 # Fortran order : llm,nx*ny,nqdyn / indices start at 1 767 # Python order : nqdyn,ny,nx,llm / indices start at 0 768 # indices below follow Fortran while x,y follow Python/C 769 index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+1 770 indexu=lambda x,y : 2*index(x,y)-1 771 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 1 777 # primal cells 778 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 cells 784 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 cells 791 # flux is positive when going from left to right 792 # up and down are adjacent dual cells (vertices) 793 # circulation is positive when going from down to up 794 795 Aiv=np.zeros(nx*ny)+dx*dy 796 Ai=Av=np.zeros(nx*ny)+dx*dy 797 798 self.llm=llm 799 self.nqdyn=nqdyn 800 self.nx=nx 801 self.ny=ny 802 self.primal_deg=primal_deg 803 self.primal_edge=primal_edge 804 self.primal_ne=primal_ne 805 self.dual_deg=dual_deg 806 self.dual_edge=dual_edge 807 self.dual_ne=dual_ne 808 self.dual_vertex=dual_vertex 809 self.primal_vertex=primal_vertex 810 self.left=left 811 self.right=right 812 self.down=down 813 self.up=up 814 self.trisk_deg=trisk_deg 815 self.trisk=trisk 816 self.Aiv=Aiv 817 self.Ai=Ai 818 self.Av=Av 819 self.le=le 820 self.de=de 821 self.angle_e=angle_e 822 self.Riv2=Riv2 823 self.wee=wee 824 self.lon_i = lon_i 825 self.lon_v = lon_v 826 self.lon_e = lon_e 827 #self.lon_ev = indices(2*nx*ny) 828 self.lat_i = lat_i 829 self.lat_v = lat_v 830 self.lat_e = lat_e 831 #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),:]=u 846 def vcomp(self,u): return u[:,range(1,2*self.nx,2),:] 847 848 -
codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py
r751 r758 1 1 from dynamico.meshes import zeros 2 from dynamico import meshes 2 3 import numpy as np 3 4 import netCDF4 as cdf 4 5 import argparse 5 6 6 #----------------------- Cartesian mesh -----------------------7 8 def concat(x,y):9 z = np.asarray([x,y])10 return z.transpose()11 12 # arrays is a list of arrays13 # vals is a list of tuples14 # each tuple is stored in each array15 def put(ij, deg, arrays, vals):16 k=017 for vv in vals: # vv is a tuple of values to be stored in arrays18 for array,v in zip(arrays,vv):19 array[ij,k]=v20 k=k+121 deg[ij]=k22 23 class Cartesian_mesh:24 def __init__(self,nx,ny,llm,nqdyn,Lx,Ly):25 dx,dy = Lx/nx, Ly/ny26 self.dx, self.dy = dx,dy27 self.nx, self.ny, self.llm, self.nqdyn = nx,ny,llm,nqdyn28 self.field_z = self.field_mass29 # 1D coordinate arrays30 x=(np.arange(nx)-nx/2.)*Lx/nx31 y=(np.arange(ny)-ny/2.)*Ly/ny32 lev=np.arange(llm)33 levp1=np.arange(llm+1)34 self.x, self.y, self.lev, self.levp1 = x,y,lev,levp135 # 3D coordinate arrays36 self.xx,self.yy,self.ll = np.meshgrid(x,y,lev, indexing='ij')37 self.xxp1,self.yyp1,self.llp1 = np.meshgrid(x,y,levp1, indexing='ij')38 # beware conventions for indexing39 # Fortran order : llm,nx*ny,nqdyn / indices start at 140 # Python order : nqdyn,ny,nx,llm / indices start at 041 # indices below follow Fortran while x,y follow Python/C42 index=lambda x,y : ((x+(nx*(y+2*ny)))%(nx*ny))+143 indexu=lambda x,y : 2*index(x,y)-144 indexv=lambda x,y : 2*index(x,y)45 indices = lambda shape : np.zeros(shape,np.int32)46 47 primal_deg = indices(nx*ny)48 primal_edge = indices((nx*ny,4))49 primal_ne = indices((nx*ny,4))50 51 dual_deg = indices(nx*ny)52 dual_edge = indices((nx*ny,4))53 dual_ne = indices((nx*ny,4))54 primal_vertex = indices((nx*ny,4))55 dual_vertex = indices((nx*ny,4))56 Riv2 = np.zeros((nx*ny,4))57 58 left = indices(2*nx*ny)59 right = indices(2*nx*ny)60 up = indices(2*nx*ny)61 down = indices(2*nx*ny)62 le_de = np.zeros(2*nx*ny)63 le = np.zeros(2*nx*ny)64 de = np.zeros(2*nx*ny)65 angle_e = np.zeros(2*nx*ny)66 67 trisk_deg = indices(2*nx*ny)68 trisk = indices((2*nx*ny,4)) # 4 TRiSK coefs per edge69 wee = np.zeros((2*nx*ny,4))70 71 lon_i = indices(nx*ny)72 lon_v = indices(nx*ny)73 lon_e = indices(2*nx*ny)74 lat_i = indices(nx*ny)75 lat_v = indices(nx*ny)76 lat_e = indices(2*nx*ny)77 78 for x in range(nx):79 for y in range(ny):80 # NB : Fortran indices start at 181 # primal cells82 put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne),83 ((indexu(x,y),index(x,y),1),84 (indexv(x,y),index(x-1,y),1),85 (indexu(x-1,y),index(x-1,y-1),-1),86 (indexv(x,y-1),index(x,y-1),-1) ))87 # dual cells88 put(index(x,y)-1,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2),89 ((indexv(x+1,y),index(x,y),1,.25),90 (indexu(x,y+1),index(x+1,y),-1,.25),91 (indexv(x,y),index(x+1,y+1),-1,.25),92 (indexu(x,y),index(x,y+1),1,.25) ))93 # edges :94 # left and right are adjacent primal cells95 # flux is positive when going from left to right96 # up and down are adjacent dual cells (vertices)97 # circulation is positive when going from down to up98 # u-points99 ij =indexu(x,y)-1100 left[ij]=index(x,y)101 left[ij]=index(x,y)102 right[ij]=index(x+1,y)103 down[ij]=index(x,y-1)104 up[ij]=index(x,y)105 #le_de[ij]=dy/dx106 le[ij]=dy107 de[ij]=dx108 le_de[ij]=le[ij]/de[ij]109 angle_e[ij]=0.110 put(ij,trisk_deg,(trisk,wee),(111 (indexv(x,y),.25),112 (indexv(x+1,y),.25),113 (indexv(x,y-1),.25),114 (indexv(x+1,y-1),.25)))115 # v-points116 ij = indexv(x,y)-1117 left[ij]=index(x,y)118 right[ij]=index(x,y+1)119 down[ij]=index(x,y)120 up[ij]=index(x-1,y)121 #le_de[ij]=dx/dy122 le[ij]=dx123 de[ij]=dy124 le_de[ij]=le[ij]/de[ij]125 angle_e[ij]=.5*np.pi126 put(ij,trisk_deg,(trisk,wee),(127 (indexu(x,y),-.25),128 (indexu(x-1,y),-.25),129 (indexu(x,y+1),-.25),130 (indexu(x-1,y+1),-.25)))131 ij = index(x,y)-1132 lon_i[ij], lat_i[ij] = x, y133 ij = index(x,y)-1134 lon_v[ij], lat_v[ij] = x+.5, y+.5135 ij = indexu(x,y)-1136 lon_e[ij], lat_e[ij] = x+.5, y137 ij = indexv(x,y)-1138 lon_e[ij], lat_e[ij] = x, y+.5139 140 Aiv=np.zeros(nx*ny)+dx*dy141 Ai=Av=np.zeros(nx*ny)+dx*dy142 143 self.llm=llm144 self.nqdyn=nqdyn145 self.nx=nx146 self.ny=ny147 self.primal_deg=primal_deg148 self.primal_edge=primal_edge149 self.primal_ne=primal_ne150 self.dual_deg=dual_deg151 self.dual_edge=dual_edge152 self.dual_ne=dual_ne153 self.dual_vertex=dual_vertex154 self.primal_vertex=primal_vertex155 self.left=left156 self.right=right157 self.down=down158 self.up=up159 self.trisk_deg=trisk_deg160 self.trisk=trisk161 self.Aiv=Aiv162 self.Ai=Ai163 self.Av=Av164 self.le=le165 self.de=de166 self.angle_e=angle_e167 self.Riv2=Riv2168 self.wee=wee169 self.lon_i = lon_i170 self.lon_v = lon_v171 self.lon_e = lon_e172 #self.lon_ev = indices(2*nx*ny)173 self.lat_i = lat_i174 self.lat_v = lat_v175 self.lat_e = lat_e176 #self.lat_ev = indices(2*nx*ny)177 178 179 def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))180 def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))181 def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))182 def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))183 def field_u(self,n=1):184 if n==1:185 return np.zeros((self.ny,2*self.nx,self.llm))186 else:187 return np.zeros((n,self.ny,2*self.nx,self.llm))188 def field_ps(self,n=1): return zeros((n,self.ny,self.nx))189 def ucomp(self,u): return u[:,range(0,2*self.nx,2),:]190 def set_ucomp(self,uv,u): uv[:,range(0,2*self.nx,2),:]=u191 def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]192 193 194 #----------------------------WRITING MESH----------------------195 def ncwrite(self, name):196 """The method writes Cartesian mesh on the disc.197 Args: Mesh parameters"""198 199 #----------------OPENING NETCDF OUTPUT FILE------------200 201 try:202 f = cdf.Dataset(name, 'w', format='NETCDF4')203 except:204 print("CartesianMesh.ncwrite : Error occurred while opening new netCDF mesh file")205 raise206 207 #----------------DEFINING DIMENSIONS--------------------208 209 for dimname, dimsize in [("primal_cell", self.primal_deg.size),210 ("dual_cell", self.dual_deg.size),211 ("edge", self.left.size),212 ("primal_edge_or_vertex", self.primal_edge.shape[1]),213 ("dual_edge_or_vertex", self.dual_edge.shape[1]),214 ("TWO", 2),215 ("trisk_edge", 4)]:216 f.createDimension(dimname,dimsize)217 218 #----------------DEFINING VARIABLES---------------------219 220 f.description = "Cartesian_mesh"221 f.nx, f.ny = self.nx, self.ny222 223 def create_vars(dimname, info):224 for vname, vtype, vdata in info:225 # print('create_vars', dimname, vname, vdata.shape)226 var = f.createVariable(vname,vtype,dimname)227 var[:] = vdata228 229 create_vars("primal_cell",230 [("primal_deg","i4",self.primal_deg),231 ("Ai","f8",self.Aiv),232 ("lon_i","f8",self.lon_i),233 ("lat_i","f8",self.lat_i)] )234 create_vars("dual_cell",235 [("dual_deg","i4",self.dual_deg),236 ("Av","f8",self.Aiv),237 ("lon_v","f8",self.lon_v),238 ("lat_v","f8",self.lat_v),239 ] )240 241 create_vars( ("primal_cell","primal_edge_or_vertex"),242 [("primal_edge", "i4", self.primal_edge),243 ("primal_ne", "i4", self.primal_ne),244 ("primal_vertex", "i4", self.primal_vertex)] )245 246 create_vars( ("dual_cell","dual_edge_or_vertex"),247 [("dual_edge", "i4", self.dual_edge),248 ("dual_vertex","i4",self.dual_vertex),249 ("dual_ne","i4",self.dual_ne),250 ("Riv2","f8",self.Riv2)] )251 252 create_vars("edge",253 [("trisk_deg","i4",self.trisk_deg),254 ("le","f8",self.le),255 ("de","f8",self.de),256 ("lon_e","f8", self.lon_e),257 ("lat_e","f8", self.lat_e),258 ("angle_e","f8", self.angle_e)] )259 260 create_vars( ("edge","TWO"),261 [("edge_left_right","i4", concat(self.left,self.right)),262 ("edge_down_up","i4", concat(self.down,self.up))] )263 264 265 create_vars( ("edge","trisk_edge"),266 [("trisk","i4",self.trisk),267 ("wee","f8",self.wee)] )268 269 f.close()270 271 272 #--------------------------------- Main program -------------------------------------273 274 7 parser = argparse.ArgumentParser() 275 8 276 9 parser.add_argument("-nx", type=int, 277 default= 128, choices=None,10 default=64, choices=None, 278 11 help="number of x points") 279 12 parser.add_argument("-ny", type=int, 280 default= 128, choices=None,13 default=64, choices=None, 281 14 help="number of y points") 282 15 parser.add_argument("-Lx", type=float, … … 293 26 294 27 dx,dy=Lx/nx,Ly/ny 295 mesh = Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly)28 mesh = meshes.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,Ly,0.) 296 29 print('Successfully initialized Cartesian mesh') 297 mesh.ncwrite('cart_ '+str(nx)+'.nc')30 mesh.ncwrite('cart_%03d_%03d.nc'%(nx,ny)) 298 31 print('Successfully written Cartesian mesh to NetCDF File')
Note: See TracChangeset
for help on using the changeset viewer.