MpiRunControl.F90

Go to the documentation of this file.

module mpiruncontrolmodule
#if defined(__WITH_PETSC__)
#include <petsc/finclude/petscksp.h>
  ! cannot use all symbols because of clash with 'use mpi'
  use petscksp, only: petsc_comm_world, petscinitialize, petscfinalize, &
                      petsc_null_character, petsc_null_options
#endif
  use mpi
  use mpiworldmodule
  use simvariablesmodule, only: proc_id, nr_procs
  use simstagesmodule
  use profilermodule
  use kindmodule, only: i4b, lgp, dp
  use stlvecintmodule
  use numericalsolutionmodule
  use runcontrolmodule, only: runcontroltype
  implicit none
  private
 
  public :: create_mpi_run_control
 
  type, public, extends(runcontroltype) :: mpiruncontroltype
  contains
    ! override
    procedure :: start => mpi_ctrl_start
    procedure :: finish => mpi_ctrl_finish
    procedure :: after_con_cr => mpi_ctrl_after_con_cr
    ! private
    procedure, private :: wait_for_debugger
  end type mpiruncontroltype
 
contains
 
  function create_mpi_run_control() result(controller)
    class(runcontroltype), pointer :: controller
    ! local
    class(mpiruncontroltype), pointer :: mpi_controller
 
    allocate (mpi_controller)
    controller => mpi_controller
 
  end function create_mpi_run_control
 
  subroutine mpi_ctrl_start(this)
    use errorutilmodule, only: pstop_alternative
    ! local
    integer(I4B) :: tmr_init_par
 
    class(mpiruncontroltype) :: this
    ! local
    integer :: ierr
    character(len=*), parameter :: petsc_db_file = '.petscrc'
    logical(LGP) :: petsc_db_exists, wait_dbg, is_parallel_mode
    type(mpiworldtype), pointer :: mpi_world
 
    ! add timed section for parallel initialization
    tmr_init_par = -1
    call g_prof%start("Initialize MPI and PETSc", tmr_init_par)
 
    ! set mpi abort function
    pstop_alternative => mpi_stop
 
    wait_dbg = .false.
    mpi_world => get_mpi_world()
 
    ! if PETSc we need their initialize
#if defined(__WITH_PETSC__)
    ! PetscInitialize calls MPI_Init only when it is not called yet,
    ! which could be through the API. If it is already called, we
    ! should assign the MPI communicator to PETSC_COMM_WORLD first
    ! (PETSc manual)
    if (mpi_world%has_comm()) then
      petsc_comm_world = mpi_world%comm
    end if
 
    inquire (file=petsc_db_file, exist=petsc_db_exists)
    if (.not. petsc_db_exists) then
      call petscinitialize(petsc_null_character, ierr)
      chkerrq(ierr)
    else
      call petscinitialize(petsc_db_file, ierr)
      chkerrq(ierr)
    end if
 
    if (.not. mpi_world%has_comm()) then
      call mpi_world%set_comm(petsc_comm_world)
    end if
 
    call petscoptionshasname(petsc_null_options, petsc_null_character, &
                             '-wait_dbg', wait_dbg, ierr)
    chkerrq(ierr)
    call petscoptionshasname(petsc_null_options, petsc_null_character, &
                             '-p', is_parallel_mode, ierr)
    chkerrq(ierr)
#else
    if (.not. mpi_world%has_comm()) then
      call mpi_init(ierr)
      call check_mpi(ierr)
      call mpi_world%set_comm(mpi_comm_world)
    end if
#endif
 
    call mpi_world%init()
 
    call mpi_comm_size(mpi_world%comm, nr_procs, ierr)
    call mpi_comm_rank(mpi_world%comm, proc_id, ierr)
 
    ! possibly wait to attach debugger here
    if (wait_dbg) call this%wait_for_debugger()
 
    ! done with parallel pre-work
    call g_prof%stop(tmr_init_par)
 
    ! start everything else by calling parent
    call this%RunControlType%start()
 
  end subroutine mpi_ctrl_start
 
  subroutine wait_for_debugger(this)
    class(mpiruncontroltype) :: this
    ! local
    integer :: ierr
    integer(I4B) :: icnt
    type(mpiworldtype), pointer :: mpi_world
 
    mpi_world => get_mpi_world()
    if (proc_id == 0) then
      icnt = 0
      write (*, *) 'Hit enter to continue...'
      read (*, *)
    end if
    call mpi_barrier(mpi_world%comm, ierr)
 
  end subroutine wait_for_debugger
 
  subroutine mpi_ctrl_finish(this)
    use errorutilmodule, only: pstop_alternative
    class(mpiruncontroltype) :: this
    ! local
    integer :: ierr
    integer(I4B) :: tmr_final_par
 
    ! timer
    tmr_final_par = -1
    call g_prof%start("Finalize MPI and PETSc", tmr_final_par)
 
    ! release MPI related memory in router before MPI_Finalize
    call this%virtual_data_mgr%router%finalize()
 
    ! finish mpi
