scale_meshfieldcomm_base.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Data / Communication base
!!
!! @par Description
!!      Base module to manage data communication for element-based methods
!!
!! @author Yuta Kawai, Team SCALE
!<
#include "scaleFElib.h"
module scale_meshfieldcomm_base
 
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc, only: &
    prc_abort
 
  use scale_element_base, only: elementbase
  use scale_mesh_base, only: meshbase
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield1d, meshfield2d, meshfield3d
  use scale_localmesh_base, only: &
    localmeshbase
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  ! 
 
  !> Derived type to manage data communication at a face between adjacent local meshes
  type, public :: localmeshcommdata
    class(localmeshbase), pointer :: lcmesh
    real(rp), allocatable :: send_buf(:,:) !< Buffer for sending data
    real(rp), allocatable :: recv_buf(:,:) !< Buffer for receiving data
    integer :: nnode_lcmeshface !< Number of nodes at a face
    integer :: s_tileid       !< Destination tile ID
    integer :: s_faceid       !< Destination face ID
    integer :: s_panelid      !< Destination panel ID
    integer :: s_rank         !< Destination MPI rank
    integer :: s_tilelocalid  !< Destination local tile ID 
    integer :: faceid         !< Own face ID
  contains
    procedure, public :: init => localmeshcommdata_init
    procedure, public :: final => localmeshcommdata_final
    procedure, public :: sendrecv => localmeshcommdata_sendrecv
    procedure, public :: pc_init_recv => localmeshcommdata_pc_init_recv
    procedure, public :: pc_init_send => localmeshcommdata_pc_init_send
  end type
 
  !> Base derived type to manage data communication
  type, abstract, public :: meshfieldcommbase
    integer :: sfield_num        !< Number of scalar fields
    integer :: hvfield_num       !< Number of horizontal vector fields
    integer :: htensorfield_num  !< Number of horizontal tensor fields
    integer :: field_num_tot     !< Total number of fields
    integer :: nfaces_comm       !< Number of faces where halo data is communicated
 
    integer :: bufsize_per_field !< Buffer size per a field data
 
    class(meshbase), pointer :: mesh
    real(rp), allocatable :: send_buf(:,:,:) !< Buffer for sending data
    real(rp), allocatable :: recv_buf(:,:,:) !< Buffer for receiving data
    integer, allocatable :: request_send(:)
    integer, allocatable :: request_recv(:)
 
    type(localmeshcommdata), allocatable :: commdata_list(:,:)
    integer, allocatable :: is_f(:,:)
    integer, allocatable :: nnode_lcmeshallface(:)
 
    logical :: mpi_pc_flag                !< Flag whether persistent communication is used
    logical :: use_mpi_pc_fujitsu_ext     !< Flag whether Fujitsu extension routines are used for persistent communication
    integer, allocatable :: request_pc(:)
 
    integer :: req_counter
    logical :: call_wait_flag_sub_get     !< Flag whether MPI_wait need to be called before getting halo data from recv_buf
 
    integer :: obj_ind
 
    !-
    logical :: use_vmap_wide_flag
    integer, allocatable :: vmapb_size(:)
    integer, allocatable :: vmapb2(:)
  contains
    procedure(meshfieldcommbase_put), public, deferred :: put   
    procedure(meshfieldcommbase_get), public, deferred :: get
    procedure(meshfieldcommbase_exchange), public, deferred :: exchange
    procedure, public :: prepare_pc => meshfieldcommbase_prepare_pc
  end type meshfieldcommbase
 
  public :: meshfieldcommbase_init, meshfieldcommbase_final
  public :: meshfieldcommbase_exchange_core
  public :: meshfieldcommbase_wait_core
  public :: meshfieldcommbase_extract_bounddata
  public :: meshfieldcommbase_extract_bounddata2
  public :: meshfieldcommbase_extract_bounddata_2
  public :: meshfieldcommbase_extract_bounddata_3
  public :: meshfieldcommbase_set_bounddata
  public :: meshfieldcommbase_set_bounddata_3
 
  !> Container to save a pointer of MeshField(1D, 2D, 3D) object
  type, public :: meshfieldcontainer
    class(meshfield1d), pointer :: field1d
    class(meshfield2d), pointer :: field2d
    class(meshfield3d), pointer :: field3d        
  end type
  
  interface
    subroutine meshfieldcommbase_put(this, field_list, varid_s)
      import meshfieldcommbase
      import meshfieldcontainer
      class(meshfieldcommbase), intent(inout) :: this
      type(meshfieldcontainer), intent(in) :: field_list(:)
      integer, intent(in) :: varid_s
    end subroutine meshfieldcommbase_put
 
    subroutine meshfieldcommbase_get(this, field_list, varid_s)
      import meshfieldcommbase
      import meshfieldcontainer
      class(meshfieldcommbase), intent(inout) :: this
      type(meshfieldcontainer), intent(inout) :: field_list(:)
      integer, intent(in) :: varid_s
    end subroutine meshfieldcommbase_get
 
    subroutine meshfieldcommbase_exchange(this, do_wait)
        import meshfieldcommbase
        class(meshfieldcommbase), intent(inout), target :: this
        logical, intent(in), optional :: do_wait
    end subroutine meshfieldcommbase_exchange
 
    subroutine meshfieldcomm_final(this)
      import meshfieldcommbase  
      class(meshfieldcommbase), intent(inout) :: this
    end subroutine meshfieldcomm_final    
  end interface
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
 
  integer            :: obj_ind       = 0
  integer, parameter :: OBJ_INDEX_MAX = 8192
 
