MigrateMesh.F90

Go to the documentation of this file.

 
program migratemesh
  implicit none
 
#include "moab/MOABConfig.h"
#ifdef MOAB_HAVE_MPI
#  include "mpif.h"
#else
#  error "enable parallel build"
#endif
 
!#define NONOVERLAP
 
    ! init the parallel partition
    integer ierr, sz, rank, i
    integer  newcomm
    integer gcomm, comm1, comm2
    integer pid1, pid2 !  this is for physics ids
    integer compid1, compid2  !  component ids are unique over all pes, and established in
                              !  advance;
    integer nghlay ! number of ghost layers for loading
    integer grouptasks(9)   !   run on at most 9 processes
    integer startg1, startg2, endg1, endg2    !   start and end for group tasks, for creation
    integer sizeg1, sizeg2           !   size of the group that gets created
    character*10 appname
    character*132 readopts
    character*132 filename
    character*132 outfile
    character*132 wopts
    integer allgroup, group1, group2 ! Corresponding to MPI_Group in C
    integer tagcomm1, tagcomm2
    integer imoab_initializefortran, imoab_registerfortranapplication
    integer imoab_loadmesh, imoab_sendmesh, imoab_receivemesh, imoab_writemesh
    integer imoab_freesenderbuffers
    integer imoab_deregisterapplication, imoab_finalize
    integer repart_scheme , context_id
 
    call mpi_init(ierr)
    call mpi_comm_dup(mpi_comm_world, gcomm, ierr)
    call mpi_comm_size(gcomm, sz, ierr)
    call mpi_comm_rank(gcomm, rank, ierr)
    if (rank .eq. 0) print *, "size:", sz
    call errorout(ierr, 'cannot get rank' )
    if ( (0 .eq. rank) .and.  (sz>9)  ) then
      print *, "size is " , sz, ". run on at most 9 tasks "
      call exit(1)
    endif
    ! create 2 overlapping groups, for generality
    ! create communicators for each group; 
    ! one group will represent the sender, the other group the receiver
    ! about one third of tasks will be on group 1 only,  and one fourth will be on group 2 only
    !  about (1-1./3 -1./4) will be overlapping, these tasks will be common to both groups
    ! the mesh will be read on the sender comm, will be sent to receiver comm 
 
    ! create new MPI groups for processors 1/3*sz, 1/3*sz+1, ..., sz-1 (group 1) and 0, 1, .., 3/4*sz-1 (group 2)
    
 
    call mpi_comm_group (gcomm, allgroup, ierr)
    call errorout(ierr, 'cannot get world group' )
    ! first group, sz/3 to sz-1
    startg1 = sz/2
    endg1 = sz-1
    sizeg1 = endg1 - startg1 + 1
    ! used for new API in iMOAB, for tag migrate, release buffers
    context_id = -1
    
    do i=1, sizeg1
      grouptasks(i) = startg1+i-1
    end do 
 
    call mpi_group_incl(allgroup, sizeg1, grouptasks, group1, ierr)
    call errorout(ierr, 'cannot create group 1' )
 
    ! second group, 0, 1, 3/4*sz
    startg2 = 0
    endg2 = 3*sz/4 -1
    if (endg2 <0) endg2 = 0  !  so it will work even for 1 task
#ifdef NONOVERLAP
    endg2 = startg1-1
#endif
    sizeg2 = endg2 - startg2 + 1
    do i=1, sizeg2
      grouptasks(i) = startg2+i-1
    enddo 
    
    call mpi_group_incl(allgroup, sizeg2, grouptasks, group2, ierr)
    call errorout(ierr, 'cannot create group 2' )
 
    if ( (0 .eq. rank) ) then
      print *, "group 1 tasks: ", (i, i=startg1, endg1)
      print *, "group 2 tasks: ", (i, i=startg2, endg2)
    endif
    ! now create both communicators
    !  when we are not on tasks in the communicator, the MPI_Comm created will be null
    tagcomm1 = 1
    call mpi_comm_create_group(gcomm, group1, tagcomm1, comm1, ierr)
    call errorout(ierr, 'cannot create communicator 1' )
 
    tagcomm2 = 2
    call mpi_comm_create_group(gcomm, group2, tagcomm2, comm2, ierr)
    call errorout(ierr, 'cannot create communicator 2' )
 
 
    ierr = imoab_initializefortran()
 
    repart_scheme = 0 !  this is for trivial partitioning
