Changeset 760


Ignore:
Timestamp:
10/15/18 01:04:44 (6 years ago)
Author:
dubos
Message:

devel/Python : block-wise partitioning

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

Legend:

Unmodified
Added
Removed

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

    r759 r760  
    167167 
    168168        primal_i, primal_j = [ x.astype(np.int32) for x in lon_i, lat_i] 
     169        edge_i, edge_j = [ x.astype(np.int32) for x in lon_e, lat_e] 
    169170        lon_i, lon_v, lon_e = [x*dx-Lx/2 for x in lon_i, lon_v, lon_e] 
    170171        lat_i, lat_v, lat_e = [y*dy-Ly/2 for y in lat_i, lat_v, lat_e] 
    … …  
    183184                         'Aiv','lon_i','lat_i','lon_v','lat_v', 
    184185                         'le_de','le','de','lon_e','lat_e','angle_e', 
    185                          'primal_i','primal_j') 
     186                         'primal_i','primal_j','edge_i','edge_j') 
    186187    def ncwrite(self, name): 
    187188        """The method writes Cartesian mesh on the disc. 
    … …  
    249250                     ("lon_e","f8", self.lon_e), 
    250251                     ("lat_e","f8", self.lat_e), 
    251                       ("angle_e","f8", self.angle_e)] ) 
     252                     ("angle_e","f8", self.angle_e), 
     253                     ("edge_i","i4",self.edge_i), 
     254                     ("edge_j","i4",self.edge_j)] ) 
    252255 
    253256        create_vars( ("edge","TWO"), 
    … …  
    473476        if gridfile.meshtype == 'curvilinear': 
    474477            self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 
     478            self.edge_i, self.edge_j = get_pvars(dim_edge, 'edge_i','edge_j') 
    475479 
    476480        # Let indices start at 0 
    … …  
    484488        self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2) 
    485489        self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee) 
    486         print 'Partitioning ...' 
    487         edge_owner = unst.partition_mesh(trisk_deg, trisk, comm.Get_size()) 
    488         edge_owner = parallel.LocPArray1D(dim_edge, edge_owner) 
    489         primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 
    490         primal_owner = parallel.LocPArray1D(dim_primal, primal_owner) 
    491490        self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge', 
    492                          'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex',  
     491                         'primal_deg', 'primal_vertex', 'primal_edge', 'trisk_deg', 'trisk', 
    493492                         'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai',  
    494493                         'le', 'de', 'lon_e', 'lat_e', 'angle_e') 
    495 #        if meshtype == 'curvilinear' : # read extra information from mesh file 
    496 #            self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 
    497              
     494    def partition(self, edge_owner): 
     495        self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner) 
     496        primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge) 
     497        self.primal_owner = parallel.LocPArray1D(self.dim_primal, primal_owner) 
     498 
     499    def partition_metis(self): 
     500        print 'Partitioning with ParMetis...' 
     501        self.partition(unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())) 
     502 
     503    def partition_curvilinear(self, ni,nj): 
     504        def owner(dim,i,j): 
     505            owner_i, owner_j =  (ni*i.data)/nx, (nj*j.data)/ny 
     506            return parallel.LocPArray1D(dim, owner_i + ni*owner_j) 
     507        nx, ny = self.gridfile.nx, self.gridfile.ny 
     508        print 'Partitioning %d x %d cells in %d x %d blocks ...'%(nx,ny,ni,nj) 
     509        n = self.comm.Get_size() 
     510        assert n == ni*nj, 'Mismatch between ni x nj = %d and MPI_SIZE = %d.'%(ni*nj, n) 
     511        self.edge_owner   = owner(self.dim_edge,   self.edge_i,   self.edge_j) 
     512        self.primal_owner = owner(self.dim_primal, self.primal_i, self.primal_j) 
     513        print 'partition_curvilinear %d :'%(self.comm.Get_rank()), self.primal_owner.data 
     514 
    498515    def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ] 
    499516 

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

    r759 r760  
    33from unstructured import ker 
    44from dynamico.meshes import radian 
    5 from cxios import finalize 
     5from mpi4py import MPI 
    66 
    77#----------------------------------------------------------------------------- 
    … …  
    99def setup_curvilinear(cat,id,nx,ny,displ,index_own_i,index_own_j,lon,lat): 
    1010    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 
     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) 
    1316    mesh.set_attr(type='curvilinear', 
    1417                  ni_glo=nx, i_index=index_own_i, 
    … …  
    4447#----------------------------------------------------------------------------- 
    4548 
    46 class XIOS_Context: 
     49class Client: 
     50    def __enter__(self): 
     51        ker.dynamico_setup_xios() 
     52        self.comm=MPI.COMM_WORLD # FIXME : we should use the communicator returned by XIOS 
     53        return self 
     54    def __exit__(self, type, value, traceback): 
     55        print 'xios_finalize()' 
     56        cxios.finalize() 
     57 
     58class Context: 
     59    def __enter__(self): 
     60        print 'cxios.context_close_definition()' 
     61        cxios.context_close_definition() 
     62        return self 
     63    def __exit__(self, type, value, traceback): 
     64        print 'xios_context_finalize()' 
     65        cxios.context_finalize() 
    4766    def init_llm(self, mesh, nqtot): 
    4867        self.mesh=mesh 
    49         ker.dynamico_setup_xios() 
    5068        # vertical axis, nq 
    5169        llm = mesh.llm 
    … …  
    6583    def update_calendar(self, i): 
    6684        ker.dynamico_xios_update_calendar(i) 
    67     def finalize(self): 
    68         cxios.context_finalize() 
     85    def send_field_primal(self,name, data): 
     86        own_loc = self.mesh.primal_own_loc 
     87        if data.ndim==1: 
     88            data = data[own_loc] 
     89        if data.ndim==2: 
     90            data = data[own_loc,:] 
     91            data = data.transpose() # XIOS expects contiguous horizontal slices 
     92        data = np.ascontiguousarray(data) 
     93        cxios.send_field(name, data) 
    6994 
    70 class XIOS_Context_Curvilinear(XIOS_Context): 
     95class Context_Curvilinear(Context): 
    7196    def __init__(self, mesh, nqtot, step): 
    7297        self.init_llm(mesh, nqtot) 
    7398        # primal mesh 
    74 #        lat_i,lon_i = np.meshgrid(mesh.y,mesh.x, indexing='ij') 
    75 #        ny,nx=lat_i.shape 
    7699        own = mesh.primal_own_loc 
    77100        lon_i, lat_i = [ x[own]*radian for x in mesh.lon_i, mesh.lat_i ] 
    78101        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 
    80102        setup_curvilinear('domaingroup','i', mesh.nx, mesh.ny, mesh.displ, 
    81103                          primal_i, primal_j, lon_i, lat_i) 
    82104        # calendar 
    83105        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) 
    91106 
    92 class XIOS_Context_Unstructured(XIOS_Context): 
     107class Context_Unstructured(Context): 
    93108    def __init__(self, pmesh, mesh,nqtot, step): 
    94109        self.init_llm(mesh, nqtot) 
    … …  
    101116        # calendar 
    102117        self.init_calendar(step) 
    103         print 'cxios.context_close_definition()' 
    104         cxios.context_close_definition() 
    105     def send_field_primal(self,name, data): 
    106         own_loc = self.mesh.primal_own_loc 
    107         if data.ndim==1: 
    108             data = data[own_loc] 
    109         if data.ndim==2: 
    110             data = data[own_loc,:] 
    111             data = data.transpose() # XIOS expects contiguous horizontal slices 
    112         data = np.ascontiguousarray(data) 
    113         cxios.send_field(name, data) 

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

    r759 r760  
    2828    cdef char* idptr = id 
    2929    cdef int idlen = len(id) 
    30     print 'cxios.write_data', id, data.shape 
     30#    print 'cxios.write_data', id, data.shape 
    3131    ndim=data.ndim 
    3232    if ndim==1: send_field1(idptr,idlen,data) 
    … …  
    6161            create_fun(byref(self.handle), id, c_int(len(id))) 
    6262        self.cat, self.prefix_set, self.id = cat, 'set_'+cat+'_', id 
    63         print 'Handle', cat, self.id, self.handle 
     63#        print 'Handle', cat, self.id, self.handle 
    6464    def set_attr(self, **kwargs): 
    6565        prefix, handle = self.prefix_set, self.handle 

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

    r701 r760  
    5858#    mesh = meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, f) 
    5959    pmesh = meshes.Unstructured_PMesh(comm,meshfile) 
     60    pmesh.partition_metis() 
    6061    mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, f) 
    6162    mesh.pmesh=pmesh 

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

    r749 r760  
    2929meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 
    3030pmesh = PMesh(comm,meshfile) 
     31pmesh.partition_metis() 
    3132mesh = Mesh(pmesh, llm, nqdyn, radius, f) 
    3233print '...Done' 

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

    r692 r760  
    8080meshfile = meshes.MPAS_Format('grids/%s.grid.nc'%grid) 
    8181pmesh = meshes.Unstructured_PMesh(comm, meshfile) 
     82pmesh.partition_metis() 
    8283 
    8384def coriolis(lon,lat): return 0.*lat 

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

    r759 r760  
    1 import dynamico.unstructured as unst 
    2 import dynamico.xios as xios 
     1from dynamico import unstructured as unst 
     2from dynamico import xios 
    33from dynamico.meshes import radian, MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh 
     4import numpy as np 
    45 
    5 from mpi4py import MPI 
    6 comm = MPI.COMM_WORLD 
    7 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
    8 print '%d/%d starting'%(mpi_rank,mpi_size) 
    9  
    10 import numpy as np 
    11 import cProfile 
    12 import re 
     6with xios.Client() as client:     
     7    comm = client.comm 
     8    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
     9    print '%d/%d starting'%(mpi_rank,mpi_size) 
    1310 
    1411#----------------------------- read MPAS mesh -------------------------------- 
    1512 
    16 grid, llm, nqdyn, nqtot = 2562, 4,1,1 # 2562, 10242, 40962   
    17 Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 
    18 N, T, courant = 40, 10800., 1.2 # simulation length = N*T 
    19 print 'Omega, planetary PV', Omega, 2*Omega/gh0 
     13    grid, llm, nqdyn, nqtot = 2562, 4,1,1 # 2562, 10242, 40962   
     14    Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 
     15    N, T, courant = 40, 10800., 1.2 # simulation length = N*T 
     16    print 'Omega, planetary PV', Omega, 2*Omega/gh0 
    2017 
    21 def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter                                                                                                                 
    22 print 'Reading MPAS mesh ...' 
    23 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 
    24 pmesh = PMesh(comm,meshfile) 
    25 mesh = Mesh(pmesh, llm, nqdyn, radius, f) 
    26 print '...Done' 
     18    def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter                                                                                                                 
     19    print 'Reading MPAS mesh ...' 
     20    meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 
     21    pmesh = PMesh(comm,meshfile) 
     22    pmesh.partition_metis() 
     23    mesh = Mesh(pmesh, llm, nqdyn, radius, f) 
     24    print '...Done' 
    2725 
    2826#--------------------------------- write some data ---------------------------------------- 
    2927 
    30 context=xios.XIOS_Context_Unstructured(pmesh,mesh,nqtot, 3600) 
     28    with xios.Context_Unstructured(pmesh,mesh,nqtot, 3600) as context: 
     29        lat_i = radian*mesh.lat_i 
     30        lat_ik, junk  = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 
    3131 
    32 lat_i = radian*mesh.lat_i 
    33 lat_ik, junk  = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 
    34  
    35 no_error=True 
    36 try: 
    37     for i in range(100): 
    38         context.update_calendar(i) 
    39         print 'send_field', i, lat_ik.shape 
    40         context.send_field_primal('ps', lat_i) 
    41         context.send_field_primal('theta', lat_ik) 
    42 except: 
    43     print 'An error has occured, trying to close XIOS properly' 
    44     no_error=False 
    45  
    46 print 'xios.context_finalize()' 
    47 context.finalize() 
    48 print 'xios.finalize()' 
    49 xios.finalize() 
    50 print 'Done' 
    51  
    52 if not no_error : raise 
     32        for i in range(100): 
     33            context.update_calendar(i) 
     34            print 'send_field', i, lat_ik.shape 
     35            context.send_field_primal('ps', lat_i) 
     36            context.send_field_primal('theta', lat_ik) 

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

    r759 r760  
    33from dynamico import unstructured as unst 
    44from dynamico import meshes 
    5 ker=unst.ker 
    6  
    7 from mpi4py import MPI 
    8 comm = MPI.COMM_WORLD 
    9 mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
    10 print '%d/%d starting'%(mpi_rank,mpi_size) 
    11  
    125import numpy as np 
    136 
    14 nqdyn,nqtot,llm = 1,1,4 
     7# client.__exit__() and context.__exit__() are called at the end of the 'with' blocks 
     8# even if an exception is raised in the body of the block 
    159 
    16 if False: 
    17     nx,ny,llm = 64,64 
    18     Lx,Ly,f = 1.,1.,1. 
    19     mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) 
    20 else: 
    21     filename, g, f, radius = 'cart_128_128.nc', 1., 1., None 
    22     filename = 'cart_004_004.nc' 
    23     unst.setvar('g',g) 
    24     print 'Reading Cartesian mesh ...' 
    25     def coriolis(lon,lat): 
    26         return f+0.*lon 
    27     meshfile = meshes.DYNAMICO_Format(filename) 
    28     nx,ny=meshfile.nx,meshfile.ny 
     10with xios.Client() as client: 
     11    comm = client.comm 
     12    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
     13    print '%d/%d starting'%(mpi_rank,mpi_size) 
    2914 
    30     pmesh = meshes.Unstructured_PMesh(comm,meshfile) 
    31     mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 
     15    nqdyn,nqtot,llm = 1,1,4 
     16 
     17    if False: 
     18        nx,ny,llm = 64,64 
     19        Lx,Ly,f = 1.,1.,1. 
     20        mesh = meshes.Cartesian_Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) 
     21    else: 
     22        filename, g, f, radius = 'cart_128_128.nc', 1., 1., None 
     23    #    filename = 'cart_008_008.nc' 
     24        unst.setvar('g',g) 
     25        print 'Reading Cartesian mesh ...' 
     26        def coriolis(lon,lat): 
     27            return f+0.*lon 
     28        meshfile = meshes.DYNAMICO_Format(filename) 
     29        nx,ny=meshfile.nx,meshfile.ny 
     30 
     31        pmesh = meshes.Unstructured_PMesh(comm,meshfile) 
     32    #    pmesh.partition_metis() 
     33        pmesh.partition_curvilinear(2,2) 
     34        mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) 
    3235 
    3336#--------------------------------- write some data ---------------------------------------- 
    3437 
    35 context=xios.XIOS_Context_Curvilinear(mesh,nqtot, 3600) 
    36 #lon_i,lat_i = mesh.xx[:,:,0].flatten(), mesh.yy[:,:,0].flatten() 
    37 lon_i,lat_i = mesh.lon_i, mesh.lat_i 
    38 lat_ik, junk  = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 
     38    with xios.Context_Curvilinear(mesh,nqtot, 3600) as context: 
     39        lon_i,lat_i = mesh.lon_i, mesh.lat_i 
     40        lat_ik, junk  = np.meshgrid(lat_i, np.arange(llm), indexing='ij') 
     41        lon_ik, junk  = np.meshgrid(lon_i, np.arange(llm), indexing='ij') 
    3942 
    40 no_error=True 
    41 try: 
    42     for i in range(3): 
    43         context.update_calendar(i) 
    44         print 'send_field', i, lat_ik.shape 
    45         context.send_field_primal('ps', lat_i) 
    46         context.send_field_primal('theta', lat_ik) 
    47 except: 
    48     print 'An error has occured, trying to close XIOS properly' 
    49     no_error=False 
    50  
    51 print 'xios.context_finalize()' 
    52 context.finalize() 
    53 print 'xios.finalize()' 
    54 xios.finalize() 
    55 print 'Done' 
    56  
    57 if not no_error : raise 
     43        for i in range(48): 
     44            context.update_calendar(i) 
     45#            print 'send_field', i, lat_ik.shape 
     46            context.send_field_primal('ps', lat_i+lon_i) 
     47            context.send_field_primal('theta', lat_ik+lon_ik) 

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

    r755 r760  
    33        <context  id="xios"> 
    44                <variable_definition> 
    5                 <variable id="print_file" type="bool"> false </variable>  
     5                <variable id="print_file" type="bool"> true </variable>  
    66                <variable_group id="buffer"> 
    77                </variable_group> 
Note: See TracChangeset for help on using the changeset viewer.