contains
 
!> Initialize a base object to manage data communication of fields
!!
!! @param sfield_num Number of scalar fields
!! @param hvfield_num Number of horizontal vector fields
!! @param htensorfield_num Number of horizontal tensor fields
!! @param bufsize_per_field Buffer size per a field
!! @param comm_face_num Number of faces on a local mesh which perform data communication
!! @param Nnode_LCMeshFace Array to store the number of nodes  
!OCL SERIAL
  subroutine meshfieldcommbase_init( this, &
    sfield_num, hvfield_num, htensorfield_num, bufsize_per_field, comm_face_num, &
    Nnode_LCMeshFace, mesh )
 
    implicit none
    
    class(meshfieldcommbase), intent(inout) :: this
    class(meshbase), intent(in), target :: mesh
    integer, intent(in) :: sfield_num
    integer, intent(in) :: hvfield_num
    integer, intent(in) :: htensorfield_num
    integer, intent(in) :: bufsize_per_field
    integer, intent(in) :: comm_face_num
    integer, intent(in) :: nnode_lcmeshface(comm_face_num,mesh%local_mesh_num)
 
    integer :: n, f
    class(localmeshbase), pointer :: lcmesh
    !-----------------------------------------------------------------------------
 
    !$acc enter data create(this)
    this%mesh => mesh
    this%sfield_num = sfield_num
    this%hvfield_num = hvfield_num
    this%htensorfield_num = htensorfield_num
    this%field_num_tot = sfield_num + hvfield_num*2 + htensorfield_num*4
    this%nfaces_comm = comm_face_num
    !$acc update device(this%sfield_num, this%hvfield_num, this%htensorfield_num, this%field_num_tot, this%nfaces_comm)
    !$acc enter data attach(this%mesh)
 
    if (this%field_num_tot > 0) then
      allocate( this%send_buf(bufsize_per_field, this%field_num_tot, mesh%LOCAL_MESH_NUM) )
      allocate( this%recv_buf(bufsize_per_field, this%field_num_tot, mesh%LOCAL_MESH_NUM) )
      allocate( this%request_send(comm_face_num*mesh%LOCAL_MESH_NUM) )
      allocate( this%request_recv(comm_face_num*mesh%LOCAL_MESH_NUM) )
      !$acc enter data create(this%send_buf, this%recv_buf, this%request_send, this%request_recv)
 
      allocate( this%commdata_list(comm_face_num,mesh%LOCAL_MESH_NUM) )
      allocate( this%is_f(comm_face_num,mesh%LOCAL_MESH_NUM) ) 
      allocate( this%Nnode_LCMeshAllFace(mesh%LOCAL_MESH_NUM) )
      !$acc enter data create(this%is_f, this%Nnode_LCMeshAllFace)
 
      do n=1, mesh%LOCAL_MESH_NUM
        this%is_f(1,n) = 1
        do f=2, this%nfaces_comm
          this%is_f(f,n) = this%is_f(f-1,n) + nnode_lcmeshface(f-1,n)
        end do
        this%Nnode_LCMeshAllFace(n) = sum(nnode_lcmeshface(:,n))
 
        call mesh%GetLocalMesh(n, lcmesh)
        do f=1, this%nfaces_comm
          call this%commdata_list(f,n)%Init( this, lcmesh, f, nnode_lcmeshface(f,n) )
        end do
      end do
      !$acc update device(this%is_f, this%Nnode_LCMeshAllFace)
      !$acc enter data copyin(this%commdata_list)
#ifdef _OPENACC
      do n=1, mesh%LOCAL_MESH_NUM
      do f=1, this%nfaces_comm
        call localmeshcommdata_setupforgpu( this%commdata_list(f,n) )
      end do
      end do
#endif
    end if 
 
    this%MPI_pc_flag = .false.
 
    obj_ind = obj_ind + 1
    if ( obj_ind > obj_index_max ) then
      log_error("MeshFieldCommBase_Init",*) 'obj_ind > OBJ_INDEX_MAX. Check!'
      call prc_abort
    end if
    this%obj_ind = obj_ind
    
    return
  end subroutine meshfieldcommbase_init
 
!> Finalize a base object to manage data communication of fields
!!
!OCL SERIAL
  subroutine meshfieldcommbase_final( this )
    implicit none
    
    class(meshfieldcommbase), intent(inout) :: this
 
    integer :: n, f
    integer :: ireq
    integer :: ierr
    !-----------------------------------------------------------------------------
 
    if (this%field_num_tot > 0) then
      !$acc exit data delete(this%send_buf, this%recv_buf)
      deallocate( this%send_buf, this%recv_buf )
      !$acc exit data delete(this%request_send, this%request_recv)
      deallocate( this%request_send, this%request_recv )
 
      do n=1, this%mesh%LOCAL_MESH_NUM
      do f=1, this%nfaces_comm
        call this%commdata_list(f,n)%Final()
      end do
      end do     
      !$acc exit data delete(this%commdata_list, this%is_f, this%Nnode_LCMeshAllFace)
      deallocate( this%commdata_list, this%is_f, this%Nnode_LCMeshAllFace )
 
      if ( this%MPI_pc_flag ) then
        do ireq=1, this%req_counter
          call mpi_request_free( this%request_pc(ireq), ierr )
        end do           
        deallocate( this%request_pc ) 
      end if
    end if
 
    if ( allocated(this%VMapB_size) ) then
      !$acc exit data delete(this%VMapB_size)
      deallocate(this%VMapB_size)
    end if
    if ( this%use_vmap_wide_flag ) then
      !$acc exit data delete(this%VMapB2)
      deallocate( this%VMapB2 )
    end if
 
    return
  end subroutine meshfieldcommbase_final
 
