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