Changeset 681 for codes/icosagcm


Ignore:
Timestamp:
02/14/18 00:06:52 (6 years ago)
Author:
dubos
Message:

devel/unstructured : local mesh setup + halo exchange

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

Legend:

Unmodified
Added
Removed

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

    r680 r681  
    1010import unstructured as unst 
    1111import parallel 
    12 from util import list_stencil, Base_class 
     12from util import list_stencil, Base_class, inverse_list 
    1313 
    1414radian=180/math.pi 
    … …  
    275275            if left[edge]==cell+1: ne[cell,iedge]+=1 
    276276            if right[edge]==cell+1: ne[cell,iedge]-=1 
    277             if ne[cell,iedge]==0 : print 'error at cell,iedge', cell, iedge 
     277#            if ne[cell,iedge]==0 : print 'compute_ne error at cell,iedge', cell, iedge 
    278278    return ne 
    279279 
    … …  
    299299        plt.colorbar() 
    300300        plt.ylim((-90,90)) 
    301         plt.xlim((0,360)) 
     301        plt.xlim((-180,180)) 
    302302    def plot_i(self,data, **kwargs): 
    303303        self.plot(self.primal,data, **kwargs) 
    … …  
    309309class Unstructured_Mesh(Abstract_Mesh): 
    310310    def __init__(self, gridfile, llm, nqdyn, radius, f): 
    311         # open mesh file, get main dimensions 
    312311        getvars = gridfile.getvars 
    313312        # get positions, lengths, surfaces and weights 
    … …  
    318317        primal_num, dual_num, edge_num = [ x.size for x in Ai, Av, le ] 
    319318        gridfile.primal_num, gridfile.dual_num = primal_num, dual_num 
    320  
    321         primal_deg, dual_deg, trisk_deg = getvars('primal_deg','dual_deg','trisk_deg') 
    322          
     319        primal_deg, dual_deg, trisk_deg = getvars('primal_deg','dual_deg','trisk_deg')         
    323320        # get indices for stencils  
    324321        # primal <-> edges 
    … …  
    349346                  left,right,down,up,trisk_deg,trisk, 
    350347                  self.Ai,self.Av,self.fv,le/de,Riv2,wee) 
    351         self.primal  = tri.Triangulation(lon_i*180./math.pi, lat_i*180./math.pi) 
    352         self.dual    = tri.Triangulation(lon_v*180./math.pi, lat_v*180./math.pi) 
    353         self.triedge = tri.Triangulation(lon_e*180./math.pi, lat_e*180./math.pi)         
     348        self.primal  = tri.Triangulation(lon_i*radian, lat_i*radian) 
     349        self.dual    = tri.Triangulation(lon_v*radian, lat_v*radian) 
     350        self.triedge = tri.Triangulation(lon_e*radian, lat_e*radian)         
    354351         
    355352class Unstructured_PMesh(Abstract_Mesh): # Mesh data distributed across MPI processes 
    … …  
    358355        dim_primal, dim_edge, dim_dual = gridfile.get_pdims(comm, 
    359356            'primal_num','edge_num','dual_num') 
    360         (primal_deg, primal_vertex, primal_edge, lon_i,  
    361          lat_i) = get_pvars(dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i' ) 
    362         (edge_deg, trisk_deg, left_right, down_up,  
    363          trisk) = get_pvars(dim_edge, 'edge_deg', 'trisk_deg', 'left_right', 'down_up', 'trisk') 
    364         dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v = get_pvars(dim_dual,  
    365             'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v') 
     357        primal_deg, primal_vertex, primal_edge, lon_i, lat_i, Ai = get_pvars( 
     358            dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' ) 
     359        edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars( 
     360             dim_edge, 'edge_deg', 'trisk_deg', 'left_right', 'down_up',  
     361             'trisk', 'le', 'de', 'wee', 'lon_e', 'lat_e', 'angle_e') 
     362        dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars( 
     363            dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'Riv2', 'lon_v', 'lat_v', 'Av') 
    366364        # Let indices start at 0 
    367365        for x in (primal_vertex, primal_edge, dual_vertex,dual_edge, 
    … …  
    371369        self.edge2vertex = Stencil_glob(edge_deg, down_up) 
    372370        self.vertex2edge = Stencil_glob(dual_deg, dual_edge) 
    373         self.edge2edge   = Stencil_glob(trisk_deg, trisk) 
     371        self.vertex2cell = Stencil_glob(dual_deg, dual_vertex, Riv2) 
     372        self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee) 
    374373        print 'Partitioning ...' 
    375374        edge_owner = unst.partition_mesh(trisk_deg, trisk, comm.Get_size()) 
    … …  
    377376        primal_owner = partition_from_stencil(edge_owner, primal_deg, primal_edge) 
    378377        primal_owner = parallel.LocPArray1D(dim_primal, primal_owner) 
    379         self.set_members(locals(), 'dim_primal', 'dim_dual', 'edge_owner',  
    380                          'primal_owner', 'primal_deg', 'primal_vertex',  
    381                          'lon_v', 'lat_v', 'lon_i', 'lat_i') 
    382         #        self.dim_primal, self.edge_owner, self.primal_owner = dim_primal, edge_owner, primal_owner 
     378        self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge', 
     379                         'edge_owner', 'primal_owner', 'primal_deg', 'primal_vertex',  
     380                         'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai',  
     381                         'le', 'de', 'lon_e', 'lat_e', 'angle_e') 
     382    def get(self, getter, *names): return [ getter(self.__dict__[x]) for x in names ] 
    383383 
    384384def chain(start, links): 
    … …  
    404404        yield polygon 
    405405 
    406 class Local_Mesh(Base_class): 
    407     def __init__(self, pmesh): 
    408         dim_primal, edge_owner, primal_owner = pmesh.dim_primal, pmesh.edge_owner, pmesh.primal_owner 
     406class Local_Mesh(Abstract_Mesh): 
     407    def __init__(self, pmesh, llm, nqdyn, radius, f): 
     408        dim_primal, primal_owner, dim_edge, edge_owner, dim_dual = pmesh.members( 
     409            'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual') 
    409410        print 'Constructing halos ...' 
    410411        edges_E0 = find_my_cells(edge_owner) 
    … …  
    414415        cells_C0, cells_C1 = progressive_list(cells_C0, cells_C1) 
    415416        print 'E2,E1,E0 ; C1,C0 : ', map(len, (edges_E2, edges_E1, edges_E0, cells_C1, cells_C0)) 
    416         mycells, halo_cells = cells_C0, cells_C1 
    417         get_mycells, get_halo_cells = dim_primal.getter(mycells), dim_primal.getter(halo_cells) 
    418         self.com_primal = parallel.Halo_Xchange(42, dim_primal, halo_cells, get_halo_cells(primal_owner)) 
    419     def plot(self, ax, clim, degree, bounds, lon, lat, data): 
     417        get_own_cells, get_all_cells = dim_primal.getter(cells_C0), dim_primal.getter(cells_C1) 
     418        get_all_duals, get_all_edges = dim_dual.getter(vertices_V1), dim_edge.getter(edges_E2) 
     419        self.com_primal = parallel.Halo_Xchange( 
     420            42, dim_primal, cells_C1, get_all_cells(primal_owner)) 
     421        self.com_edges = parallel.Halo_Xchange( 
     422            73, dim_edge, edges_E2, get_all_edges(edge_owner)) 
     423         
     424        self.lon_i, self.lat_i, self.Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai') 
     425        self.lon_v, self.lat_v, self.Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av') 
     426        self.fv = f(self.lon_v,self.lat_v) # Coriolis parameter 
     427        self.le, self.de, self.lon_e, self.lat_e, self.angle_e = pmesh.get( 
     428            get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 
     429        primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2)) 
     430         
     431        dict_E2, dict_C1, dict_V1 = map(inverse_list, (edges_E2, cells_C1, vertices_V1)) 
     432        down_up = pmesh.edge2vertex(edges_E2, dict_V1) 
     433        left_right = pmesh.edge2cell(edges_E2, dict_C1) 
     434        primal_edge = pmesh.cell2edge(cells_C1, dict_E2) 
     435        dual_edge = pmesh.vertex2edge(vertices_V1, dict_E2) 
     436        trisk = pmesh.edge2edge(edges_E2, dict_E2) 
     437        Riv2  = pmesh.vertex2cell(vertices_V1, dict_C1) 
     438        print 'Edge-vertex degree min/max %d %d'%(down_up.degree.min(), down_up.degree.max() ) 
     439        print 'Edge-cell degree min/max %d %d'%(left_right.degree.min(), left_right.degree.max() ) 
     440        print 'Primal degree min/max : %d %d'%(primal_edge.degree.min(), primal_edge.degree.max() ) 
     441        print 'Dual degree min/max : %d %d'%(dual_edge.degree.min(), dual_edge.degree.max() ) 
     442        print 'TRiSK degree min/max : %d %d'%(trisk.degree.min(), trisk.degree.max() ) 
     443        print 'Riv2 degree min/max : %d %d'%(Riv2.degree.min(), Riv2.degree.max() ) 
     444        primal_deg, primal_edge = primal_edge.degree, primal_edge.neigh_loc 
     445        dual_deg, dual_edge = dual_edge.degree, dual_edge.neigh_loc 
     446        trisk_deg, trisk, wee = trisk.degree, trisk.neigh_loc, trisk.weight 
     447        dual_deg, dual_vertex, Riv2 = Riv2.degree, Riv2.neigh_loc, Riv2.weight 
     448        edge_deg, left, right = left_right.degree, left_right.neigh_loc[:,0], left_right.neigh_loc[:,1] 
     449        edge_deg, down, up = down_up.degree, down_up.neigh_loc[:,0], down_up.neigh_loc[:,1] 
     450        for edge in range(edge_num):  
     451            if edge_deg[edge]<2: up[edge]=down[edge] 
     452        max_primal_deg, max_dual_deg, max_trisk_deg = [x.shape[1] for x in primal_edge, dual_edge, trisk] 
     453        # compute_ne(...) and normalize(...) expect indices starting at 1 
     454        trisk, dual_vertex, primal_edge, dual_edge =  trisk+1, dual_vertex+1, primal_edge+1, dual_edge+1 
     455        left, right, down, up = left+1, right+1, down+1, up+1 
     456        primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right) 
     457        dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down) 
     458        self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num', 
     459                         'trisk_deg', 'trisk', 'wee', 'dual_deg', 'dual_vertex', 'Riv2') 
     460 
     461        pmesh.gridfile.normalize(self,radius) 
     462        unst.init_mesh(llm, nqdyn, self.edge_num, self.primal_num, self.dual_num, 
     463                  max_trisk_deg, max_primal_deg, max_dual_deg, 
     464                  primal_deg,primal_edge,primal_ne, 
     465                  dual_deg,dual_edge,dual_ne,dual_vertex, 
     466                  left,right,down,up,trisk_deg,trisk, 
     467                  self.Ai,self.Av,self.fv,self.le/self.de,Riv2,wee) 
     468        # llm must be set before we call set_dynamico_transfer 
     469        self.com_primal.set_dynamico_transfer('primal') 
     470        self.com_edges.set_dynamico_transfer('edge') 
     471 
     472        self.primal  = tri.Triangulation(self.lon_i*radian, self.lat_i*radian) 
     473        self.dual    = tri.Triangulation(self.lon_v*radian, self.lat_v*radian) 
     474        self.triedge = tri.Triangulation(self.lon_e*radian, self.lat_e*radian)         
     475 
     476    def plot_patches(self, ax, clim, degree, bounds, lon, lat, data): 
    420477        nb_vertex = lon.size # global 
    421478        p = list(patches(degree, bounds, lon, lat)) 
    … …  
    486543    return owner1 
    487544             
    488 def find_my_cells(owner): # a PArray1D containing the data returned by partition_mesh 
     545def find_my_cells(owner): # owner : a PArray1D containing the data returned by partition_mesh 
    489546    dim, comm, owner = owner.dim, owner.dim.comm, owner.data 
    490547    mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
    … …  
    511568    for i in range(degree.size): 
    512569        for j in range(degree[i]): 
    513             bounds[i,j] = vertex_dict[bounds[i,j]]     
     570            bounds[i,j] = vertex_dict[bounds[i,j]] 
     571    return bounds 
    514572 
    515573class Stencil_glob: 
    … …  
    546604        myneigh_loc = reindex(neigh_dict, mydegree, myneigh) 
    547605        self.degree, self.neigh_glob, self.neigh_loc = mydegree, myneigh, myneigh_loc 
    548  
     606        if weight is not None: self.weight = myweight 
     607             
    549608def progressive_iter(mylist, cell_lists): 
    550609    for thelist in cell_lists:  

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

    r615 r681  
    99        order = [order[i] for i in range(len(lst))] # dict->list 
    1010        return lst,order   
    11      
     11 
     12from unstructured import ker 
     13 
    1214#----------------------- Parallel access to distributed array ------------------------------# 
    1315 
    … …  
    143145        # asssociate local index to global index 
    144146        own_dict = { glob:loc for glob,loc in zip(own_global, own_local) } 
    145         self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank]                     
    146          
     147        self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank] 
     148    def set_dynamico_transfer(self,name): 
     149        def list_to_args(lst): 
     150            ranks, lst = zip(*lst) # list of tuples => tuple of lists 
     151            lens = map(len, lst) 
     152            return len(ranks), ranks, lens, sum(lens), sum(lst,[]) 
     153        index = { 'primal':1, 'edge':2, 'dual':3 }[name] 
     154        self.index = index 
     155        send_num, send_rank, send_len, send_size, send_list = list_to_args(self.send_list)  
     156        recv_num, recv_rank, recv_len, recv_size, recv_list = list_to_args(self.recv_list)  
     157        send_rank, send_len, send_list, recv_rank, recv_len, recv_list = [ np.asarray(x, dtype=np.int32) 
     158            for x in send_rank, send_len, send_list, recv_rank, recv_len, recv_list ] 
     159        send_list, recv_list = send_list+1, recv_list+1 # Fortran expects that indices start at 1 
     160        ker.dynamico_init_transfer(index, send_num, send_size, send_rank, send_len, send_list, 
     161                                       recv_num, recv_size, recv_rank, recv_len, recv_list) 
    147162class LocalArray: # a base class for arrays with halos 
    148163    def read_own(self, parray): self.put(self.halo.own_local, self.halo.get_own(parray)) 

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

    r680 r681  
    1313        for name in names : 
    1414            self.__dict__[name]=loc[name] 
    15  
     15    def members(self, *names): 
     16        return [ self.__dict__[name] for name in names ] 

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

    r680 r681  
    6060                     ['dynamico_caldyn_unstructured', c_double, c_void_p,20], 
    6161                     ['dynamico_partition_graph', c_int,2, c_void_p,3, c_int, c_void_p], 
     62                     ['dynamico_init_transfer', c_int, c_int,2,c_void_p,3, c_int,2,c_void_p,3], 
     63                     ['dynamico_update_halo', c_int,3,c_void_p], 
    6264                     ['dynamico_morton_encode', c_int,c_void_p,4] 
    6365                     ]) 

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

    r680 r681  
     1print 'Starting' 
     2 
     3from mpi4py import MPI 
     4comm = MPI.COMM_WORLD 
     5mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() 
     6print '%d/%d starting'%(mpi_rank,mpi_size) 
     7prefix='fig_RSW2_MPAS_W02/%02d'%mpi_rank 
     8 
    19print 'Loading DYNAMICO modules ...' 
    210from dynamico import unstructured as unst 
    3 from dynamico.meshes import MPAS_Format, Unstructured_Mesh as Mesh 
     11from dynamico.meshes import MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh 
    412from dynamico import time_step 
    513print '...Done' 
    … …  
    1119print '...Done' 
    1220 
    13 grid, llm, nqdyn = 10242, 1,1 # 2562, 10242, 40962 
     21grid, llm, nqdyn = 2562, 1,1 # 2562, 10242, 40962 
    1422Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 
    15 N, T, courant = 40, 10800., 1.2 # simulation length = N*T 
     23N, T, courant = 40, 10400., 1.2 # simulation length = N*T 
    1624 
    1725print 'Omega, planetary PV', Omega, 2*Omega/gh0 
    … …  
    2028print 'Reading MPAS mesh ...' 
    2129meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid) 
    22 mesh = Mesh(meshfile, llm, nqdyn, radius, f) 
     30pmesh = PMesh(comm,meshfile) 
     31mesh = Mesh(pmesh, llm, nqdyn, radius, f) 
    2332print '...Done' 
    2433lon, lat = mesh.lon_i, mesh.lat_i 
    … …  
    3039dx = mesh.de.min() 
    3140dt = courant*dx/c0 
     41print dx, dt 
    3242nt = int(math.ceil(T/dt)) 
    3343dt = T/nt 
    … …  
    6171flow=gh,u 
    6272 
     73unst.ker.dynamico_update_halo(mesh.com_edges.index, 1, u.shape[0], u)  
     74                                   
    6375for i in range(20): 
    6476    if True: 
    … …  
    6779        print i, gh.shape, gh.min(), gh.max(), u.min(), u.max() 
    6880        mesh.plot_i(gh-Phi0) ; plt.title('err(gh)'); 
    69         plt.savefig('fig_RSW2_MPAS_W02/err_gh_%02d.png'%i) 
     81        plt.savefig('%s_err_gh_%02d.png'%(prefix,i)) 
    7082        plt.close() 
    7183    else: 
    … …  
    100112        ax8.scatter(u_nowref,u_now-u_nowref) ; ax8.set_title('u_now'); 
    101113        f.subplots_adjust(hspace=1.) 
    102         plt.savefig('fig_RSW2_MPAS_W02/scatter_%02d.png'%i) 
     114        plt.savefig('%s_scatter_%02d.png'%(prefix,i)) 
    103115        plt.close() 
    104116 

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

    r675 r681  
    1818       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, & 
    1919       caldyn_vert_cons=1, max_nb_stage=5 
    20   INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, & 
     20  INDEX,  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, & 
    2121       caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0 
    2222  LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE. 
    … …  
    5454       'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', & 
    5555       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ' /) 
     56 
     57  INTEGER, PARAMETER ::transfer_primal=1, transfer_edge=2, transfer_dual=3, transfer_max=3 
     58  TYPE Halo_transfer 
     59     INTEGER :: ranks ! size of arrays rank, len 
     60     INTEGER, ALLOCATABLE :: rank(:), & ! MPI ranks to communicate with 
     61          num(:), & ! number of cells to send to / receive from other MPI ranks 
     62          cells(:) ! local indices of cells to send/receive 
     63     DBL, ALLOCATABLE :: buf1(:), buf2(:,:), buf3(:,:,:) 
     64  END TYPE Halo_transfer 
     65  TYPE(Halo_transfer), TARGET :: send_info(transfer_max), recv_info(transfer_max) 
    5666 
    5767CONTAINS 

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

    r677 r681  
    66#endif 
    77  USE caldyn_unstructured_mod 
     8  USE transfer_unstructured_mod 
    89  IMPLICIT NONE 
    910  PRIVATE 
    … …  
    151152       DO stage=1, nb_stage 
    152153           
     154          CALL update_halo(transfer_primal, rhodz) 
     155          CALL update_halo(transfer_primal, theta_rhodz(:,:,1)) ! FIXME 
     156 
    153157          IF(hydrostatic) THEN 
    154158              
    … …  
    163167              
    164168             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u) 
    165               
     169             CALL update_halo(transfer_edge, u) 
    166170             CALL compute_pvort_only(rhodz,u,qv,qu) 
     171             CALL update_halo(transfer_edge, qu) 
     172 
    167173             CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage)) 
    168174             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage)) 
Note: See TracChangeset for help on using the changeset viewer.