advection.F90

Go to the documentation of this file.

! advection: do a one time step transport, using c binding module
! example of how to do it in parallel in fortran
! This program shows how to load in parallel in moab from Fortran90, and how
! to call the transport wrapper. It is a fortran equivalent of intx_imesh
! driver
!
! Usage: advection
 
#define CHECK(a) if (ierr .ne. 0) print *, a
program advection
 
 use iso_c_binding
 implicit none
 
#include "moab/MOABConfig.h"
#ifdef MOAB_HAVE_MPI
# include "mpif.h"
# include "iMeshP_f.h"
#else
# include "iMesh_f.h"
#endif
 
! extern void update_tracer( iMesh_Instance instance, iBase_EntitySetHandle * opEulerSet, int * ierr);
 INTERFACE
 SUBROUTINE update_tracer ( instance , opEulerSet, ierr) bind(C)
 use iso_c_binding
 implicit none
 imesh_instance, INTENT(IN) , VALUE :: instance
 ibase_entitysethandle, INTENT(IN), VALUE :: opeulerset
 integer(c_int) , INTENT (OUT) :: ierr
 END SUBROUTINE update_tracer
 END INTERFACE
 
 ! declarations
 ! imesh is the instance handle
 imesh_instance imesh
 
 ! cells will be storing 2d cells
 TYPE(c_ptr) cells_ptr
 ibase_entityhandle, pointer :: cells(:)
 INTEGER ents_alloc, ents_size
 
 ibase_entitysethandle root_set
 ibase_entitysethandle opeulerset
 CHARACTER (LEN=200) options
 CHARACTER (LEN=200) filename
 CHARACTER (LEN=200) optionswrite
 CHARACTER (LEN=200) outname
 ! TYPE(C_PTR) :: vertsPtr, entsPtr
 
 integer rank, sz, ierr
 integer lenopt, lenname
 integer islist
 
#ifdef MOAB_HAVE_MPI
 ! local variables for parallel runs
 imeshp_partitionhandle imeshp
 ibase_handle_t mpi_comm_c
#endif
 
 ! init the parallel partition
 call mpi_init(ierr)
 check("fail to initialize MPI")
 call mpi_comm_size(mpi_comm_world, sz, ierr)
 check("fail to get MPI size")
 call mpi_comm_rank(mpi_comm_world, rank, ierr)
 check("fail to get MPI rank")
 
 ! now load the mesh; this also initializes parallel sharing
 imesh = 0
 imeshp = 0
 call imesh_newmesh("MOAB", imesh, ierr)
 check("fail to initialize imesh")
 
 call imesh_getrootset(%VAL(imesh), root_set, ierr)
 check("fail to get root set")
 
 call imeshp_getcommunicator(%VAL(imesh), mpi_comm_world, mpi_comm_c, ierr)
 
 call imeshp_createpartitionall(%VAL(imesh), %VAL(mpi_comm_c), imeshp, ierr)
 check("fail to create parallel partition ")
 options = " moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION" // &
 " moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE "
! " moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION " &
! " moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE ", & ! options
 if (0 .eq. rank) then
 print *, "load in parallel file HN16DP.h5m"
 endif
 lenname = 11;
 lenopt = 123
 filename = "HN16DP.h5m"
 call imeshp_loadall(%VAL(imesh), &
 %VAL(imeshp), &
 %VAL(root_set), &
 filename, & ! filename
 options, & !options
 ierr, &
 %VAL(lenname), & ! strlen(filename),
 %VAL(lenopt) ) !119) !strlen(options));
 check("fail to load mesh in parallel ")
 
 islist = 0
 call imesh_createentset(%VAL(imesh), %VAL(islist), opeulerset, ierr)
 check("fail to create euler set ")
 
 ents_alloc = 0
 ents_size = 0
 
 call imesh_getentities(%VAL(imesh),%VAL(root_set), &
 %VAL(ibase_face), &
 %VAL(imesh_all_topologies), cells_ptr, &
 ents_alloc, ents_size, ierr);
 check("fail to get 2d entities ")
 
 call c_f_pointer(cells_ptr, cells, [ents_size])
 
 call imesh_addentarrtoset(%VAL(imesh), cells, %VAL(ents_size), &
 %VAL(opeulerset), ierr)
 
 call update_tracer(imesh, opeulerset, ierr)
 check("fail to update tracer ")
 
 outname = "outF.h5m";
 optionswrite = " moab:PARALLEL=WRITE_PART " ;
 lenname = 8
 lenopt = 27
 call imeshp_saveall( %VAL(imesh), &
 %VAL(imeshp), &
 %VAL(opeulerset), &
 outname, &
 optionswrite, &
 ierr, &
 %VAL(lenname), & ! strlen(filename),
 %VAL(lenopt) ) !119) !strlen(options));
 check(" can't save ")
 
 if (0==rank) then
 print *, "Done"
 endif
 
 call mpi_finalize(ierr)
 stop
end program advection