boundaries.F90

Go to the documentation of this file.

!! Copyright (C) 2005-2010 Florian Lorenzen, Heiko Appel, X. Andrade
!! Copyright (C) 2021 Sebastian Ohlmann
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module boundaries_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use global_oct_m
  use index_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mpi_oct_m
  use mpi_debug_oct_m
  use namespace_oct_m
  use par_vec_oct_m
  use partition_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  !
  type boundaries_t
    private
    logical              :: periodic = .false.       
    logical              :: fully_periodic = .false. 
    integer              :: nper = 0
    integer, allocatable :: per_points(:, :)
    integer, allocatable :: per_send(:, :)
    integer, allocatable :: per_recv(:, :)
    logical              :: per_contiguous
    integer, allocatable :: nsend(:)
    integer, allocatable :: nrecv(:)
    type(accel_mem_t)    :: buff_per_points
    type(accel_mem_t)    :: buff_per_send
    type(accel_mem_t)    :: buff_per_recv
    type(accel_mem_t)    :: buff_nsend
    type(accel_mem_t)    :: buff_nrecv
    logical, public      :: spiralBC = .false. 
    logical, public      :: spiral   = .false. 
    real(real64),  allocatable, public :: spiral_q(:)
  end type boundaries_t
 
  public ::                        &
    boundaries_t,                  &
    boundaries_init,               &
    boundaries_end,                &
    boundaries_set
 
  public ::                        &
    par_vec_handle_batch_t,        &
    dpar_vec_ghost_update,         &
    zpar_vec_ghost_update,         &
    dghost_update_batch_start,     &
    zghost_update_batch_start,     &
    dghost_update_batch_finish,    &
    zghost_update_batch_finish
 
  integer, parameter, public ::    &
    POINT_BOUNDARY = 1,            &
    point_inner    = 2
 
  type par_vec_handle_batch_t
    private
    type(batch_t)        :: ghost_send
    type(MPI_Request), allocatable :: requests(:)
    integer              :: nnb
    ! these are needed for CL
    real(real64), pointer       :: drecv_buffer(:)
    complex(real64), pointer       :: zrecv_buffer(:)
    real(real64), pointer       :: dsend_buffer(:)
    complex(real64), pointer       :: zsend_buffer(:)
    type(batch_t),   pointer :: v_local
    type(par_vec_t), pointer :: pv
  end type par_vec_handle_batch_t
 
  interface boundaries_set
    module procedure boundaries_set_batch
    module procedure dboundaries_set_single
    module procedure zboundaries_set_single
  end interface boundaries_set
 
