global_interpolation.F90

Go to the documentation of this file.

! Copyright (c) 2020-2023, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module global_interpolation
  use num_types, only: rp
  use space, only: space_t
  use dofmap, only: dofmap_t
  use mesh, only: mesh_t
  use logger, only: neko_log, log_size
  use utils, only: neko_error, neko_warning
  use local_interpolation
  use comm
  use math, only: copy
  use neko_mpi_types
  use structs, only: array_ptr_t
  use, intrinsic :: iso_c_binding
  implicit none
  private
  type, public :: global_interpolation_t
     type(array_ptr_t) :: x
     type(array_ptr_t) :: y
     type(array_ptr_t) :: z
     integer :: gdim
     integer :: nelv
     type(space_t), pointer :: xh
     type(local_interpolator_t) :: local_interp
     logical :: all_points_local = .false.
     !! Gslib handle.
     !! @note: Remove when we remove gslib.
     integer :: gs_handle
     logical :: gs_init = .false.
     integer :: n_points
     real(kind=rp), allocatable :: xyz(:,:)
     integer, allocatable :: proc_owner(:)
     integer, allocatable :: el_owner(:)
     type(c_ptr) :: el_owner_d = c_null_ptr
     real(kind=rp), allocatable :: rst(:,:)
     real(kind=rp), allocatable :: dist2(:)
     integer, allocatable :: error_code(:)
     real(kind=rp) :: tol = 5d-13
   contains
     procedure, pass(this) :: init_xyz => global_interpolation_init_xyz
     procedure, pass(this) :: init_dof => global_interpolation_init_dof
     procedure, pass(this) :: free => global_interpolation_free
     procedure, pass(this) :: free_points => global_interpolation_free_points
     procedure, pass(this) :: find_points_and_redist => &
                              global_interpolation_find_and_redist
     procedure, pass(this) :: find_points_coords => &
          global_interpolation_find_coords
     procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
     generic :: find_points => find_points_xyz, find_points_coords
     procedure, pass(this) :: evaluate => global_interpolation_evaluate
 
     generic :: init => init_dof, init_xyz
 
  end type global_interpolation_t
 
contains
 
  subroutine global_interpolation_init_dof(this, dof, tol)
    class(global_interpolation_t), intent(inout) :: this
    type(dofmap_t), target :: dof
    real(kind=rp), optional :: tol
 
    ! NOTE: Passing dof%x(:,1,1,1), etc in init_xyz passes down the entire
    ! dof%x array and not a slice. It is done this way for
    ! this%x%ptr to point to dof%x (see global_interpolation_init_xyz).
    call this%init_xyz(dof%x(:,1,1,1), dof%y(:,1,1,1), dof%z(:,1,1,1), &
         dof%msh%gdim, dof%msh%nelv, dof%Xh, tol = tol)
 
  end subroutine global_interpolation_init_dof
 
  subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, tol)
    class(global_interpolation_t), intent(inout) :: this
    real(kind=rp), intent(in), target :: x(:)
    real(kind=rp), intent(in), target :: y(:)
    real(kind=rp), intent(in), target :: z(:)
    integer, intent(in) :: gdim
    integer, intent(in) :: nelv
    type(space_t), intent(in), target :: Xh
    real(kind=rp), intent(in), optional :: tol
    integer :: lx, ly, lz, max_pts_per_iter
 
#ifdef HAVE_GSLIB
 
    this%x%ptr => x
    this%y%ptr => y
    this%z%ptr => z
    this%gdim = gdim
    this%nelv = nelv
    this%Xh => xh
    if (present(tol)) this%tol = tol
 
    ! Number of points to iterate on simultaneosuly
    max_pts_per_iter = 128
    lx = xh%lx
    ly = xh%ly
    lz = xh%lz
 
    call fgslib_findpts_setup(this%gs_handle, &
         neko_comm, pe_size, &
         this%gdim, &
         this%x%ptr, this%y%ptr, this%z%ptr, & ! Physical nodal values
         lx, ly, lz, nelv, & ! Mesh dimensions
         2*lx, 2*ly, 2*lz, & ! Mesh size for bounding box computation
         0.01, & ! relative size to expand bounding boxes by
         lx*ly*lz*nelv, lx*ly*lz*nelv, & ! local/global hash mesh sizes
         max_pts_per_iter, this%tol)
    this%gs_init = .true.