!> Prepare persistent communication
!!
!! @param use_mpi_pc_fujitsu_ext Flag whether the extension routines of Fujitsu MPI are used
!!
!OCL SERIAL
  subroutine meshfieldcommbase_prepare_pc( this, &
    use_mpi_pc_fujitsu_ext )
    implicit none    
    class(meshfieldcommbase), intent(inout) :: this
    logical, intent(in), optional :: use_mpi_pc_fujitsu_ext
 
    integer :: n, f
    !-----------------------------------------------------------------------------
 
    this%MPI_pc_flag = .true.
 
#ifdef __FUJITSU
    this%use_mpi_pc_fujitsu_ext = .true.
#else
    this%use_mpi_pc_fujitsu_ext = .false.
#endif
    if ( present(use_mpi_pc_fujitsu_ext) ) then
      this%use_mpi_pc_fujitsu_ext = use_mpi_pc_fujitsu_ext
#ifndef __FUJITSU
      if ( use_mpi_pc_fujitsu_ext ) then
        log_error("MeshFieldCommBase_prepare_PC",*) 'use_mpi_pc_fujitsu_ext=.true., but Fujitsu MPI is unavailable. Check!'
        call prc_abort
      end if
#endif
    end if
 
    if (this%field_num_tot > 0) then
      allocate( this%request_pc(2*this%nfaces_comm*this%mesh%LOCAL_MESH_NUM) )
    end if
 
    this%req_counter = 0
    do n=1, this%mesh%LOCAL_MESH_NUM
    do f=1, this%nfaces_comm
      call this%commdata_list(f,n)%PC_Init_recv( &
        this%req_counter, this%request_pc(:),     & ! (inout)
        this%obj_ind, this%use_mpi_pc_fujitsu_ext ) ! (in)
    end do
    end do    
    do n=1, this%mesh%LOCAL_MESH_NUM
    do f=1, this%nfaces_comm
      call this%commdata_list(f,n)%PC_Init_send( &
        this%req_counter, this%request_pc(:),     & ! (inout)
        this%obj_ind, this%use_mpi_pc_fujitsu_ext ) ! (in)
    end do
    end do
 
    return
  end subroutine meshfieldcommbase_prepare_pc
 
!> Exchange halo data
!!
!! @param commdata_list Array of LocalMeshCommData objects which manage information and halo data 
!! @param do_wait Flag whether MPI_waitall is called and move tmp data of LocalMeshCommData object to a recv buffer
!OCL SERIAL
  subroutine meshfieldcommbase_exchange_core( this, commdata_list, do_wait )
#ifdef __FUJITSU
     use mpi_ext, only: &
       fjmpi_prequest_startall
#endif
!    use mpi, only: &
!      MPI_startall
    use scale_prof
    implicit none
  
    class(meshfieldcommbase), intent(inout) :: this
    type(localmeshcommdata), intent(inout), target :: commdata_list(this%nfaces_comm, this%mesh%local_mesh_num)
    logical, intent(in), optional :: do_wait
  
    integer :: n, f    
    integer :: ierr
    logical :: do_wait_
 
    integer :: i
    type(localmeshcommdata), pointer :: lcommdata
    !-----------------------------------------------------------------------------
 
    if ( present(do_wait) ) then
      do_wait_ = do_wait
    else
      do_wait_ = .true.
    end if    
    
!    call PROF_rapstart( 'meshfiled_comm_ex_core', 3)
    !
    if ( this%MPI_pc_flag ) then
#ifdef _OPENACC
      do n=1, this%mesh%LOCAL_MESH_NUM      
      do f=1, this%nfaces_comm
        lcommdata => commdata_list(f,n)
        if ( lcommdata%s_rank /= lcommdata%lcmesh%PRC_myrank ) then
          !$acc update host(lcommdata%send_buf) async(1)
        end if
      end do
      end do
      !$acc wait(1)
 
#endif
      !$omp parallel
      !$omp master
      if ( this%use_mpi_pc_fujitsu_ext ) then
#ifdef __FUJITSU
        ! Use Fujitsu MPI extension for persistent communication
        call fjmpi_prequest_startall( this%req_counter, this%request_pc(1:this%req_counter), ierr )
#endif
      else
        call mpi_startall( this%req_counter, this%request_pc(1:this%req_counter), ierr )
      end if
      !$omp end master
      
      !$omp do collapse(2) private(lcommdata)
      do n=1, this%mesh%LOCAL_MESH_NUM      
      do f=1, this%nfaces_comm
        lcommdata => commdata_list(f,n)
        if ( lcommdata%s_rank == lcommdata%lcmesh%PRC_myrank ) then
#ifdef _OPENACC
          call set_recvbuf_from_sendbuf_lc( commdata_list(abs(lcommdata%s_faceID), lcommdata%s_tilelocalID)%recv_buf, &
                  lcommdata%send_buf, size(lcommdata%send_buf) )          
#else
          commdata_list(abs(lcommdata%s_faceID), lcommdata%s_tilelocalID)%recv_buf(:,:) &
            = lcommdata%send_buf(:,:)