#ifdef MOAB_HAVE_ZOLTAN
    repart_scheme = 1 !  use the graph partitioner in that case
#endif
    ! give some dummy values to component ids, just to differentiate between them
    ! the par comm graph is unique between components
    compid1 = 4
    compid2 = 7
    call errorout(ierr, 'did not initialize fortran' )
    if (rank == 0) print *, "initialize iMOAB fortran applications"
 
    if (comm1 /= mpi_comm_null) then
       appname='phis1'//char(0)
       ierr = imoab_registerfortranapplication(trim(appname), comm1, compid1, pid1)
       print *, ' register ', appname, " on rank ", rank, " pid1 ", pid1
    endif
    if (comm2 /= mpi_comm_null) then
       appname = 'phis2'//char(0)
       ierr = imoab_registerfortranapplication(trim(appname), comm2, compid2, pid2)
       print *, ' register ', appname, " on rank ", rank, " pid2 ", pid2
    endif
    
 
    if (comm1 /= mpi_comm_null) then
       filename = 'spherecube.h5m'//char(0)
       readopts = 'PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS'//char(0)
       if (rank .eq. sz-2 ) print *, "loading " , trim(filename) , " with options " , trim(readopts)
       nghlay = 0
 
       ierr = imoab_loadmesh(pid1, trim(filename), trim(readopts), nghlay)
       if (rank .eq. sz-1 ) print *, "loaded in parallel ", trim(filename), " error: ", ierr
       ierr = imoab_sendmesh(pid1, gcomm, group2, compid2, repart_scheme); ! send to component 2
       call errorout(ierr, 'cannot send elements' )
    endif
 
    if (comm2 /= mpi_comm_null) then
       ierr = imoab_receivemesh(pid2, gcomm, group1, compid1); ! receive from component 1
       call errorout(ierr, 'cannot receive elements' )
    endif
 
    ! we can now free the sender buffers
    if (comm1 /= mpi_comm_null) then
 
       ierr = imoab_freesenderbuffers(pid1, context_id)
    endif
    call mpi_barrier(gcomm, ierr)
    call errorout(ierr, 'cannot stop at barrier' )
 
if (comm2 /= mpi_comm_null) then
       outfile = 'receivedMesh.h5m'//char(0)
       wopts   = 'PARALLEL=WRITE_PART;'//char(0)
!      write out the mesh file to disk
       ierr = imoab_writemesh(pid2, trim(outfile), trim(wopts))
       call errorout(ierr, 'cannot write received mesh' )
    endif
 
 
    if (comm2 /= mpi_comm_null) then
       ierr = imoab_deregisterapplication(pid2)
       call errorout(ierr, 'cannot deregister app 2 receiver' )
    endif
    if (comm1 /= mpi_comm_null) then
       ierr = imoab_deregisterapplication(pid1)
         call errorout(ierr, 'cannot deregister app 1 sender' )
    endif
 
    ierr = imoab_finalize()
    call errorout(ierr, 'did not finalize iMOAB' )
 
    if (mpi_comm_null /= comm1) call mpi_comm_free(comm1, ierr)
    call errorout(ierr, 'did not free comm1' )
 
    if (mpi_comm_null /= comm2) call mpi_comm_free(comm2, ierr)
    call errorout(ierr, 'did not free comm2' )
 
    call mpi_group_free(allgroup, ierr)
    call mpi_group_free(group1, ierr)
    call mpi_group_free(group2, ierr)
    call mpi_comm_free(gcomm, ierr)
 
    call mpi_finalize(ierr)
    call errorout(ierr, 'did not finalize MPI' )
contains
  SUBROUTINE errorout(ierr, message)
  integer ierr
  character*(*) message
  if (ierr.ne.0) then
    print *, message
    call exit (1)
  end if
  end
 
end program migratemesh