Changeset 787 for codes/icosagcm
- Timestamp:
- 11/23/18 13:06:47 (6 years ago)
- Location:
- codes/icosagcm/devel/Python
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/dynamico/meshes.py
r779 r787 478 478 if gridfile.meshtype == 'curvilinear': 479 479 self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 480 self.dual_i, self.dual_j = get_pvars(dim_dual, 'primal_i','primal_j') 480 481 self.edge_i, self.edge_j = get_pvars(dim_edge, 'edge_i','edge_j') 481 482 … … 494 495 'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 495 496 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 496 def partition(self, edge_owner): 497 def partition_metis(self): 498 print 'Partitioning with ParMetis...' 499 edge_owner = unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size()) 497 500 self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner) 498 501 primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge) 499 502 self.primal_owner = parallel.LocPArray1D(self.dim_primal, primal_owner) 500 501 def partition_metis(self): 502 print 'Partitioning with ParMetis...' 503 self.partition(unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())) 503 dual_owner = partition_from_stencil(self.edge_owner, self.dual_deg, self.dual_edge) 504 self.dual_owner = parallel.LocPArray1D(self.dim_dual, dual_owner) 504 505 505 506 def partition_curvilinear(self, ni,nj): … … 513 514 self.edge_owner = owner(self.dim_edge, self.edge_i, self.edge_j) 514 515 self.primal_owner = owner(self.dim_primal, self.primal_i, self.primal_j) 516 self.dual_owner = self.primal_owner 515 517 print 'partition_curvilinear %d :'%(self.comm.Get_rank()), self.primal_owner.data 516 518 … … 541 543 class Local_Mesh(Abstract_Mesh): 542 544 def __init__(self, pmesh, llm, nqdyn, radius, f): 543 comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual = pmesh.members(544 'comm', 'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual' )545 comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual, dual_owner = pmesh.members( 546 'comm', 'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual', 'dual_owner' ) 545 547 print 'Constructing halos ...' 546 548 edges_E0 = find_my_cells(edge_owner) … … 564 566 self.nx, self.ny = pmesh.gridfile.nx, pmesh.gridfile.ny 565 567 self.primal_i, self.primal_j = pmesh.get(get_all_cells, 'primal_i', 'primal_j') 568 self.dual_i, self.dual_j = pmesh.get(get_all_duals, 'dual_i', 'dual_j') 566 569 567 570 … … 599 602 primal_own_all = comm.allgather(primal_own_num) 600 603 displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS 601 602 604 print 'local primal mesh :', primal_vertex.min(), primal_vertex.max(), primal_own_all, displ 605 606 # construct own dual mesh for XIOS output 607 dual_own_glo = find_my_cells(dual_owner) 608 dual_own_glo = np.asarray(primal_own_glo, dtype=np.int32) 609 dual_own_loc = [dict_V1[i] for i in dual_own_glo] 603 610 604 611 # compute_ne(...) and normalize(...) expect indices starting at 1 … … 610 617 self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num', 611 618 'primal_own_glo', 'primal_own_loc', 'primal_own_deg', 612 'primal_deg', 'primal_vertex', 'displ', 619 'primal_deg', 'primal_vertex', 'displ', 'dual_own_loc', 613 620 'trisk_deg', 'trisk', 'wee', 614 621 'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2', -
codes/icosagcm/devel/Python/dynamico/xios.py
r783 r787 4 4 from dynamico.meshes import radian 5 5 from mpi4py import MPI 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 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)16 mesh.set_attr(type='curvilinear',17 ni_glo=nx, i_index=index_own_i,18 nj_glo=ny, j_index=index_own_j)19 mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten())20 21 def cf_boundaries(degree,points,lon,lat):22 n, nvertex = len(degree), degree.max()23 bnd_lon, bnd_lat = np.zeros((n,nvertex)),np.zeros((n,nvertex))24 for ij in range(n):25 nb=degree[ij]26 for k in range(nb):27 vertex = points[ij,k]28 bnd_lon[ij,k]=lon[vertex]29 bnd_lat[ij,k]=lat[vertex]30 for k in range(nb,nvertex): # repeat last vertex if nb<nvertex31 bnd_lon[ij,k]=bnd_lon[ij,nb-1]32 bnd_lat[ij,k]=bnd_lat[ij,nb-1]33 return bnd_lon, bnd_lat34 35 def setup_unstructured(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_own_glo, displ): # cat in 'domain','domaingroup'36 # ncell_glo (int) : global number of cells37 # index_own_glo : global indices of cells own by this MPI process38 # displ : min of index_own across all MPI processes39 bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v)40 ncell, nvertex = bnd_lon.shape41 mesh = cxios.Handle(cat,id)42 mesh.set_attr(type='unstructured')43 mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_own_glo)44 mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat)45 mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat)46 6 47 7 #----------------------------------------------------------------------------- … … 55 15 print 'xios_finalize()' 56 16 cxios.finalize() 57 17 def min(self,data): return self.comm.allreduce(data, op=MPI.MIN) 18 def max(self,data): return self.comm.allreduce(data, op=MPI.MAX) 19 58 20 class Context: 59 21 def __enter__(self): … … 84 46 ker.dynamico_xios_update_calendar(i) 85 47 def send_field_primal(self,name, data): 86 own_loc = self.mesh.primal_own_loc 48 self.send_field(self.mesh.primal_own_loc, name, data) 49 def send_field_dual(self,name, data): 50 self.send_field(self.mesh.dual_own_loc, name, data) 51 def send_field(self,own_loc, name, data): 87 52 if data.ndim==1: 88 53 data = data[own_loc] … … 93 58 cxios.send_field(name, data) 94 59 60 #----------------------------------------------------------------------------- 61 62 def setup_curvilinear(cat,id,nx,ny,own,i_index,j_index,lon,lat): 63 mesh = cxios.Handle(cat,id) 64 def log(name,data) : print 'setup_curvilinear : %s shape min max'%name, data.shape, data.min(), data.max() 65 i_index, j_index, lon, lat = [ x[own] for x in i_index, j_index, lon, lat ] 66 log('i_index',i_index) 67 log('j_index',j_index) 68 log('lon',lon) 69 log('lat',lat) 70 mesh.set_attr(type='curvilinear', ni_glo=nx, i_index=i_index, nj_glo=ny, j_index=j_index) 71 mesh.set_attr(lonvalue_1d=lon.flatten(), latvalue_1d=lat.flatten()) 72 95 73 class Context_Curvilinear(Context): 96 74 def __init__(self, mesh, nqtot, step): 97 75 self.init_llm(mesh, nqtot) 98 76 # primal mesh 99 own = mesh.primal_own_loc100 lon_i, lat_i = [ x[own] for x in mesh.lon_i, mesh.lat_i ]101 primal_i, primal_j = [ x[own] for x in mesh.primal_i, mesh.primal_j ]102 setup_curvilinear('domaingroup',' i', mesh.nx, mesh.ny, mesh.displ,103 primal_i, primal_j, lon_i, lat_i)77 setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.primal_own_loc, 78 mesh.primal_i, mesh.primal_j, mesh.lon_i, mesh.lat_i) 79 # dual mesh 80 setup_curvilinear('domaingroup','v', mesh.nx, mesh.ny, mesh.dual_own_loc, 81 mesh.dual_i, mesh.dual_j, mesh.lon_v, mesh.lat_v) 104 82 # calendar 105 83 self.init_calendar(step) 84 85 #----------------------------------------------------------------------------- 86 87 def cf_boundaries(degree,points,lon,lat): 88 n, nvertex = len(degree), degree.max() 89 bnd_lon, bnd_lat = np.zeros((n,nvertex)),np.zeros((n,nvertex)) 90 for ij in range(n): 91 nb=degree[ij] 92 for k in range(nb): 93 vertex = points[ij,k] 94 bnd_lon[ij,k]=lon[vertex] 95 bnd_lat[ij,k]=lat[vertex] 96 for k in range(nb,nvertex): # repeat last vertex if nb<nvertex 97 bnd_lon[ij,k]=bnd_lon[ij,nb-1] 98 bnd_lat[ij,k]=bnd_lat[ij,nb-1] 99 return bnd_lon, bnd_lat 100 101 def setup_unstructured(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_own_glo, displ): # cat in 'domain','domaingroup' 102 # ncell_glo (int) : global number of cells 103 # index_own_glo : global indices of cells own by this MPI process 104 # displ : min of index_own across all MPI processes 105 bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v) 106 ncell, nvertex = bnd_lon.shape 107 mesh = cxios.Handle(cat,id) 108 mesh.set_attr(type='unstructured') 109 mesh.set_attr(ni_glo=ncell_glo, ni=ncell, ibegin=displ, i_index=index_own_glo) 110 mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat) 111 mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat) 106 112 107 113 class Context_Unstructured(Context): -
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r785 r787 31 31 parser.add_argument("--Davies_N1", type=int, default=3) 32 32 parser.add_argument("--Davies_N2", type=int, default=3) 33 33 parser.add_argument("--llm", type=int, default=50) 34 parser.add_argument("--nx", type=int, default=200) 35 parser.add_argument("--ny", type=int, default=30) 36 parser.add_argument("--dt", type=float, default=360., help='Time step in seconds') 37 parser.add_argument("--hydrostatic", action='store_true') 38 parser.add_argument("--betaplane", action='store_true') 34 39 35 40 args = parser.parse_args() … … 44 49 a = 6.371229e6 # Radius of the Earth (m) 45 50 ptop = 2000. 46 y0 = Ly*0.551 y0 = .5*Ly 47 52 Cpd = 1004.5 48 53 p0 = 1e5 … … 51 56 52 57 omega = 7.292e-5 # Angular velocity of the Earth (s^-1) 53 phi0 = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg)58 phi0 = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg) 54 59 f0 = 2*omega*np.sin(phi0) 55 beta0 = 2*omega*np.cos(phi0)/a 56 fb = 2*omega*np.sin(phi0) - y0*2*omega*np.cos(phi0)/a60 beta0 = 2*omega*np.cos(phi0)/a if args.betaplane else 0. 61 fb = f0 - beta0*y0 57 62 58 63 def Phi_xy(y): … … 64 69 def ulon(x,y,eta): 65 70 u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 66 u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))71 # u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) 67 72 return u 68 73 … … 71 76 def p(eta): return p0*eta # eta = p/p_s 72 77 73 def eta(alpha) : return (1 -(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels74 def coriolis(x,y): return f0+beta0*y 78 def eta(alpha) : return (1.-(lap*ztop*alpha/T0))**(g/(Rd*lap)) # roughly equispaced levels 79 def coriolis(x,y): return f0+beta0*y # here y is 0 at domain center 75 80 76 81 filename = 'cart_%03d_%03d.nc'%(nx,ny) … … 82 87 mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 83 88 84 85 89 alpha_k = (np.arange(llm) +.5)/llm 86 90 alpha_l = (np.arange(llm+1)+ 0.)/llm 87 91 x_ik, alpha_ik = np.meshgrid(mesh.lon_i, alpha_k, indexing='ij') 88 y_ik, alpha_ik = np.meshgrid(mesh.lat_i+ .5*Ly, alpha_k, indexing='ij')92 y_ik, alpha_ik = np.meshgrid(mesh.lat_i+y0, alpha_k, indexing='ij') 89 93 x_il, alpha_il = np.meshgrid(mesh.lon_i, alpha_l, indexing='ij') 90 y_il, alpha_il = np.meshgrid(mesh.lat_i+ .5*Ly, alpha_l, indexing='ij')94 y_il, alpha_il = np.meshgrid(mesh.lat_i+y0, alpha_l, indexing='ij') 91 95 x_ek, alpha_ek = np.meshgrid(mesh.lon_e, alpha_k, indexing='ij') 92 y_ek, alpha_ek = np.meshgrid(mesh.lat_e+ .5*Ly, alpha_k, indexing='ij')96 y_ek, alpha_ek = np.meshgrid(mesh.lat_e+y0, alpha_k, indexing='ij') 93 97 94 98 print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g))) … … 103 107 Phi_il = Phi_xyeta(y_il, eta_il) 104 108 Phi_ik = Phi_xyeta(y_ik, eta_ik) 105 p_ik = p(eta_ik)109 p_ik, p_il = p(eta_ik), p(eta_il) 106 110 T_ik = T(y_ik, eta_ik) #ik full level(40), il(41) 107 111 … … 110 114 for l in range(llm): 111 115 mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l]) 116 # mass_ik[:,l]=(p_il[:,l]-p_il[:,l+1])/g 112 117 Sik, Wil = gas.s*mass_ik, mesh.field_w() 113 118 … … 224 229 225 230 g, Lx, Ly = 9.81, 4e7, 6e6 226 nx, ny, llm = 200, 30, 30231 nx, ny, llm = args.nx, args.ny, args.llm 227 232 dx,dy = Lx/nx,Ly/ny 228 233 … … 236 241 237 242 T = args.T 238 dt = 100.#180.#360.243 dt = args.dt 239 244 dz = flow0[3].max()/(params.g*llm) 240 245 nt = int(math.ceil(T/dt)) … … 242 247 logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) 243 248 249 dt=0. # FIXME 250 244 251 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 245 252 … … 252 259 else: # time stepping in Fortran 253 260 scheme = time_step.ARK2(None, dt) 254 caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 261 if args.hydrostatic: 262 caldyn_step = unst.caldyn_step_HPE(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 263 else: 264 caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) 255 265 davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 256 266 print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) … … 292 302 z=mesh.field_mass() 293 303 294 295 304 #T, nslice, dt = 3600., 1, 3600. 296 305 Nslice = args.N … … 305 314 # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W 306 315 gas, w, z = diagnose(Phi,S,m,W) 316 ps = caldyn_step.ps 307 317 curl_vk = curl(mesh, u) 318 du_fast, du_slow = caldyn_step.du_fast[0,:,:], caldyn_step.du_slow[0,:,:] 319 div_fast, div_slow = div(mesh,du_fast), div(mesh,du_slow) 308 320 # write to disk 309 context.send_field_primal('ps', davies.beta_i)321 context.send_field_primal('ps', ps) 310 322 context.send_field_primal('temp', gas.T) 311 323 context.send_field_primal('p', gas.p) … … 313 325 context.send_field_primal('uz', w) 314 326 context.send_field_primal('z', z) 315 context.send_field_primal('curl', curl_vk) 327 context.send_field_primal('div_fast', div_fast) 328 context.send_field_primal('div_slow', div_slow) 329 print curl_vk.size 330 context.send_field_dual('curl', curl_vk) 316 331 317 332 print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')) -
codes/icosagcm/devel/Python/test/xml/context_icosa_simple.xml
r622 r787 20 20 21 21 <domain_group id="i"> 22 <domain id="i" name="mesh"/> 22 <domain id="i" name="primal"/> 23 </domain_group> 24 25 <domain_group id="v"> 26 <domain id="v" name="dual"/> 23 27 </domain_group> 24 28 -
codes/icosagcm/devel/Python/test/xml/field_def_simple.xml
r785 r787 9 9 <field id="mid_ap" axis_ref="lev" long_name="hybrid A coefficient at midpoints" /> 10 10 <field id="mid_bp" axis_ref="lev" long_name="hybrid B coefficient at midpoints" /> 11 12 <field_group domain_ref="v"> 13 <field_group axis_ref="lev"> 14 <field id="curl" /> 15 </field_group> 16 </field_group> 11 17 12 18 <field_group domain_ref="i"> … … 23 29 <field id="z" /> 24 30 <field id="theta" /> 25 <field id="curl" /> 31 <field id="div_fast" /> 32 <field id="div_slow" /> 26 33 </field_group> 27 34 -
codes/icosagcm/devel/Python/test/xml/filedef_simple.xml
r785 r787 21 21 <field field_ref="theta" /> 22 22 <field field_ref="curl" /> 23 <field field_ref="div_fast" /> 24 <field field_ref="div_slow" /> 23 25 <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS" standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/> 24 26 </field_group>
Note: See TracChangeset
for help on using the changeset viewer.