#else
    call neko_error('Neko needs to be built with GSLIB support')
#endif
 
  end subroutine global_interpolation_init_xyz
 
 
  subroutine global_interpolation_free(this)
    class(global_interpolation_t), intent(inout) :: this
 
    nullify(this%x%ptr)
    nullify(this%y%ptr)
    nullify(this%z%ptr)
    nullify(this%Xh)
 
    this%nelv = 0
    this%gdim = 0
 
    call this%free_points()
 
#ifdef HAVE_GSLIB
    if (this%gs_init) then
       call fgslib_findpts_free(this%gs_handle)
       this%gs_init = .false.
    end if
#endif
 
  end subroutine global_interpolation_free
 
  subroutine global_interpolation_free_points(this)
    class(global_interpolation_t), intent(inout) :: this
 
    this%n_points = 0
    this%all_points_local = .false.
 
    if (allocated(this%xyz)) deallocate(this%xyz)
    if (allocated(this%rst)) deallocate(this%rst)
    if (allocated(this%proc_owner)) deallocate(this%proc_owner)
    if (allocated(this%el_owner)) deallocate(this%el_owner)
    if (allocated(this%dist2)) deallocate(this%dist2)
    if (allocated(this%error_code)) deallocate(this%error_code)
 
    if (c_associated(this%el_owner_d)) then
       call device_free(this%el_owner_d)
    end if
 
  end subroutine global_interpolation_free_points
 
  subroutine global_interpolation_find_common(this)
    class(global_interpolation_t), intent(inout) :: this
    !!Perhaps this should be kind dp
    real(kind=rp) :: xdiff, ydiff, zdiff
    character(len=8000) :: log_buf
    real(kind=rp), allocatable :: x_check(:)
    real(kind=rp), allocatable :: y_check(:)
    real(kind=rp), allocatable :: z_check(:)
    logical :: isdiff
    integer :: i
 
 
#ifdef HAVE_GSLIB
 
    ! gslib find points, which element they belong, to process etc.
    call fgslib_findpts(this%gs_handle, &
         this%error_code, 1, &
         this%proc_owner, 1, &
         this%el_owner, 1, &
         this%rst, this%gdim, &
         this%dist2, 1, &
         this%xyz(1,1), this%gdim, &
         this%xyz(2,1), this%gdim, &
         this%xyz(3,1), this%gdim, this%n_points)
 
    do i = 1 , this%n_points
 
       !
       ! Check validity of points
       !
       if (this%error_code(i) .eq. 1) then
          if (this%dist2(i) .gt. this%tol) then
             write(*,*) 'Point with coords: ',&
                this%xyz(1,i),&
                this%xyz(2,i),&
                this%xyz(3,i),&
                'Did not converge to tol. Absolute differences squared: ',&
                this%dist2(i), 'PE rank', pe_rank
          end if
       end if
 
       if (this%error_code(i) .eq. 2) &
             write(*,*) 'Point with coords: ',&
                this%xyz(1,i), this%xyz(2,i), this%xyz(3,i),&
                'Outside the mesh!',&
                ' Interpolation on these points will return 0.0. dist2: ', &
                this%dist2(i),&
                'el_owner, rst coords, pe: ',&
                this%el_owner(i), this%rst(1,i), this%rst(2,i), &
                this%rst(3,i), pe_rank
 
    end do
 
    allocate(x_check(this%n_points))
    allocate(y_check(this%n_points))
    allocate(z_check(this%n_points))
 
    call fgslib_findpts_eval(this%gs_handle, x_check, &
                             1, this%error_code, 1, &
                             this%proc_owner, 1, this%el_owner, 1, &
                             this%rst, this%gdim, &
                             this%n_points, this%x%ptr)
 
    call fgslib_findpts_eval(this%gs_handle, y_check, &
                             1, this%error_code, 1, &
                             this%proc_owner, 1, this%el_owner, 1, &
                             this%rst, this%gdim, &
                             this%n_points, this%y%ptr)
 
    call fgslib_findpts_eval(this%gs_handle, z_check, &
                             1, this%error_code, 1, &
                             this%proc_owner, 1, this%el_owner, 1, &
                             this%rst, this%gdim, &
                             this%n_points, this%z%ptr)
 
    do i = 1 , this%n_points
 
       !
       ! Check validity of points
       !
       isdiff = .false.
       xdiff = (x_check(i)-this%xyz(1,i))**2
       if ( xdiff .gt. this%tol) isdiff = .true.
       ydiff = (y_check(i)-this%xyz(2,i))**2
       if ( ydiff .gt. this%tol) isdiff = .true.
       zdiff = (z_check(i)-this%xyz(3,i))**2
       if ( zdiff .gt. this%tol) isdiff = .true.
 
       if (isdiff) then
          write(*,*) 'Points with coords: ', &
                this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
                'Differ from interpolated coords: ', &
                x_check(i), y_check(i), z_check(i), &
                'Distance squared: ', &
                xdiff, ydiff, zdiff
       end if
 
    end do
 
    deallocate(x_check)
    deallocate(y_check)
    deallocate(z_check)
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%el_owner, this%el_owner_d, &
            this%n_points, host_to_device, sync = .true.)
    end if
#else
    call neko_error('Neko needs to be built with GSLIB support')
#endif
  end subroutine global_interpolation_find_common
 
  subroutine global_interpolation_find_coords(this, x, y, z, n_points)
    class(global_interpolation_t), intent(inout) :: this
    integer :: n_points
    !!Perhaps this should be kind dp
    !!this is to get around that x,y,z often is 4 dimensional...
    !!Should maybe add interface for 1d aswell
    real(kind=rp) :: x(n_points,1,1,1)
    real(kind=rp) :: y(n_points,1,1,1)
    real(kind=rp) :: z(n_points,1,1,1)
    integer :: i
 
    call this%free_points()
 
    this%n_points = n_points
 
    call global_interpolation_init_point_arrays(this)
 
    do i = 1, n_points
       this%xyz(1, i) = x(i,1,1,1)
       this%xyz(2, i) = y(i,1,1,1)
       this%xyz(3, i) = z(i,1,1,1)
    end do
 
    call global_interpolation_find_common(this)
 
  end subroutine global_interpolation_find_coords
 
  subroutine global_interpolation_init_point_arrays(this)
    class(global_interpolation_t) :: this
 
    allocate(this%xyz(3, this%n_points))
    allocate(this%rst(3, this%n_points))
    allocate(this%proc_owner(this%n_points))
    allocate(this%el_owner(this%n_points))
    allocate(this%dist2(this%n_points))
    allocate(this%error_code(this%n_points))
 
    if (neko_bcknd_device .eq. 1) &
       call device_map(this%el_owner, this%el_owner_d, this%n_points)
 
  end subroutine global_interpolation_init_point_arrays
 
  subroutine global_interpolation_find_xyz(this, xyz, n_points)
    class(global_interpolation_t), intent(inout) :: this
    integer, intent(in) :: n_points
    !!Perhaps this should be kind dp
    real(kind=rp), intent(inout) :: xyz(3, n_points)
 
 
    call this%free_points()
 
    this%n_points = n_points
 
    call global_interpolation_init_point_arrays(this)
 
    call copy(this%xyz, xyz, 3 * n_points)
 
    call global_interpolation_find_common(this)
 
  end subroutine global_interpolation_find_xyz
 
  subroutine global_interpolation_find_and_redist(this, xyz, n_points)
    class(global_interpolation_t), intent(inout) :: this
    integer, intent(inout) :: n_points
    !!Perhaps this should be kind dp
    real(kind=rp), allocatable, intent(inout) :: xyz(:,:)
    integer :: i
 
 
    call this%free_points()
 
    this%n_points = n_points
    call global_interpolation_init_point_arrays(this)
 
    call copy(this%xyz, xyz, 3 * n_points)
 
    call global_interpolation_find_common(this)
    call global_interpolation_redist(this)
    call global_interpolation_find_common(this)
 
    do i = 1, this%n_points
       if (this%proc_owner(i) .ne. pe_rank) then
          write(*,*) 'Redistribution failed on rank: ', pe_rank, &
                     'for point with coord: ', &
                     this%xyz(1, i), this%xyz(2, i), this%xyz(3, i)
          exit
       end if
    end do
 
    n_points = this%n_points
    if (allocated(xyz)) then
       deallocate(xyz)
    end if
    allocate(xyz(3, n_points))
    call copy(xyz, this%xyz, 3*n_points)
 
    call this%local_interp%init(this%Xh, this%rst(1,:),&
                                this%rst(2,:), this%rst(3,:), n_points)
    this%all_points_local = .true.
 
 
  end subroutine global_interpolation_find_and_redist
 
  subroutine global_interpolation_redist(this)
    class(global_interpolation_t), intent(inout) :: this
    integer, allocatable :: n_points_per_pe(:)
    integer, allocatable :: n_points_from_pe(:)
    integer, allocatable :: n_point_offset_from_pe(:)
    real(kind=rp), allocatable :: xyz_send_to_pe(:,:)
    real(kind=rp), allocatable :: new_xyz(:,:)
    integer :: i, j, k, ierr, n_new_points, max_n_points_to_send
 
    n_new_points = 0
 
    allocate(n_points_per_pe(0:(pe_size-1)))
    allocate(n_points_from_pe(0:(pe_size-1)))
    n_points_per_pe = 0
    n_points_from_pe = 0
    do i = 1, this%n_points
       n_points_per_pe(this%proc_owner(i)) = &
            n_points_per_pe(this%proc_owner(i)) + 1
    end do
    do i = 0, (pe_size - 1)
       call mpi_reduce(n_points_per_pe(i), n_new_points, 1, mpi_integer, &
            mpi_sum, i, neko_comm, ierr)
       call mpi_gather(n_points_per_pe(i), 1, mpi_integer,&
                      n_points_from_pe, 1, mpi_integer, i, neko_comm, ierr)
    end do
 
    allocate(n_point_offset_from_pe(0:(pe_size-1)))
    n_point_offset_from_pe(0) = 0
    do i = 1, (pe_size - 1)
       n_point_offset_from_pe(i) = n_points_from_pe(i-1)&
                                 + n_point_offset_from_pe(i-1)
    end do
 
    allocate(new_xyz(3, n_new_points))
    max_n_points_to_send = maxval(n_points_per_pe)
    allocate(xyz_send_to_pe(3, max_n_points_to_send))
    do i = 0, (pe_size - 1)
       k = 0
       do j = 1, this%n_points
          if (this%proc_owner(j) .eq. i) then
             k = k + 1
             xyz_send_to_pe(:,k) = this%xyz(:,j)
          end if
       end do
       if (k .ne. n_points_per_pe(i)) then
          write(*,*) 'PE: ', pe_rank, ' has k= ', k, &
                     'points for PE:', i,' but should have: ', &
                     n_points_per_pe(i)
          call neko_error('Error in redistribution of points')
       end if
       call mpi_gatherv(xyz_send_to_pe,3*n_points_per_pe(i), &
                        mpi_double_precision, new_xyz,3*n_points_from_pe, &
                        3*n_point_offset_from_pe, &
                        mpi_double_precision, i, neko_comm, ierr)
 
    end do
 
    call this%free_points()
 
    this%n_points = n_new_points
    call global_interpolation_init_point_arrays(this)
    call copy(this%xyz, new_xyz, 3 * n_new_points)
 
    deallocate(n_point_offset_from_pe)
    deallocate(n_points_from_pe)
    deallocate(n_points_per_pe)
    deallocate(xyz_send_to_pe)
 
  end subroutine global_interpolation_redist
 
 
 
  subroutine global_interpolation_evaluate(this, interp_values, field)
    class(global_interpolation_t), intent(inout) :: this
    real(kind=rp), intent(inout) :: interp_values(this%n_points)
    real(kind=rp), intent(inout) :: field(this%nelv*this%Xh%lxyz)
 
 
#ifdef HAVE_GSLIB
    if (.not. this%all_points_local) then
       call fgslib_findpts_eval(this%gs_handle, interp_values, &
                                1, this%error_code, 1, &
                                this%proc_owner, 1, this%el_owner, 1, &
                                this%rst, this%gdim, &
                                this%n_points, field)
    else
       if (this%n_points .gt. 0) &
          call this%local_interp%evaluate(interp_values, this%el_owner, &
                                          field, this%nelv)
    end if
#else
    call neko_error('Neko needs to be built with GSLIB support')
#endif
 
  end subroutine global_interpolation_evaluate
 
end module global_interpolation