contains
 
  ! ---------------------------------------------------------
  subroutine boundaries_init(this, namespace, space, mesh, qvector)
    type(boundaries_t),   intent(inout) :: this
    type(namespace_t),    intent(in)    :: namespace
    class(space_t),       intent(in)    :: space
    type(mesh_t), target, intent(in)    :: mesh
    real(real64), optional,      intent(in)    :: qvector(:)
 
    integer :: sp, ip, ip_inner, iper, ip_bnd
    integer(int64) :: ip_inner_global
    integer :: ipart
    integer(int64), allocatable :: recv_rem_points(:, :), points(:), per_send(:, :)
    integer, allocatable :: part(:), points_local(:)
    integer :: nper_recv, iper_recv
    integer, allocatable :: rdispls(:), sdispls(:)
 
    push_sub(boundaries_init)
 
 
    this%periodic = space%is_periodic()
    this%fully_periodic = space%periodic_dim == space%dim
 
    if (space%is_periodic()) then
 
      !%Variable SpiralBoundaryCondition
      !%Type logical
      !%Default no
      !%Section Mesh
      !%Description
      !% (Experimental) If set to yes, Octopus will apply spin-spiral boundary conditions.
      !% The momentum of the spin spiral is defined by the variable
      !% <tt>TDMomentumTransfer</tt>
      !%End
      call parse_variable(namespace, 'SpiralBoundaryCondition', .false., this%spiralBC)
      if (this%spiralBC) then
        call messages_experimental("SpiralBoundaryCondition")
        if(.not. present(qvector)) then
          message(1) = "TDMomentumTransfer or TDReducedMomentumTransfer must be defined if SpiralBoundaryCondition=yes"
          call messages_fatal(1, namespace=namespace)
        end if
        safe_allocate(this%spiral_q(1:space%dim))
        this%spiral_q(1:space%dim) = qvector(1:space%dim)
      end if
 
      sp = mesh%np
      if (mesh%parallel_in_domains) sp = mesh%np + mesh%pv%np_ghost
 
      ! count the number of points that are periodic
      this%nper = 0
      nper_recv = 0
      do ip = sp + 1, mesh%np_part
        ip_inner_global = mesh_periodic_point(mesh, space, ip)
        ip_inner = mesh_global2local(mesh, ip_inner_global)
 
        ! it is the same point, can happen for mixed periodicity
        if (ip == ip_inner) cycle
        ! the point maps to the boundary, can happen for mixed periodicity
        ! in this case the point is already set to zero, so we can ignore it
        ! for different mixed boundary conditions, we would need to be careful here
        if (ip_inner_global > mesh%np_global) cycle
        ! now check if point is local or if it needs to be communicated
        if (ip_inner /= 0 .and. ip_inner <= mesh%np) then
          this%nper = this%nper + 1
        else
          nper_recv = nper_recv + 1
        end if
      end do
      if (.not. mesh%parallel_in_domains) then
        assert(nper_recv == 0)
      end if
 
      safe_allocate(this%per_points(1:2, 1:max(this%nper, 1)))
      !$omp parallel do
      do ip = 1, this%nper
        this%per_points(1:2, ip) = -1
      end do
 
      if (mesh%parallel_in_domains) then
        safe_allocate(this%per_recv(1:max(nper_recv, 1), 1:max(mesh%pv%npart, 1)))
        safe_allocate(this%nrecv(1:mesh%pv%npart))
        safe_allocate(recv_rem_points(1:nper_recv, 1:mesh%pv%npart))
        safe_allocate(points(1:nper_recv))
        safe_allocate(points_local(1:nper_recv))
        safe_allocate(part(1:nper_recv))
        this%nrecv = 0
      end if
 
      iper = 0
      iper_recv = 0
      do ip = sp + 1, mesh%np_part
        ip_inner_global = mesh_periodic_point(mesh, space, ip)
        ip_inner = mesh_global2local(mesh, ip_inner_global)
 
        ! it is the same point, can happen for mixed periodicity
        if (ip == ip_inner) cycle
        ! the point maps to the boundary, can happen for mixed periodicity
        ! in this case the point is already set to zero, so we can ignore it
        if (ip_inner_global > mesh%np_global) cycle
        ! now check if point is local or if it needs to be communicated
        if (ip_inner /= 0 .and. ip_inner <= mesh%np) then
          iper = iper + 1
          this%per_points(point_boundary, iper) = ip
          this%per_points(point_inner, iper) = ip_inner
        else
          ! this can only happen if parallel in domain
          ! the point is on another node
          iper_recv = iper_recv + 1
          points(iper_recv) = ip_inner_global
          points_local(iper_recv) = ip
        end if
      end do
      if (mesh%parallel_in_domains) then
        ! find the points in the other partitions
        call partition_get_partition_number(mesh%partition, nper_recv, points, part)
        do iper_recv = 1, nper_recv
          ipart = part(iper_recv)
          assert(this%nrecv(ipart) + 1 < huge(0_int32))
          ! count the points to receive from each node
          this%nrecv(ipart) = this%nrecv(ipart) + 1
          ! and store the number of the point
          this%per_recv(this%nrecv(ipart), ipart) = points_local(iper_recv)
          ! and its global index
          recv_rem_points(this%nrecv(ipart), ipart) = points(iper_recv)
        end do
      end if
 
      if (mesh%parallel_in_domains) then
        ! communicate the number of points to receive/send
        safe_allocate(this%nsend(1:mesh%pv%npart))
        call mesh%mpi_grp%alltoall(this%nrecv, 1, mpi_integer, this%nsend, 1, mpi_integer)
 
        safe_allocate(sdispls(1:mesh%pv%npart))
        safe_allocate(rdispls(1:mesh%pv%npart))
        safe_allocate(this%per_send(1:max(maxval(this%nsend), 1), 1:mesh%pv%npart))
        safe_allocate(per_send(1:maxval(this%nsend), 1:mesh%pv%npart))
        ! compute displacements
        assert(int(nper_recv, int64)*mesh%pv%npart < huge(0_int32))
        do ipart = 1, mesh%pv%npart
          rdispls(ipart) = nper_recv * (ipart - 1)
          sdispls(ipart) = maxval(this%nsend) * (ipart - 1)
        end do
 
        ! exchange indices
        call mesh%mpi_grp%alltoallv(recv_rem_points, this%nrecv, rdispls, mpi_integer8, &
          per_send, this%nsend, sdispls, mpi_integer8)
 
        do ipart = 1, mesh%pv%npart
          ! get local index of the points to send
          do ip = 1, this%nsend(ipart)
            this%per_send(ip, ipart) = mesh_global2local(mesh, per_send(ip, ipart))
            ! make sure we have local points here
            assert(this%per_send(ip, ipart) > 0)
            assert(this%per_send(ip, ipart) <= mesh%np)
          end do
        end do
 
        safe_deallocate_a(per_send)
        safe_deallocate_a(sdispls)
        safe_deallocate_a(rdispls)
        safe_deallocate_a(recv_rem_points)
        safe_deallocate_a(points)
        safe_deallocate_a(points_local)
        safe_deallocate_a(part)
      end if
 
      ! We now determine if the boundary points are contiguous or not
      ! If so, the first mapping boundaries%per_points(POINT_BOUNDARY, ip) is not necessary
      this%per_contiguous = .true.
      do ip = 1, this%nper-1
        ip_bnd = this%per_points(point_boundary, ip)
        if (ip_bnd + 1 /= this%per_points(point_boundary, ip+1)) then
          this%per_contiguous = .false.
          exit
        end if
      end do
 
      if (accel_is_enabled()) then
        call accel_create_buffer(this%buff_per_points, accel_mem_read_only, type_integer, 2*this%nper)
        call accel_write_buffer(this%buff_per_points, 2*this%nper, this%per_points)
 
        if (mesh%parallel_in_domains) then
          call accel_create_buffer(this%buff_per_send, accel_mem_read_only, type_integer, product(ubound(this%per_send)))
          call accel_write_buffer(this%buff_per_send, product(ubound(this%per_send)), this%per_send)
 
          call accel_create_buffer(this%buff_per_recv, accel_mem_read_only, type_integer, product(ubound(this%per_recv)))
          call accel_write_buffer(this%buff_per_recv, product(ubound(this%per_recv)), this%per_recv)
 
          call accel_create_buffer(this%buff_nsend, accel_mem_read_only, type_integer, mesh%pv%npart)
          call accel_write_buffer(this%buff_nsend, mesh%pv%npart, this%nsend)
 
          call accel_create_buffer(this%buff_nrecv, accel_mem_read_only, type_integer, mesh%pv%npart)
          call accel_write_buffer(this%buff_nrecv, mesh%pv%npart, this%nrecv)
        end if
      end if
 
    end if
 
    pop_sub(boundaries_init)
  end subroutine boundaries_init
 
  ! ---------------------------------------------------------
 
  subroutine boundaries_end(this)
    type(boundaries_t),  intent(inout) :: this
 
    push_sub(boundaries_end)
 
    if (this%periodic) then
      safe_deallocate_a(this%spiral_q)
 
      safe_deallocate_a(this%per_send)
      safe_deallocate_a(this%per_recv)
      safe_deallocate_a(this%nsend)
      safe_deallocate_a(this%nrecv)
 
      if (accel_is_enabled()) then
        call accel_release_buffer(this%buff_per_send)
        call accel_release_buffer(this%buff_per_recv)
        call accel_release_buffer(this%buff_nsend)
        call accel_release_buffer(this%buff_nrecv)
      end if
 
      if (accel_is_enabled()) call accel_release_buffer(this%buff_per_points)
 
      safe_deallocate_a(this%per_points)
    end if
 
    pop_sub(boundaries_end)
  end subroutine boundaries_end
 
  ! -------------------------------------------------------
 
  subroutine boundaries_set_batch(this, mesh, ffb, phase_correction, buff_phase_corr, offset, async)
    type(boundaries_t), intent(in)    :: this
    class(mesh_t),      intent(in)    :: mesh
    class(batch_t),     intent(inout) :: ffb
    complex(real64), optional,    intent(in)    :: phase_correction(:)
    type(accel_mem_t), optional,intent(in)    :: buff_phase_corr
    integer, optional,          intent(in)    :: offset
    logical, optional,          intent(in)    :: async
 
 
    push_sub(boundaries_set_batch)
 
    if (ffb%type() == type_float) then
      call dboundaries_set_batch(this, mesh, ffb, phase_correction, buff_phase_corr, offset, async)
    else if (ffb%type() == type_cmplx) then
      call zboundaries_set_batch(this, mesh, ffb, phase_correction, buff_phase_corr, offset, async)
    else
      assert(.false.)
    end if
 
    pop_sub(boundaries_set_batch)
  end subroutine boundaries_set_batch
 
#include "undef.F90"
#include "complex.F90"
#include "boundaries_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "boundaries_inc.F90"
 
end module boundaries_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: