Changeset 787


Ignore:
Timestamp:
11/23/18 13:06:47 (5 years ago)
Author:
dubos
Message:

devel/Python : XIOS output on dual mesh

Location:
codes/icosagcm/devel/Python
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/devel/Python/dynamico/meshes.py

    r779 r787  
    478478        if gridfile.meshtype == 'curvilinear': 
    479479            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') 
    480481            self.edge_i, self.edge_j = get_pvars(dim_edge, 'edge_i','edge_j') 
    481482 
     
    494495                         'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai',  
    495496                         '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()) 
    497500        self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner) 
    498501        primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge) 
    499502        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) 
    504505 
    505506    def partition_curvilinear(self, ni,nj): 
     
    513514        self.edge_owner   = owner(self.dim_edge,   self.edge_i,   self.edge_j) 
    514515        self.primal_owner = owner(self.dim_primal, self.primal_i, self.primal_j) 
     516        self.dual_owner = self.primal_owner 
    515517        print 'partition_curvilinear %d :'%(self.comm.Get_rank()), self.primal_owner.data 
    516518 
     
    541543class Local_Mesh(Abstract_Mesh): 
    542544    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' ) 
    545547        print 'Constructing halos ...' 
    546548        edges_E0 = find_my_cells(edge_owner) 
     
    564566            self.nx, self.ny = pmesh.gridfile.nx, pmesh.gridfile.ny 
    565567            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') 
    566569 
    567570 
     
    599602        primal_own_all = comm.allgather(primal_own_num) 
    600603        displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS 
    601  
    602604        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] 
    603610         
    604611        # compute_ne(...) and normalize(...) expect indices starting at 1 
     
    610617        self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num', 
    611618                         'primal_own_glo', 'primal_own_loc', 'primal_own_deg', 
    612                          'primal_deg', 'primal_vertex', 'displ', 
     619                         'primal_deg', 'primal_vertex', 'displ', 'dual_own_loc',  
    613620                         'trisk_deg', 'trisk', 'wee',  
    614621                         'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2', 
  • codes/icosagcm/devel/Python/dynamico/xios.py

    r783 r787  
    44from dynamico.meshes import radian 
    55from 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<nvertex 
    31             bnd_lon[ij,k]=bnd_lon[ij,nb-1] 
    32             bnd_lat[ij,k]=bnd_lat[ij,nb-1] 
    33     return bnd_lon, bnd_lat 
    34  
    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 cells 
    37     # index_own_glo : global indices of cells own by this MPI process 
    38     # displ : min of index_own across all MPI processes 
    39     bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v) 
    40     ncell, nvertex = bnd_lon.shape 
    41     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) 
    466 
    477#----------------------------------------------------------------------------- 
     
    5515        print 'xios_finalize()' 
    5616        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     
    5820class Context: 
    5921    def __enter__(self): 
     
    8446        ker.dynamico_xios_update_calendar(i) 
    8547    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): 
    8752        if data.ndim==1: 
    8853            data = data[own_loc] 
     
    9358        cxios.send_field(name, data) 
    9459 
     60#----------------------------------------------------------------------------- 
     61 
     62def 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 
    9573class Context_Curvilinear(Context): 
    9674    def __init__(self, mesh, nqtot, step): 
    9775        self.init_llm(mesh, nqtot) 
    9876        # primal mesh 
    99         own = mesh.primal_own_loc 
    100         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) 
    10482        # calendar 
    10583        self.init_calendar(step) 
     84 
     85#----------------------------------------------------------------------------- 
     86 
     87def 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 
     101def 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) 
    106112 
    107113class Context_Unstructured(Context): 
  • codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py

    r785 r787  
    3131parser.add_argument("--Davies_N1", type=int, default=3) 
    3232parser.add_argument("--Davies_N2", type=int, default=3) 
    33  
     33parser.add_argument("--llm", type=int, default=50) 
     34parser.add_argument("--nx", type=int, default=200) 
     35parser.add_argument("--ny", type=int, default=30) 
     36parser.add_argument("--dt", type=float, default=360., help='Time step in seconds') 
     37parser.add_argument("--hydrostatic", action='store_true') 
     38parser.add_argument("--betaplane", action='store_true') 
    3439 
    3540args = parser.parse_args() 
     
    4449    a     = 6.371229e6 # Radius of the Earth (m) 
    4550    ptop  = 2000. 
    46     y0    = Ly*0.5 
     51    y0    = .5*Ly 
    4752    Cpd   = 1004.5 
    4853    p0    = 1e5 
     
    5156 
    5257    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) 
    5459    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)/a 
     60    beta0  = 2*omega*np.cos(phi0)/a if args.betaplane else 0. 
     61    fb     = f0 - beta0*y0 
    5762 
    5863    def Phi_xy(y): 
     
    6469    def ulon(x,y,eta): 
    6570        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)) 
    6772        return  u 
    6873 
     
    7176    def p(eta): return p0*eta     # eta  = p/p_s 
    7277 
    73     def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels 
    74     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 
    7580 
    7681    filename = 'cart_%03d_%03d.nc'%(nx,ny) 
     
    8287    mesh  = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 
    8388 
    84  
    8589    alpha_k = (np.arange(llm)  +.5)/llm 
    8690    alpha_l = (np.arange(llm+1)+ 0.)/llm 
    8791    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') 
    8993    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') 
    9195    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') 
    9397 
    9498    print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g)))  
     
    103107    Phi_il = Phi_xyeta(y_il, eta_il) 
    104108    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) 
    106110    T_ik = T(y_ik, eta_ik)  #ik full level(40), il(41) 
    107111 
     
    110114    for l in range(llm):  
    111115        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 
    112117    Sik, Wil  = gas.s*mass_ik, mesh.field_w() 
    113118     
     
    224229 
    225230    g, Lx, Ly = 9.81, 4e7, 6e6 
    226     nx, ny, llm = 200, 30, 30 
     231    nx, ny, llm = args.nx, args.ny, args.llm 
    227232    dx,dy = Lx/nx,Ly/ny 
    228233 
     
    236241 
    237242    T  = args.T 
    238     dt = 100.#180.#360. 
     243    dt = args.dt 
    239244    dz = flow0[3].max()/(params.g*llm) 
    240245    nt = int(math.ceil(T/dt)) 
     
    242247    logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) 
    243248     
     249    dt=0. # FIXME 
     250 
    244251    caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass 
    245252 
     
    252259    else: # time stepping in Fortran 
    253260        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) 
    255265        davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) 
    256266        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) 
     
    292302    z=mesh.field_mass() 
    293303 
    294  
    295304    #T, nslice, dt = 3600., 1, 3600. 
    296305    Nslice = args.N 
     
    305314            # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W 
    306315            gas, w, z = diagnose(Phi,S,m,W) 
     316            ps = caldyn_step.ps 
    307317            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) 
    308320            # write to disk 
    309             context.send_field_primal('ps', davies.beta_i) 
     321            context.send_field_primal('ps', ps) 
    310322            context.send_field_primal('temp', gas.T) 
    311323            context.send_field_primal('p', gas.p) 
     
    313325            context.send_field_primal('uz', w) 
    314326            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) 
    316331 
    317332            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  
    2020 
    2121     <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"/> 
    2327     </domain_group> 
    2428 
  • codes/icosagcm/devel/Python/test/xml/field_def_simple.xml

    r785 r787  
    99    <field id="mid_ap" axis_ref="lev"   long_name="hybrid A coefficient at midpoints" /> 
    1010    <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> 
    1117     
    1218    <field_group domain_ref="i"> 
     
    2329        <field id="z" /> 
    2430        <field id="theta" /> 
    25         <field id="curl" /> 
     31        <field id="div_fast" /> 
     32        <field id="div_slow" /> 
    2633      </field_group> 
    2734      
  • codes/icosagcm/devel/Python/test/xml/filedef_simple.xml

    r785 r787  
    2121      <field field_ref="theta" /> 
    2222      <field field_ref="curl" /> 
     23      <field field_ref="div_fast" /> 
     24      <field field_ref="div_slow" /> 
    2325      <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS"   standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/>       
    2426    </field_group> 
Note: See TracChangeset for help on using the changeset viewer.