#endif
        end if
      end do 
      end do
      !$acc wait(1)
      !$omp end parallel
    else
      this%req_counter = 0
      do n=1, this%mesh%LOCAL_MESH_NUM      
      do f=1, this%nfaces_comm
        call commdata_list(f,n)%SendRecv( &
          this%req_counter, this%request_send, this%request_recv, & ! (inout)
          commdata_list                                           ) ! (inout) 
      end do
      end do
      !$acc wait(1)
    end if
 
!    call PROF_rapend( 'meshfiled_comm_ex_core', 3)
 
    if ( do_wait_ ) then
      call meshfieldcommbase_wait_core( this, commdata_list )
      this%call_wait_flag_sub_get = .false.
    else
      this%call_wait_flag_sub_get = .true.
    end if
 
    return
  end subroutine meshfieldcommbase_exchange_core
 
!> Wait data communication and move tmp data of LocalMeshCommData object to a recv buffer
!! 
!! @param commdata_list Array of LocalMeshCommData objects which manage information and halo data 
!OCL SERIAL
  subroutine meshfieldcommbase_wait_core( this, commdata_list, &
    field_list, dim, varid_s, lcmesh_list )
    use mpi, only: &
!      MPI_waitall,    &
      mpi_status_size
    use scale_prof
    implicit none
  
    class(meshfieldcommbase), intent(inout) :: this
    type(localmeshcommdata), intent(inout), target :: commdata_list(this%nfaces_comm, this%mesh%local_mesh_num)
    type(meshfieldcontainer), intent(inout), optional :: field_list(:)
    integer, intent(in), optional :: dim
    integer, intent(in), optional :: varid_s
    class(localmeshbase), intent(in), optional, target :: lcmesh_list(:)
 
    integer :: ierr
    integer :: stat_send(mpi_status_size, this%req_counter)
    integer :: stat_recv(mpi_status_size, this%req_counter)
    integer :: stat_pc(mpi_status_size, this%req_counter)
 
    integer :: n, f
    integer :: i, var_id
    integer :: irs(this%nfaces_comm,this%mesh%local_mesh_num), ire(this%nfaces_comm,this%mesh%local_mesh_num)
 
    class(localmeshbase), pointer :: lcmesh
    integer :: val_size(this%mesh%local_mesh_num)
    !----------------------------
 
!   call PROF_rapstart( 'meshfiled_comm_wait_core', 2)
    if ( this%MPI_pc_flag ) then
      if (this%req_counter > 0) then
        call mpi_waitall( this%req_counter, this%request_pc(1:this%req_counter), stat_pc, ierr )
      end if
    else
      if ( this%req_counter > 0 ) then
        call mpi_waitall( this%req_counter, this%request_recv(1:this%req_counter), stat_recv, ierr )
        call mpi_waitall( this%req_counter, this%request_send(1:this%req_counter), stat_send, ierr )
      end if
    end if
 
#ifdef _OPENACC
    do n=1, this%mesh%LOCAL_MESH_NUM
    do f=1, this%nfaces_comm
        if ( commdata_list(f,n)%s_rank /= commdata_list(f,n)%lcmesh%PRC_myrank ) then
          !$acc update device( commdata_list(f,n)%recv_buf ) async(1)
        end if
    end do
    end do
#endif
 
!   call PROF_rapend( 'meshfiled_comm_wait_core', 2)
  
    !---------------------
 
!   call PROF_rapstart( 'meshfiled_comm_wait_post', 2)
 
    if ( present(field_list) ) then
      do n=1, this%mesh%LOCAL_MESH_NUM
        lcmesh => lcmesh_list(n)
        val_size(n) = lcmesh%refElem%Np * lcmesh%NeA
        irs(1,n) = lcmesh%refElem%Np * lcmesh%Ne + 1
        do f=1, this%nfaces_comm
          ire(f,n) = irs(f,n) + commdata_list(f,n)%Nnode_LCMeshFace - 1
          if (f<this%nfaces_comm) irs(f+1,n) = ire(f,n) + 1
        end do
      end do
      
      !$omp parallel do private(var_id,n,i,f) collapse(3)
      !$acc parallel loop gang collapse(3) present(field_list, commdata_list) copyin(irs, ire, val_size) async(1)
      do n=1, this%mesh%LOCAL_MESH_NUM
      do i=1, size(field_list)
      do f=1, this%nfaces_comm
        var_id = varid_s + i - 1
 
        if (dim==1) then
          call set_bounddata( field_list(var_id)%field1d%local(n)%val, val_size(n), irs(f,n), ire(f,n), commdata_list(f,n)%recv_buf(:,var_id) )
        else if (dim==2) then
          call set_bounddata( field_list(var_id)%field2d%local(n)%val, val_size(n), irs(f,n), ire(f,n), commdata_list(f,n)%recv_buf(:,var_id) )
        else if (dim==3) then
          call set_bounddata( field_list(var_id)%field3d%local(n)%val, val_size(n), irs(f,n), ire(f,n), commdata_list(f,n)%recv_buf(:,var_id) )
        end if
      end do ! end loop for face
      end do
      end do
