hecmw_comm_contact_f.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module m_hecmw_contact_comm
 
  use hecmw_util
  use m_hecmw_comm_f
 
  implicit none
 
  private
  public :: hecmwst_contact_comm
  public :: hecmw_contact_comm_init
  public :: hecmw_contact_comm_finalize
  public :: hecmw_contact_comm_reduce_r
  public :: hecmw_contact_comm_bcast_r
  public :: hecmw_contact_comm_reduce_i
  public :: hecmw_contact_comm_bcast_i
  public :: hecmw_contact_comm_allreduce_i
 
  type hecmwst_contact_comm
    private
    integer(kind=kint) :: n_neighbor_pe
    integer(kind=kint), pointer :: neighbor_pe(:)
    integer(kind=kint) :: MPI_COMM
    integer(kind=kint), pointer :: ext_index(:)
    integer(kind=kint), pointer :: ext_item(:)
    integer(kind=kint), pointer :: int_index(:)
    integer(kind=kint), pointer :: int_item(:)
  end type hecmwst_contact_comm
 
  integer(kind=kint), parameter :: op_overwrite = 46810
 
  integer(kind=kint), parameter :: DBG = 0
 
contains
 
  subroutine hecmw_contact_comm_init(conComm, hecMESH, ndof, n_contact_dof, contact_dofs)
    implicit none
    type (hecmwst_contact_comm), intent(inout) :: concomm
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: ndof, n_contact_dof
    integer(kind=kint), intent(in) :: contact_dofs(:)
    integer(kind=kint), pointer :: ext_index(:), ext_item(:), int_index(:), int_item(:)
    integer(kind=kint), allocatable :: n_ext_per_dom(:), n_int_per_dom(:), ext_item_remote(:)
    integer(kind=kint), allocatable :: statuses(:,:), requests(:)
    integer(kind=kint) :: nn_int, np, ilag, icontact, irow, irank, idom, tag, idof, idx, irow_remote
    integer(kind=kint) :: n_send, is, ie, len
    if (hecmesh%n_neighbor_pe == 0) then
      concomm%n_neighbor_pe = 0
      return
    endif
    nn_int = hecmesh%nn_internal
    np = hecmesh%n_node
    ! count external contact dofs
    allocate(n_ext_per_dom(hecmesh%n_neighbor_pe))
    n_ext_per_dom(:) = 0
    do ilag = 1, n_contact_dof
      icontact = contact_dofs(ilag)
      if (icontact <= nn_int*ndof) cycle   ! skip internal contact dof
      irow = (icontact+ndof-1) / ndof
      irank = hecmesh%node_ID(2*irow)
      call rank_to_idom(hecmesh, irank, idom)
      n_ext_per_dom(idom) = n_ext_per_dom(idom) + 1
    enddo
    ! send external / recv internal contact dofs
    allocate(statuses(hecmw_status_size, hecmesh%n_neighbor_pe))
    allocate(requests(hecmesh%n_neighbor_pe))
    do idom = 1, hecmesh%n_neighbor_pe
      irank = hecmesh%neighbor_pe(idom)
      tag = 5001
      call hecmw_isend_int(n_ext_per_dom(idom), 1, irank, tag, &
           hecmesh%MPI_COMM, requests(idom))
    enddo
    allocate(n_int_per_dom(hecmesh%n_neighbor_pe))
    do idom = 1, hecmesh%n_neighbor_pe
      irank = hecmesh%neighbor_pe(idom)
      tag = 5001
      call hecmw_recv_int(n_int_per_dom(idom), 1, irank, tag, &
           hecmesh%MPI_COMM, statuses(:,1))
    enddo
    call hecmw_waitall(hecmesh%n_neighbor_pe, requests, statuses)
    ! make index
    allocate(ext_index(0:hecmesh%n_neighbor_pe))
    allocate(int_index(0:hecmesh%n_neighbor_pe))
    ext_index(0) = 0
    int_index(0) = 0
    do idom = 1, hecmesh%n_neighbor_pe
      ext_index(idom) = ext_index(idom-1) + n_ext_per_dom(idom)
      int_index(idom) = int_index(idom-1) + n_int_per_dom(idom)
    enddo
    ! make ext_item
    allocate(ext_item(ext_index(hecmesh%n_neighbor_pe)))
    allocate(ext_item_remote(ext_index(hecmesh%n_neighbor_pe)))
    n_ext_per_dom(:) = 0
    do ilag = 1, n_contact_dof
      icontact = contact_dofs(ilag)
      if (icontact <= nn_int*ndof) cycle   ! skip internal contact dof
      irow = (icontact+ndof-1) / ndof
      idof = icontact - ndof*(irow-1)
      irank = hecmesh%node_ID(2*irow)
      call rank_to_idom(hecmesh, irank, idom)
      n_ext_per_dom(idom) = n_ext_per_dom(idom) + 1
      idx = ext_index(idom-1)+n_ext_per_dom(idom)
      ext_item(idx) = icontact
      irow_remote = hecmesh%node_ID(2*irow-1)
      ext_item_remote(idx) = ndof*(irow_remote-1)+idof
    enddo
    deallocate(n_ext_per_dom)
    deallocate(n_int_per_dom)
    ! send ext_item_remote and recv int_item
    n_send = 0
    do idom = 1, hecmesh%n_neighbor_pe
      irank = hecmesh%neighbor_pe(idom)
      is = ext_index(idom-1)+1
      ie = ext_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      n_send = n_send + 1
      tag = 5002
      call hecmw_isend_int(ext_item_remote(is:ie), len, irank, tag, &
           hecmesh%MPI_COMM, requests(n_send))
    enddo
    allocate(int_item(int_index(hecmesh%n_neighbor_pe)))
    do idom = 1, hecmesh%n_neighbor_pe
      irank = hecmesh%neighbor_pe(idom)
      is = int_index(idom-1)+1
      ie = int_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      tag = 5002
      call hecmw_recv_int(int_item(is:ie), len, irank, tag, &
           hecmesh%MPI_COMM, statuses(:,1))
    enddo
    call hecmw_waitall(n_send, requests, statuses)
    deallocate(statuses, requests)
    if (dbg >= 2) then
      write(0,*) '  DEBUG2: ext_index',ext_index(:)
      write(0,*) '  DEBUG2: ext_item',ext_item(:)
      write(0,*) '  DEBUG2: ext_item_remote',ext_item_remote(:)
      write(0,*) '  DEBUG2: int_index',int_index(:)
      write(0,*) '  DEBUG2: int_item',int_item(:)
    endif
    deallocate(ext_item_remote)
    !
    concomm%n_neighbor_pe = hecmesh%n_neighbor_pe
    allocate(concomm%neighbor_pe(concomm%n_neighbor_pe))
    concomm%neighbor_pe(:) = hecmesh%neighbor_pe(:)
    concomm%MPI_COMM = hecmesh%MPI_COMM
    concomm%ext_index => ext_index
    concomm%ext_item => ext_item
    concomm%int_index => int_index
    concomm%int_item => int_item
  end subroutine hecmw_contact_comm_init
 
  subroutine hecmw_contact_comm_finalize(conComm)
    implicit none
    type (hecmwst_contact_comm), intent(inout) :: concomm
    if (concomm%n_neighbor_pe == 0) return
    if (associated(concomm%neighbor_pe)) deallocate(concomm%neighbor_pe)
    if (associated(concomm%ext_index)) deallocate(concomm%ext_index)
    if (associated(concomm%ext_item)) deallocate(concomm%ext_item)
    if (associated(concomm%int_index)) deallocate(concomm%int_index)
    if (associated(concomm%int_item)) deallocate(concomm%int_item)
    concomm%n_neighbor_pe = 0
    concomm%MPI_COMM = 0
  end subroutine hecmw_contact_comm_finalize
 
  subroutine hecmw_contact_comm_reduce_r(conComm, vec, op)
    implicit none
    type (hecmwst_contact_comm), intent(in) :: concomm
    real(kind=kreal), intent(inout) :: vec(:)
    integer(kind=kint), intent(in) :: op
    if (concomm%n_neighbor_pe == 0) return
    call send_recv_contact_info_r(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
         concomm%ext_index, concomm%ext_item, concomm%int_index, concomm%int_item, vec, op)
  end subroutine hecmw_contact_comm_reduce_r
 
  subroutine hecmw_contact_comm_bcast_r(conComm, vec)
    implicit none
    type (hecmwst_contact_comm), intent(in) :: concomm
    real(kind=kreal), intent(inout) :: vec(:)
    integer(kind=kint) :: op
    if (concomm%n_neighbor_pe == 0) return
    op = op_overwrite
    call send_recv_contact_info_r(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
         concomm%int_index, concomm%int_item, concomm%ext_index, concomm%ext_item, vec, op)
  end subroutine hecmw_contact_comm_bcast_r
 
  subroutine hecmw_contact_comm_reduce_i(conComm, vec, op)
    implicit none
    type (hecmwst_contact_comm), intent(in) :: concomm
    integer(kind=kint), intent(inout) :: vec(:)
    integer(kind=kint), intent(in) :: op
    if (concomm%n_neighbor_pe == 0) return
    call send_recv_contact_info_i(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
         concomm%ext_index, concomm%ext_item, concomm%int_index, concomm%int_item, vec, op)
  end subroutine hecmw_contact_comm_reduce_i
 
  subroutine hecmw_contact_comm_bcast_i(conComm, vec)
    implicit none
    type (hecmwst_contact_comm), intent(in) :: concomm
    integer(kind=kint), intent(inout) :: vec(:)
    integer(kind=kint) :: op
    if (concomm%n_neighbor_pe == 0) return
    op = op_overwrite
    call send_recv_contact_info_i(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
         concomm%int_index, concomm%int_item, concomm%ext_index, concomm%ext_item, vec, op)
  end subroutine hecmw_contact_comm_bcast_i
 
  subroutine hecmw_contact_comm_allreduce_i(conComm, vec, op)
    implicit none
    type (hecmwst_contact_comm), intent(in) :: concomm
    integer(kind=kint), intent(inout) :: vec(:)
    integer(kind=kint), intent(in) :: op
    call hecmw_contact_comm_reduce_i(concomm, vec, op)
    call hecmw_contact_comm_bcast_i(concomm, vec)
  end subroutine hecmw_contact_comm_allreduce_i
 
  !
  ! private subroutines
  !
 
  subroutine rank_to_idom(hecMESH, rank, idom)
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in) :: rank
    integer(kind=kint), intent(out) :: idom
    integer(kind=kint) :: i
    do i = 1, hecmesh%n_neighbor_pe
      if (hecmesh%neighbor_pe(i) == rank) then
        idom = i
        return
      endif
    enddo
    stop 'ERROR: exp_rank not found in neighbor_pe'
  end subroutine rank_to_idom
 
  subroutine send_recv_contact_info_r(n_neighbor_pe, neighbor_pe, MPI_COMM, &
       send_index, send_item, recv_index, recv_item, vec, op)
    implicit none
    integer(kind=kint), intent(in) :: n_neighbor_pe
    integer(kind=kint), intent(in) :: neighbor_pe(:)
    integer(kind=kint), intent(in) :: mpi_comm
    integer(kind=kint), pointer, intent(in) :: send_index(:), send_item(:), recv_index(:), recv_item(:)
    real(kind=kreal), intent(inout) :: vec(:)
    integer(kind=kint), intent(in) :: op
    real(kind=kreal), allocatable :: send_buf(:), recv_buf(:)
    integer(kind=kint) :: i, n_send, idom, irank, is, ie, len, tag
    integer(kind=kint), allocatable :: requests(:), statuses(:,:)
    if (n_neighbor_pe == 0) return
    allocate(requests(n_neighbor_pe))
    allocate(statuses(hecmw_status_size, n_neighbor_pe))
    allocate(send_buf(send_index(n_neighbor_pe)))
    allocate(recv_buf(recv_index(n_neighbor_pe)))
    do i = 1, send_index(n_neighbor_pe)
      send_buf(i) = vec(send_item(i))
    enddo
    n_send = 0
    do idom = 1, n_neighbor_pe
      irank = neighbor_pe(idom)
      is = send_index(idom-1)+1
      ie = send_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      n_send = n_send + 1
      tag = 5011
      call hecmw_isend_r(send_buf(is:ie), len, irank, tag, &
           mpi_comm, requests(n_send))
    enddo
    do idom = 1, n_neighbor_pe
      irank = neighbor_pe(idom)
      is = recv_index(idom-1)+1
      ie = recv_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      tag = 5011
      call hecmw_recv_r(recv_buf(is:ie), len, irank, tag, &
           mpi_comm, statuses(:,1))
    enddo
    call hecmw_waitall(n_send, requests, statuses)
    if (op == hecmw_sum) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = vec(recv_item(i)) + recv_buf(i)
      enddo
    elseif (op == hecmw_prod) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = vec(recv_item(i)) * recv_buf(i)
      enddo
    elseif (op == hecmw_max) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = max(vec(recv_item(i)), recv_buf(i))
      enddo
    elseif (op == hecmw_min) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = min(vec(recv_item(i)), recv_buf(i))
      enddo
    else  ! overwrite
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = recv_buf(i)
      enddo
    endif
    deallocate(requests)
    deallocate(statuses)
    if (dbg >= 2) then
      write(0,*) '  DEBUG2: send_buf',send_buf(:)
      write(0,*) '  DEBUG2: recv_buf',recv_buf(:)
    endif
    deallocate(send_buf)
    deallocate(recv_buf)
  end subroutine send_recv_contact_info_r
 
  subroutine send_recv_contact_info_i(n_neighbor_pe, neighbor_pe, MPI_COMM, &
       send_index, send_item, recv_index, recv_item, vec, op)
    implicit none
    integer(kind=kint), intent(in) :: n_neighbor_pe
    integer(kind=kint), intent(in) :: neighbor_pe(:)
    integer(kind=kint), intent(in) :: mpi_comm
    integer(kind=kint), pointer, intent(in) :: send_index(:), send_item(:), recv_index(:), recv_item(:)
    integer(kind=kint), intent(inout) :: vec(:)
    integer(kind=kint), intent(in) :: op
    integer(kind=kint), allocatable :: send_buf(:), recv_buf(:)
    integer(kind=kint) :: i, n_send, idom, irank, is, ie, len, tag
    integer(kind=kint), allocatable :: requests(:), statuses(:,:)
    if (n_neighbor_pe == 0) return
    allocate(requests(n_neighbor_pe))
    allocate(statuses(hecmw_status_size, n_neighbor_pe))
    allocate(send_buf(send_index(n_neighbor_pe)))
    allocate(recv_buf(recv_index(n_neighbor_pe)))
    do i = 1, send_index(n_neighbor_pe)
      send_buf(i) = vec(send_item(i))
    enddo
    n_send = 0
    do idom = 1, n_neighbor_pe
      irank = neighbor_pe(idom)
      is = send_index(idom-1)+1
      ie = send_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      n_send = n_send + 1
      tag = 5011
      call hecmw_isend_int(send_buf(is:ie), len, irank, tag, &
           mpi_comm, requests(n_send))
    enddo
    do idom = 1, n_neighbor_pe
      irank = neighbor_pe(idom)
      is = recv_index(idom-1)+1
      ie = recv_index(idom)
      len = ie-is+1
      if (len == 0) cycle
      tag = 5011
      call hecmw_recv_int(recv_buf(is:ie), len, irank, tag, &
           mpi_comm, statuses(:,1))
    enddo
    call hecmw_waitall(n_send, requests, statuses)
    if (op == hecmw_sum) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = vec(recv_item(i)) + recv_buf(i)
      enddo
    elseif (op == hecmw_prod) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = vec(recv_item(i)) * recv_buf(i)
      enddo
    elseif (op == hecmw_max) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = max(vec(recv_item(i)), recv_buf(i))
      enddo
    elseif (op == hecmw_min) then
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = min(vec(recv_item(i)), recv_buf(i))
      enddo
    else  ! overwrite
      do i = 1, recv_index(n_neighbor_pe)
        vec(recv_item(i)) = recv_buf(i)
      enddo
    endif
    deallocate(requests)
    deallocate(statuses)
    if (dbg >= 2) then
      write(0,*) '  DEBUG2: send_buf',send_buf(:)
      write(0,*) '  DEBUG2: recv_buf',recv_buf(:)
    endif
    deallocate(send_buf)
    deallocate(recv_buf)
  end subroutine send_recv_contact_info_i
 
end module m_hecmw_contact_comm