Changeset 759 for codes/icosagcm


Ignore:
Timestamp:
10/11/18 18:27:15 (6 years ago)
Author:
dubos
Message:

devel/unstructured : towards XIOS output with curvilinear mesh

Location:
codes/icosagcm/devel
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/Python/dynamico/meshes.py

    r758 r759  
    182182                         'trisk_deg','trisk','Riv2','wee', 
    183183                         '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') 
    185186    def ncwrite(self, name): 
    186187        """The method writes Cartesian mesh on the disc. 
    … …  
    209210         
    210211        f.description = "Cartesian_mesh" 
    211         f.nx, f.ny = self.nx, self.ny 
     212        f.meshtype, f.nx, f.ny = 'curvilinear', self.nx, self.ny 
    212213 
    213214        def create_vars(dimname, info): 
    … …  
    218219        create_vars("primal_cell",  
    219220                    [("primal_deg","i4",self.primal_deg),  
     221                     ("primal_i","i4",self.primal_i), 
     222                     ("primal_j","i4",self.primal_j), 
    220223                     ("Ai","f8",self.Aiv), 
    221224                     ("lon_i","f8",self.lon_i), 
    222225                     ("lat_i","f8",self.lat_i)] ) 
     226 
    223227        create_vars("dual_cell",  
    224228                    [("dual_deg","i4",self.dual_deg),  
    … …  
    293297    start_index=1 # indexing starts at 1 as in Fortran 
    294298    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 
    296303        print('#NC4: Opening file:',gridfilename) 
    297304    def get_pdim(self,comm,name): return parallel.PDim(self.nc.dimensions[name], comm) 
    … …  
    310317    """ Helps class Unstructured_Mesh to read a MPAS mesh file.""" 
    311318    start_index=1 # indexing starts at 1 as in Fortran 
     319    meshtype='unstructured' 
    312320    translate= { 
    313         'primal_num':'nCells', 'edge_num':'nEdges', 'dual_num':'nVertices', 
     321        'primal_cell':'nCells', 'edge':'nEdges', 'dual_cell':'nVertices', 
    314322        'primal_deg':'nEdgesOnCell', 'trisk_deg':'nEdgesOnEdge', 
    315323        'primal_edge':'edgesOnCell', 'dual_edge':'edgesOnVertex',  
    … …  
    452460         
    453461class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes 
    454     def __init__(self,comm, gridfile, meshtype='unstructured'): 
     462    def __init__(self,comm, gridfile): 
    455463        self.comm, self.gridfile, get_pvars = comm, gridfile, gridfile.get_pvars 
    456464        dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, 
    … …  
    463471        dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars( 
    464472            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 
    465476        # Let indices start at 0 
    466477        for x in (primal_vertex, primal_edge, dual_vertex,dual_edge, 
    … …  
    478489        primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 
    479490        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', 
    481492                         'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex',  
    482493                         'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai',  
    … …  
    530541        primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2)) 
    531542 
    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 
    534548 
    535549        # construct local stencils from global stencils 
    … …  
    566580        primal_own_all = comm.allgather(primal_own_num) 
    567581        displ = sum(primal_own_all[0:comm.Get_rank()]) # required by XIOS 
     582 
    568583        print 'local primal mesh :', primal_vertex.min(), primal_vertex.max(), primal_own_all, displ 
    569584         
    … …  
    742757    # example : edge_list = progressive_list(E0,E1,E2) with E0,E1,E2 increasing lists of cell edges 
    743758    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/dynamico/xios.py

    r696 r759  
    44from dynamico.meshes import radian 
    55from cxios import finalize 
     6 
     7#----------------------------------------------------------------------------- 
     8 
     9def 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()) 
    617 
    718def cf_boundaries(degree,points,lon,lat): 
    … …  
    1930    return bnd_lon, bnd_lat 
    2031 
    21 def setup_domain(cat,id, degree,vertex,lon,lat,lon_v,lat_v, ncell_glo, index_glo, displ): # cat in 'domain','domaingroup' 
     32def 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 
    2236    bnd_lon, bnd_lat = cf_boundaries(degree, vertex, lon_v, lat_v) 
    2337    ncell, nvertex = bnd_lon.shape 
    2438    mesh = cxios.Handle(cat,id) 
    2539    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) 
    2741    mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat) 
    2842    mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat) 
    2943 
     44#----------------------------------------------------------------------------- 
     45 
    3046class XIOS_Context: 
    31     def __init__(self, pmesh, mesh,nqtot, step): 
     47    def init_llm(self, mesh, nqtot): 
    3248        self.mesh=mesh 
    3349        ker.dynamico_setup_xios() 
    … …  
    3854        levp1.set_attr(n_glo=llm+1, value=np.arange(llm+1)) 
    3955        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): 
    4757        calendar=cxios.Handle('calendar_wrapper') 
    4858        dtime = cxios.Duration(second=step) 
    … …  
    5363        freq_offset = cxios.Duration(second=0) 
    5464        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 
     70class 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 
     92class 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) 
    55103        print 'cxios.context_close_definition()' 
    56104        cxios.context_close_definition() 
    … …  
    64112        data = np.ascontiguousarray(data) 
    65113        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  
        11build 
        2 xios.c 
         2cxios.c 
        33unstructured.c 

  • codes/icosagcm/devel/Python/src/cxios.pyx

    r694 r759  
    6767            fun = cxios.vardict[prefix+key] 
    6868            extra=value 
    69             if type(value) in (int,c_int): 
     69            if type(value) in (int,c_int,np.int64,np.int32): 
    7070                fun(handle,c_int(value)) 
    7171            elif type(value) is np.ndarray: 
    … …  
    102102                       ['update_calendar',c_int] ]) 
    103103 
    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']) 
    106107    import_set_attr(['field','fieldgroup'], [Duration,'freq_op','freq_offset']) 
    107108    import_set_attr(['axis'], [c_int,'n_glo'], [np.ndarray,'value']) 

  • codes/icosagcm/devel/Python/test/py/RSW_2D_mesh.py

    r756 r759  
    1515dx,dy=Lx/nx,Ly/ny 
    1616 
    17 filename, llm, nqdyn, g, f, radius = 'cart_128.nc', 1, 1, 1., 1., None 
     17filename, llm, nqdyn, g, f, radius = 'cart_%03d_%03d.nc'%(nx,ny), 1, 1, 1., 1., None 
    1818unst.setvar('g',g) 
    1919 
    … …  
    2727caldyn = unst.Caldyn_RSW(mesh) 
    2828 
    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 
     31xx,yy = mesh.lat_i, mesh.lon_i 
    3132 
    3233x1,x2,yy = xx-1., xx+1., yy 

  • codes/icosagcm/devel/Python/test/py/test_xios.py

    r696 r759  
    2828#--------------------------------- write some data ---------------------------------------- 
    2929 
    30 context=xios.XIOS_Context(pmesh,mesh,nqtot, 3600) 
     30context=xios.XIOS_Context_Unstructured(pmesh,mesh,nqtot, 3600) 
    3131 
    3232lat_i = radian*mesh.lat_i 

  • codes/icosagcm/devel/Python/test/python.sh

    r717 r759  
    88{ 
    99    LD_PRELOAD=$PYTHON_PRELOAD python -u $* 
     10} 
     11 
     12function cmd_gdb() 
     13{ 
     14    set -x 
     15    LD_PRELOAD=$PYTHON_PRELOAD gdb --args python -u $* 
    1016} 
    1117 

  • codes/icosagcm/devel/Python/test/xml/filedef_simple.xml

    r755 r759  
    4444    <variable name="time_frequency" type="string" > 24h </variable>     
    4545 
    46   </file> 
    47 --> 
     46  </file> --> 
    4847 
    4948</file_definition> 

  • codes/icosagcm/devel/Python/test/xml/link.sh

    r755 r759  
    77        ln -s xml/${name}_$1.xml ${name}.xml 
    88    done 
    9     ln -s iodef.xml 
     9    ln -s xml/iodef.xml 
    1010} 
    1111 

  • codes/icosagcm/devel/src/unstructured/data_unstructured.F90

    r749 r759  
    164164    dual_ne(:,:) = dual_ne_(:,:) 
    165165    dual_vertex(:,:) = dual_vertex_(:,:) 
    166     IF(MINVAL(dual_deg)<3) THEN 
    167        STOP 'At least one dual cell has less than 3 vertices' 
    168     END IF 
    169     IF(MINVAL(primal_deg)<3) THEN 
    170        STOP 'At least one primal cell has less than 3 vertices' 
     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' 
    171171    END IF 
    172172    left(:)=left_(:) 
Note: See TracChangeset for help on using the changeset viewer.