#ifdef _OPENACC
      do n=1, this%mesh%LOCAL_MESH_NUM
      do i=1, size(field_list)
      do f=1, this%nfaces_comm
        var_id = varid_s + i - 1
        if (dim==1) then
          !$acc update device( field_list(var_id)%field1d%local(n)%val(irs(f,n):ire(f,n)) ) async(1)
        else if (dim==2) then
          !$acc update device( field_list(var_id)%field2d%local(n)%val(irs(f,n):ire(f,n)) ) async(1)
        else if (dim==3) then
          !$acc update device( field_list(var_id)%field3d%local(n)%val(irs(f,n):ire(f,n)) ) async(1)
        end if
      end do ! end loop for face
      end do
      end do
#endif
 
    else
      
      do n=1, this%mesh%LOCAL_MESH_NUM
        irs(1,n) = 1
        do f=1, this%nfaces_comm
          ire(f,n) = irs(f,n) + commdata_list(f,n)%Nnode_LCMeshFace - 1
          if (f<this%nfaces_comm) irs(f+1,n) = ire(f,n) + 1
        end do
      end do
 
      !$omp parallel do private(n,var_id,f) collapse(3)
      !$acc parallel loop gang collapse(3) present(this%recv_buf, commdata_list) copyin(irs, ire) async(1)
      do n=1, this%mesh%LOCAL_MESH_NUM
      do var_id=1, this%field_num_tot
      do f=1, this%nfaces_comm
#ifdef _OPENACC
        call set_bounddata( this%recv_buf(:,var_id,n), size(this%recv_buf(:,var_id,n)), irs(f,n), ire(f,n), commdata_list(f,n)%recv_buf(:,var_id) )
#else        
        this%recv_buf(irs(f,n):ire(f,n),var_id,n) = commdata_list(f,n)%recv_buf(:,var_id)
#endif
      end do ! end loop for face
      end do
      end do
    end if
 
    !$acc wait(1)
!   call PROF_rapend( 'meshfiled_comm_wait_post', 2)
    return
  contains
!OCL SERIAL
    subroutine set_bounddata( var, IA, irs_, ire_, recv_buf )
      implicit none
      integer, intent(in) :: IA
      real(RP), intent(inout) :: var(IA)
      integer, intent(in) :: irs_, ire_
      real(RP), intent(in) :: recv_buf(ire_-irs_+1)
 
      integer :: ii
      !-----------------------------
      !$acc routine vector
 
      !$acc loop vector
      do ii=1, size(recv_buf)
        var(irs_+ii-1) = recv_buf(ii)
      end do
      return
    end subroutine set_bounddata
  end subroutine meshfieldcommbase_wait_core
 
!> Extract halo data from data array with MeshField object and set it to the receiving buffer
!OCL SERIAL
  subroutine meshfieldcommbase_extract_bounddata(var, refElem, mesh, buf)
    implicit none
  
    class(elementbase), intent(in) :: refelem
    class(localmeshbase), intent(in) :: mesh
    real(rp), intent(in) :: var(refelem%np * mesh%nea)
    real(rp), intent(out) :: buf(size(mesh%vmapb))
 
    integer :: i
    !-----------------------------------------------------------------------------
    !$omp parallel do
    !$acc parallel loop present(var, mesh%VmapB, buf)
!OCL PREFETCH
    do i=1, size(buf)
      buf(i) = var(mesh%VmapB(i))
    end do
    return
  end subroutine meshfieldcommbase_extract_bounddata
 
!> Extract halo data from data array with MeshField object and set it to the receiving buffer
!OCL SERIAL
  subroutine meshfieldcommbase_extract_bounddata2(var, VMapB, VMapB_size, NpxNeA, buf)
    implicit none
    integer, intent(in) :: vmapb_size
    integer, intent(in) :: npxnea
    real(rp), intent(in) :: var(npxnea)
    integer, intent(in) :: vmapb(vmapb_size)
    real(rp), intent(out) :: buf(vmapb_size)
 
    integer :: i
    !-----------------------------------------------------------------------------
    !$omp parallel do
    !$acc parallel loop present(var, VmapB, buf) async(1)
!OCL PREFETCH
    do i=1, size(buf)
      buf(i) = var(vmapb(i))
    end do
    return
  end subroutine meshfieldcommbase_extract_bounddata2  
 