#if defined(__WITH_PETSC__)
    ! NB: PetscFinalize calls MPI_Finalize only when MPI_Init
    ! was called before PetscInitialize
    call petscfinalize(ierr)
    chkerrq(ierr)
#else
    call mpi_finalize(ierr)
    call check_mpi(ierr)
#endif
 
    pstop_alternative => null()
 
    call g_prof%stop(tmr_final_par)
 
    ! finish everything else by calling parent
    call this%RunControlType%finish()
 
  end subroutine mpi_ctrl_finish
 
  !> @brief Actions after creating connections
  !<
  subroutine mpi_ctrl_after_con_cr(this)
    use virtualdatalistsmodule
    use virtualmodelmodule, only: virtualmodeltype, &
                                  get_virtual_model_from_list, &
                                  get_virtual_model
    use virtualexchangemodule, only: virtualexchangetype, &
                                     get_virtual_exchange_from_list, &
                                     get_virtual_exchange
    class(mpiruncontroltype) :: this
    ! local
    integer(I4B) :: i, j, id, irank
    integer(I4B) :: nr_models, nr_exgs, nr_remotes, max_nr_remotes
    type(stlvecint) :: remote_models, remote_exgs
    integer(I4B), dimension(:, :), pointer :: remote_models_per_process
    integer(I4B), dimension(:, :), pointer :: remote_exgs_per_process
    class(virtualmodeltype), pointer :: vm
    class(virtualexchangetype), pointer :: ve
    type(mpiworldtype), pointer :: mpi_world
    integer :: ierr
 
    mpi_world => get_mpi_world()
 
    ! activate halo through base
    call this%RunControlType%after_con_cr()
 
    ! compose list of remote models/exchanges to receive
    call remote_models%init()
    nr_models = virtual_model_list%Count()
    do i = 1, nr_models
      vm => get_virtual_model_from_list(virtual_model_list, i)
      if (vm%is_active .and. .not. vm%is_local) then
        ! remote and active
        call remote_models%push_back(vm%id)
      end if
    end do
    call remote_exgs%init()
    nr_exgs = virtual_exchange_list%Count()
    do i = 1, nr_exgs
      ve => get_virtual_exchange_from_list(virtual_exchange_list, i)
      if (ve%is_active .and. .not. ve%is_local) then
        ! remote and active
        call remote_exgs%push_back(ve%id)
      end if
    end do
 
    ! Models: find max for allocation
    nr_remotes = remote_models%size
    call mpi_allreduce(nr_remotes, max_nr_remotes, 1, mpi_integer, mpi_max, &
                       mpi_world%comm, ierr)
    call check_mpi(ierr)
 
    allocate (remote_models_per_process(max_nr_remotes, nr_procs))
    remote_models_per_process = 0
 
    ! Models: fill local portion and reduce
    do i = 1, remote_models%size
      remote_models_per_process(i, proc_id + 1) = remote_models%at(i)
    end do
    call mpi_allreduce(mpi_in_place, remote_models_per_process, &
                       max_nr_remotes * nr_procs, mpi_integer, mpi_max, &
                       mpi_world%comm, ierr)
    call check_mpi(ierr)
 
    ! Models: set remotes to virtual models
    do i = 1, nr_procs
      do j = 1, max_nr_remotes
        id = remote_models_per_process(j, i)
        if (id > 0) then
          ! assign zero-based rank number to virtual model
          vm => get_virtual_model(id)
          if (vm%is_local) then
            ! only for local models
            irank = i - 1
            call vm%rcv_ranks%push_back_unique(irank)
          end if
        end if
      end do
    end do
 
    ! Exchanges: find max for allocation
    nr_remotes = remote_exgs%size
    call mpi_allreduce(nr_remotes, max_nr_remotes, 1, mpi_integer, mpi_max, &
                       mpi_world%comm, ierr)
    call check_mpi(ierr)
 
    allocate (remote_exgs_per_process(max_nr_remotes, nr_procs))
    remote_exgs_per_process = 0
 
    ! Exchanges: fill local portion and reduce
    do i = 1, remote_exgs%size
      remote_exgs_per_process(i, proc_id + 1) = remote_exgs%at(i)
    end do
    call mpi_allreduce(mpi_in_place, remote_exgs_per_process, &
                       max_nr_remotes * nr_procs, mpi_integer, mpi_max, &
                       mpi_world%comm, ierr)
    call check_mpi(ierr)
 
    ! Exchanges: set remotes to virtual exchanges
    do i = 1, nr_procs
      do j = 1, max_nr_remotes
        id = remote_exgs_per_process(j, i)
        if (id > 0) then
          ! assign zero-based rank number to virtual exchange
          ve => get_virtual_exchange(id)
          if (ve%is_local) then
            ! only for local exchanges
            irank = i - 1
            call ve%rcv_ranks%push_back_unique(irank)
          end if
        end if
      end do
    end do
 
    ! clean up
    call remote_models%destroy()
    call remote_exgs%destroy()
 
    deallocate (remote_models_per_process)
    deallocate (remote_exgs_per_process)
 
  end subroutine mpi_ctrl_after_con_cr
 
end module mpiruncontrolmodule