!> Extract halo data from data array with MeshField object and set it to the receiving buffer
!!
!! Note: We assume that the size of VmapB is the same for all local meshes. 
!! For the future, this subroutine should be modified to handle the case with different sizes of VmapB among local meshes.
!!
!OCL SERIAL
  subroutine meshfieldcommbase_extract_bounddata_2(field_list, dim, varid_s, lcmesh_list, vmapB_size, buf)
    implicit none
    type(meshfieldcontainer), intent(in), target :: field_list(:)
    integer, intent(in) :: varid_s                                               !< Starting variable ID to extract halo data
    integer, intent(in) :: dim                                                   !< Number of dimensions of the field (1, 2, or 3)
    class(localmeshbase), intent(in), target :: lcmesh_list(:)                   !< Array of local mesh objects to extract halo data
    integer, intent(in) :: vmapb_size                                            !< Size of VmapB (we assume that the size of VmapB is the same for all local meshes)
    real(rp), intent(out) :: buf(vmapb_size,size(field_list),size(lcmesh_list))
 
    class(localmeshbase), pointer :: lcmesh
    integer :: varid
    integer :: i, n
 
    integer :: field_num
    !-----------------------------------------------------------------------------
 
    field_num = size(field_list)
 
    do n=1, size(lcmesh_list)
      lcmesh => lcmesh_list(n)
      i = 1
      do while( i <= field_num )
        varid = varid_s + i - 1
        if ( i+1 <= field_num ) then
          if (dim==1) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), field_list(varid)%field1d%local(n)%val, field_list(varid+1)%field1d%local(n)%val, lcmesh, lcmesh%refElem )
          else if(dim==2) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), field_list(varid)%field2d%local(n)%val, field_list(varid+1)%field2d%local(n)%val, lcmesh, lcmesh%refElem )
          else if(dim==3) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), field_list(varid)%field3d%local(n)%val, field_list(varid+1)%field3d%local(n)%val, lcmesh, lcmesh%refElem )
          end if
          i = i + 2
        else
          if (dim==1) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field1d%local(n)%val, lcmesh%VMapB, size(lcmesh%VMapB), lcmesh%refElem%Np*lcmesh%NeA,  buf(:,varid,n) )
          else if(dim==2) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field2d%local(n)%val, lcmesh%VMapB, size(lcmesh%VMapB), lcmesh%refElem%Np*lcmesh%NeA,  buf(:,varid,n) )
          else if(dim==3) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field3d%local(n)%val, lcmesh%VMapB, size(lcmesh%VMapB), lcmesh%refElem%Np*lcmesh%NeA,  buf(:,varid,n) )
          end if
          i = i + 1
        end if
      end do
      !!$acc update host( buf(:,varid_s:field_num,n) )
    end do
    !$acc wait(1)
 
    return
  contains
!OCL SERIAL
    subroutine extract_bounddata_var2( buf1_, buf2_, var1, var2, lmesh, elem )
      implicit none
      class(localmeshbase), intent(in) :: lmesh
      class(elementbase), intent(in) :: elem
      real(rp), intent(out) :: buf1_(size(lmesh%vmapb))
      real(rp), intent(out) :: buf2_(size(lmesh%vmapb))
      real(rp), intent(inout) :: var1(elem%np*lmesh%nea)
      real(rp), intent(inout) :: var2(elem%np*lmesh%nea)
 
      integer :: ii, b
      !-----------------------------
      !$omp parallel do
      !$acc parallel loop present(buf1_, buf2_, var1, var2, lmesh%vmapB) async(1)
!OCL PREFETCH
      do ii=1, size(buf1_)
        b = lmesh%vmapB(ii)
        buf1_(ii) = var1(b)
        buf2_(ii) = var2(b)
      end do
      return
    end subroutine extract_bounddata_var2
  end subroutine meshfieldcommbase_extract_bounddata_2
 
!> Extract halo data from data array with MeshField object and set it to the recieving buffer
!OCL SERIAL
  subroutine meshfieldcommbase_extract_bounddata_3(field_list, dim, varid_s, lcmesh_list, VMapB2, VMapB2_size, buf)
    implicit none
    class(meshfieldcontainer), intent(in), target :: field_list(:)
    integer, intent(in) :: varid_s
    integer, intent(in) :: dim
    class(localmeshbase), intent(in), target :: lcmesh_list(:)
    integer, intent(in) :: vmapb2_size
    integer, intent(in) :: vmapb2(vmapb2_size)
    real(rp), intent(out) :: buf(vmapb2_size,size(field_list),size(lcmesh_list))
 
    class(localmeshbase), pointer :: lcmesh
    integer :: varid
    integer :: i, n
 
    integer :: field_num
    !-----------------------------------------------------------------------------
 
    field_num = size(field_list)
 
    do n=1, size(lcmesh_list)
      lcmesh => lcmesh_list(n)
      i = 1
      do while( i <= field_num )
        varid = varid_s + i - 1
        if ( i+1 <= field_num ) then
          if (dim==1) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), &
              vmapb2, vmapb2_size, field_list(varid)%field1d%local(n)%val, field_list(varid+1)%field1d%local(n)%val, lcmesh%refElem%Np * lcmesh%NeA )
          else if(dim==2) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), &
              vmapb2, vmapb2_size, field_list(varid)%field2d%local(n)%val, field_list(varid+1)%field2d%local(n)%val, lcmesh%refElem%Np * lcmesh%NeA )
          else if(dim==3) then
            call extract_bounddata_var2( buf(:,varid,n), buf(:,varid+1,n), &
              vmapb2, vmapb2_size, field_list(varid)%field3d%local(n)%val, field_list(varid+1)%field3d%local(n)%val, lcmesh%refElem%Np * lcmesh%NeA )
          end if
          i = i + 2
        else
          if (dim==1) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field1d%local(n)%val, vmapb2, vmapb2_size, lcmesh%refElem%Np*lcmesh%NeA,  &
              buf(:,varid,n) )
          else if(dim==2) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field2d%local(n)%val, vmapb2, vmapb2_size, lcmesh%refElem%Np*lcmesh%NeA,  &
              buf(:,varid,n) )
          else if(dim==3) then
            call meshfieldcommbase_extract_bounddata2( field_list(varid)%field3d%local(n)%val, vmapb2, vmapb2_size, lcmesh%refElem%Np*lcmesh%NeA,  &
              buf(:,varid,n) )
          end if
          i = i + 1
        end if
      end do
    end do
    return
  contains
!OCL SERIAL
    subroutine extract_bounddata_var2( buf1_, buf2_, VMapB, VMapB_size, var1, var2, varSize )
      implicit none
      integer, intent(in) :: vmapb_size
      real(rp), intent(out) :: buf1_(vmapb_size)
      real(rp), intent(out) :: buf2_(vmapb_size)
      integer, intent(in) :: vmapb(vmapb_size)
      integer, intent(in) :: varsize
      real(rp), intent(inout) :: var1(varsize)
      real(rp), intent(inout) :: var2(varsize)
 
      integer :: ii
      !-----------------------------
      !$omp parallel do
!OCL PREFETCH
      do ii=1, vmapb_size
        buf1_(ii) = var1(vmapb(ii))
        buf2_(ii) = var2(vmapb(ii))
      end do
      return
    end subroutine extract_bounddata_var2
  end subroutine meshfieldcommbase_extract_bounddata_3
 
!> Extract halo data from the receiving buffer and set it to data array with MeshField object
  subroutine meshfieldcommbase_set_bounddata(buf, refElem, mesh, var)
    implicit none
    
    class(elementbase), intent(in) :: refelem
    class(localmeshbase), intent(in) :: mesh
    real(rp), intent(in) :: buf(size(mesh%vmapb))
    real(rp), intent(inout) :: var(refelem%np * mesh%nea)
 
    integer :: ii, iis
    !-----------------------------------------------------------------------------
 
#ifdef _OPENACC
     iis = refelem%Np*mesh%NeE
     !$acc parallel loop present(buf, var) async(1)
     do ii=1, size(buf)
       var(iis+ii) = buf(ii)
     end do
#else
    var(refelem%Np*mesh%NeE+1:refelem%Np*mesh%NeE+size(buf)) = buf(:)
#endif
    return
  end subroutine meshfieldcommbase_set_bounddata  
 
!> Extract halo data from the recieving buffer and set it to data array with MeshField object
  subroutine meshfieldcommbase_set_bounddata_3(buf, refElem, mesh, vmapB, vmapB_size, var)
    implicit none
    
    class(elementbase), intent(in) :: refelem
    class(localmeshbase), intent(in) :: mesh
    integer, intent(in) :: vmapb_size
    integer, intent(in) :: vmapb(vmapb_size)
    real(rp), intent(in) :: buf(vmapb_size)
    real(rp), intent(inout) :: var(refelem%np * mesh%nea)
    !-----------------------------------------------------------------------------
 
    var(refelem%Np*mesh%NeE+1:refelem%Np*mesh%NeE+size(buf)) = buf(:)
    return
  end subroutine meshfieldcommbase_set_bounddata_3
  
  !-------------------------------------------
  
  !> Initialize an object to manage data communication of fields for a face on a local mesh
  subroutine localmeshcommdata_init( this, comm, lcmesh, faceID, Nnode_LCMeshFace )
    implicit none
 
    class(localmeshcommdata), intent(inout) :: this
    class(meshfieldcommbase), intent(in) :: comm
    class(localmeshbase), intent(in), target :: lcmesh
    integer, intent(in) :: faceid
    integer, intent(in) :: nnode_lcmeshface
    !-----------------------------------------------------------------------------
 
    this%lcmesh => lcmesh
    this%Nnode_LCMeshFace = nnode_lcmeshface
 
    !-
    this%s_tileID  = comm%mesh%tileID_globalMap(faceid, lcmesh%tileID)
    this%s_faceID  = comm%mesh%tileFaceID_globalMap(faceid, lcmesh%tileID)
    this%s_panelID = comm%mesh%tilePanelID_globalMap(faceid, lcmesh%tileID)      
    this%s_rank    = comm%mesh%PRCrank_globalMap(this%s_tileID)
    this%s_tilelocalID = comm%mesh%tileID_global2localMap(this%s_tileID)
    this%faceID    = faceid
 
    allocate( this%send_buf(nnode_lcmeshface, comm%field_num_tot))
    allocate( this%recv_buf(nnode_lcmeshface, comm%field_num_tot))
    return
  end subroutine localmeshcommdata_init
 
  subroutine localmeshcommdata_setupforgpu( this )
    implicit none
    type(localmeshcommdata), intent(inout) :: this
    !-----------------------------------------------------------------------------
    !$acc enter data create(this%send_buf, this%recv_buf)
    !$acc enter data attach(this%lcmesh)
    return
  end subroutine localmeshcommdata_setupforgpu
 
  !> Send and receive halo data for a face on a local mesh
  subroutine localmeshcommdata_sendrecv( this, &
    req_counter, req_send, req_recv,           &
    lccommdat_list )
    
    use scale_prc, only: &
      prc_local_comm_world
    use mpi, only: &
      !MPI_isend, MPI_irecv, &
      mpi_double_precision
    implicit none
 
    class(localmeshcommdata), intent(inout) :: this
    integer, intent(inout) :: req_counter
    integer, intent(inout) :: req_send(:)
    integer, intent(inout) :: req_recv(:)
    type(localmeshcommdata), intent(inout) :: lccommdat_list(:,:)
  
    integer :: tag
    integer :: bufsize
    integer :: ierr
    !-------------------------------------------
 
    if ( this%s_rank /= this%lcmesh%PRC_myrank ) then
 
      !$acc update host(this%send_buf) async(1)
      !$acc wait(1)
 
      req_counter = req_counter + 1
 
      tag = 10 * this%lcmesh%tileID + this%faceID
      bufsize = size(this%recv_buf)
      call mpi_irecv( this%recv_buf(1,1), bufsize, mpi_double_precision, &
        this%s_rank, tag, prc_local_comm_world,                          &
        req_recv(req_counter), ierr)
 
      bufsize = size(this%send_buf)
      tag = 10 * this%s_tileID + abs(this%s_faceID)
      call mpi_isend( this%send_buf(1,1), bufsize, mpi_double_precision, &
       this%s_rank, tag, prc_local_comm_world,                           &
       req_send(req_counter), ierr )
      
    else if ( this%s_rank == this%lcmesh%PRC_myrank ) then
#ifdef _OPENACC
      call set_recvbuf_from_sendbuf_lc( lccommdat_list(abs(this%s_faceID), this%s_tilelocalID)%recv_buf, &
        this%send_buf, size(this%send_buf) )
#else      
      lccommdat_list(abs(this%s_faceID), this%s_tilelocalID)%recv_buf(:,:) &
        = this%send_buf(:,:)
#endif
    end if         
    
    return
  end subroutine localmeshcommdata_sendrecv
#ifdef _OPENACC
  subroutine set_recvbuf_from_sendbuf_lc( recv_buf, send_buf, buf_size )
      implicit none
      integer, intent(in) :: buf_size
      real(rp), intent(out) :: recv_buf(buf_size)
      real(rp), intent(in) :: send_buf(buf_size)
 
      integer :: i
      !-----------------------------
      !$acc parallel loop present(recv_buf, send_buf) async(1)
      do i=1, buf_size
        recv_buf(i) = send_buf(i)
      end do
      return
  end subroutine set_recvbuf_from_sendbuf_lc
#endif
 
  !> Initialize persistent communication for sending halo data
  subroutine localmeshcommdata_pc_init_send( this, &
    req_counter, req, obj_ind_ ,      &
    use_mpi_pc_fujisu_ext             )
    
    use scale_prc, only: &
      prc_local_comm_world
    use mpi, only: &
      mpi_double_precision!, &
!      MPI_send_init
#ifdef __FUJITSU
    use mpi_ext, only: &
      fjmpi_prequest_send_init
#endif
    implicit none
 
    class(localmeshcommdata), intent(inout) :: this
    integer, intent(inout) :: req_counter
    integer, intent(inout) :: req(:)
    integer, intent(in) :: obj_ind_
    logical, intent(in) :: use_mpi_pc_fujisu_ext
 
    integer :: tag
    integer :: bufsize
    integer :: ierr
    !-------------------------------------------
 
    if ( this%s_rank /= this%lcmesh%PRC_myrank ) then
      req_counter = req_counter + 1
 
      bufsize = size(this%send_buf)
!      tag = 1000 * this%s_tileID + 10 * obj_ind_ + abs(this%s_faceID)
      tag = obj_ind_
 
      if ( use_mpi_pc_fujisu_ext ) then
#ifdef __FUJITSU
        call fjmpi_prequest_send_init( this%send_buf(1,1), bufsize, mpi_double_precision, &
          this%s_rank, tag, prc_local_comm_world,                                         &
          req(req_counter), ierr )
#endif
      else
        call mpi_send_init( this%send_buf(1,1), bufsize, mpi_double_precision, &
          this%s_rank, tag, prc_local_comm_world,                              &
          req(req_counter), ierr )
      end if
    end if         
 
    return
  end subroutine localmeshcommdata_pc_init_send
  !> Initialize persistent communication for receiving halo data
  subroutine localmeshcommdata_pc_init_recv( this, &
    req_counter, req, obj_ind_,       &
    use_mpi_pc_fujisu_ext             )
    
    use scale_prc, only: &
      prc_local_comm_world
    use mpi, only: &
      mpi_double_precision!, &
!      MPI_recv_init
#ifdef __FUJITSU
    use mpi_ext, only: &
      fjmpi_prequest_recv_init
#endif
    implicit none
 
    class(localmeshcommdata), intent(inout) :: this
    integer, intent(inout) :: req_counter
    integer, intent(inout) :: req(:)
    integer, intent(in) :: obj_ind_
    logical, intent(in) :: use_mpi_pc_fujisu_ext
  
    integer :: tag
    integer :: bufsize
    integer :: ierr
    !-------------------------------------------
 
    if ( this%s_rank /= this%lcmesh%PRC_myrank ) then
      req_counter = req_counter + 1
 
      bufsize = size(this%recv_buf)
!      tag = 1000 * this%lcmesh%tileID + 10 * obj_ind_ + this%faceID
      tag = obj_ind_
 
      if ( use_mpi_pc_fujisu_ext ) then
#ifdef __FUJITSU
        call fjmpi_prequest_recv_init( this%recv_buf(1,1), bufsize, mpi_double_precision, &
          this%s_rank, tag, prc_local_comm_world,                                         &
          req(req_counter), ierr )
#endif
      else 
        call mpi_recv_init( this%recv_buf(1,1), bufsize, mpi_double_precision, &
          this%s_rank, tag, prc_local_comm_world,                              &
          req(req_counter), ierr )
      end if
    end if         
 
    return
  end subroutine localmeshcommdata_pc_init_recv
 
  !> Finalize an object to manage data communication of fields for a face on a local mesh
  subroutine localmeshcommdata_final( this )
    implicit none
    
    class(localmeshcommdata), intent(inout) :: this
    !-----------------------------------------------------------------------------
 
    !$acc exit data delete(this%send_buf, this%recv_buf)
    deallocate( this%send_buf )
    deallocate( this%recv_buf )
    return
  end subroutine localmeshcommdata_final
 
end module scale_meshfieldcomm_base