global_interpolation.f90

Go to the documentation of this file.

! Copyright (c) 2020-2025, 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, dp, xp
  use neko_config, only : neko_bcknd_device
  use space, only : space_t
  use stack, only : stack_i4_t
  use dofmap, only : dofmap_t
  use logger, only : neko_log, log_size
  use json_utils, only : json_get_or_lookup_or_default
  use json_module, only : json_file
  use utils, only : neko_error
  use local_interpolation, only : local_interpolator_t
  use device, only : device_free, device_map, device_memcpy, &
       device_unmap, host_to_device, device_to_host, &
       device_get_ptr
  use aabb_pe_finder, only : aabb_pe_finder_t
  use aabb_el_finder, only : aabb_el_finder_t
  use cartesian_el_finder, only : cartesian_el_finder_t
  use cartesian_pe_finder, only : cartesian_pe_finder_t
  use legendre_rst_finder, only : legendre_rst_finder_t
  use el_finder, only : el_finder_t
  use pe_finder, only : pe_finder_t
  use comm, only : neko_comm
  use mpi_f08, only : mpi_sum, mpi_comm, mpi_comm_rank, &
       mpi_comm_size, mpi_wtime, mpi_allreduce, mpi_in_place, mpi_integer, &
       mpi_min, mpi_barrier, mpi_reduce_scatter_block, mpi_alltoall, &
       mpi_isend, mpi_irecv
  use glb_intrp_comm, only : glb_intrp_comm_t
  use vector, only : vector_t
  use vector_math, only : vector_masked_gather_copy
  use matrix, only : matrix_t
  use math, only : copy, neko_eps
  use mask, only : mask_t
  use structs, only : array_ptr_t
  use, intrinsic :: iso_c_binding, only : c_ptr, c_null_ptr, c_associated
  implicit none
  private
 
  integer, public, parameter :: glob_map_size = 4096
  real(kind=dp), public, parameter :: glob_interp_tol = neko_eps*1e3_dp
  real(kind=dp), public, parameter :: glob_interp_pad = 1e-2_dp
 
  type, public :: global_interpolation_settings_t
     integer :: glob_map_size = glob_map_size
     real(kind=dp) :: tolerance = glob_interp_tol
     real(kind=dp) :: padding = glob_interp_pad
  end type global_interpolation_settings_t
 
  type, public :: global_interpolation_t
     type(vector_t) :: x
     type(vector_t) :: y
     type(vector_t) :: z
     integer :: gdim
     integer :: nelv
     integer :: glb_nelv
     type(mpi_comm) :: comm
     integer :: pe_rank
     integer :: pe_size
     type(space_t) :: xh
     integer :: n_points
     real(kind=rp), allocatable :: xyz(:,:)
     real(kind=rp), allocatable :: rst(:,:)
     integer, allocatable :: pe_owner(:)
     type(stack_i4_t), allocatable :: points_at_pe(:)
     integer, allocatable :: el_owner0(:)
     type(c_ptr) :: el_owner0_d = c_null_ptr
 
     integer :: n_points_local
     integer, allocatable :: el_owner0_local(:)
     type(c_ptr) :: el_owner0_local_d = c_null_ptr
     real(kind=rp), allocatable :: rst_local(:,:)
     real(kind=rp), allocatable :: xyz_local(:,:)
     type(local_interpolator_t) :: local_interp
     logical :: all_points_local = .false.
     real(kind=dp) :: tolerance = glob_interp_tol
     real(kind=dp) :: padding = glob_interp_pad
 
     integer, allocatable :: n_points_pe(:)
     integer, allocatable :: n_points_offset_pe(:)
     integer, allocatable :: n_points_pe_local(:)
     integer, allocatable :: n_points_offset_pe_local(:)
     class(pe_finder_t), allocatable :: pe_finder
     class(el_finder_t), allocatable :: el_finder
     type(legendre_rst_finder_t) :: rst_finder
     type(glb_intrp_comm_t) :: glb_intrp_comm
     type(vector_t) :: temp_local, temp
     integer :: n_dof = -1
     type(vector_t) :: masked_field
 
   contains
     procedure, pass(this) :: init_json_xyz => &
          global_interpolation_init_json_xyz
     procedure, pass(this) :: init_json_dof => &
          global_interpolation_init_json_dof
     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) :: free_points_local => &
          global_interpolation_free_points_local
     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_coords1d => &
          global_interpolation_find_coords1d
     procedure, pass(this) :: check_points => &
          global_interpolation_check_points
     procedure, pass(this) :: find_points_xyz => global_interpolation_find_xyz
     generic :: find_points => find_points_xyz, find_points_coords, &
          find_points_coords1d
     procedure, pass(this) :: evaluate => global_interpolation_evaluate
     procedure, pass(this) :: evaluate_masked => &
          global_interpolation_evaluate_masked
     procedure, pass(this) :: init_redist_comm => &
          global_interpolation_init_redist_comm
 
     generic :: init => init_dof, init_xyz, init_json_xyz, init_json_dof
 
  end type global_interpolation_t
 
contains
 
  subroutine global_interpolation_init_json_xyz(this, x, y, z, gdim, nelv, Xh, &
       params_subdict, comm)
    class(global_interpolation_t), target, intent(inout) :: this
    real(kind=rp), intent(in) :: x(:)
    real(kind=rp), intent(in) :: y(:)
    real(kind=rp), intent(in) :: z(:)
    integer, intent(in) :: gdim
    integer, intent(in) :: nelv
    type(space_t), intent(in) :: Xh
    type(json_file), intent(inout) :: params_subdict
    type(mpi_comm), intent(in), optional :: comm
 
    real(kind=dp) :: tol, pad
 
    call json_get_or_lookup_or_default(params_subdict, 'tolerance', &
         tol, glob_interp_tol)
    call json_get_or_lookup_or_default(params_subdict, 'padding', &
         pad, glob_interp_pad)
 
    call this%init_xyz(x, y, z, gdim, nelv, xh, comm = comm, tol = tol, &
         pad = pad)
 
  end subroutine global_interpolation_init_json_xyz
 
  subroutine global_interpolation_init_json_dof(this, dof, params_subdict, &
       comm, mask)
    class(global_interpolation_t), target, intent(inout) :: this
    type(dofmap_t) :: dof
    type(json_file), intent(inout) :: params_subdict
    type(mpi_comm), optional, intent(in) :: comm
    type(mask_t), intent(in), optional :: mask
 
    real(kind=dp) :: tol, pad
 
    call json_get_or_lookup_or_default(params_subdict, 'tolerance', &
         tol, glob_interp_tol)
    call json_get_or_lookup_or_default(params_subdict, 'padding', &
         pad, glob_interp_pad)
 
    call this%init_dof(dof, comm = comm, tol = tol, pad = pad, mask = mask)
 
  end subroutine global_interpolation_init_json_dof
 
  subroutine global_interpolation_init_dof(this, dof, comm, tol, pad, mask)
    class(global_interpolation_t), target, intent(inout) :: this
    type(dofmap_t) :: dof
    type(mpi_comm), optional, intent(in) :: comm
    real(kind=dp), optional :: tol
    real(kind=dp), optional :: pad
    type(mask_t), intent(in), optional :: mask
 
    integer :: temp_nelv
 
    ! Store the number of dofs
    this%n_dof = dof%size()
    ! NOTE: Passing dof%x%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
    ! to get the right dimension (see global_interpolation_init_xyz).
    if (.not. present(mask)) then
       call this%init_xyz(dof%x%x(:,1,1,1), dof%y%x(:,1,1,1), dof%z%x(:,1,1,1), &
            dof%msh%gdim, dof%msh%nelv, dof%Xh, comm = comm, &
            tol = tol, pad = pad)
    else
 
       ! Verify that the mask size is compatible with the dofmap
       temp_nelv = mask%size() / (dof%Xh%lx*dof%Xh%ly*dof%Xh%lz)
       if (mod(mask%size(), dof%Xh%lx*dof%Xh%ly*dof%Xh%lz) /= 0) then
          call neko_error("Mask size must be a multiple of the number of" // &
               " elements in the mesh.")
       end if
       ! Initialize with the masked coordinates
       call this%init_xyz(dof%x%x(mask%get(),1,1,1), dof%y%x(mask%get(),1,1,1), &
            dof%z%x(mask%get(),1,1,1), dof%msh%gdim, temp_nelv, dof%Xh, &
            comm = comm, tol = tol, pad = pad)
       ! Initialize a helper field with the size of the mask, after
       ! init_xyz since that frees the object
       call this%masked_field%init(mask%size())
    end if
 
  end subroutine global_interpolation_init_dof
 
  subroutine global_interpolation_init_xyz(this, x, y, z, gdim, nelv, Xh, &
       comm, tol, pad)
    class(global_interpolation_t), target, intent(inout) :: this
    real(kind=rp), intent(in) :: x(:)
    real(kind=rp), intent(in) :: y(:)
    real(kind=rp), intent(in) :: z(:)
    integer, intent(in) :: gdim
    integer, intent(in) :: nelv
    type(space_t), intent(in) :: Xh
    type(mpi_comm), intent(in), optional :: comm
    real(kind=dp), intent(in), optional :: tol
    real(kind=dp), intent(in), optional :: pad
 
    integer :: lx, ly, lz, ierr, i, n
    character(len=8000) :: log_buf
    !real(kind=dp) :: padding ! <-- note that this padding is in dp
    real(kind=rp) :: time1, time_start
    character(len=255) :: mode_str
    integer :: boxdim, envvar_len
 
    call neko_log%section('Global Interpolation')
    call neko_log%message('Initializing global interpolation')
 
    call this%free()
 
    ! Set communicator
    if (present(comm)) then
       this%comm = comm
    else
       this%comm = neko_comm
    end if
 
    ! Set point search parameters
    this%padding = glob_interp_pad
    if (present(pad)) this%padding = pad
 
    this%tolerance = glob_interp_tol
    if (present(tol)) this%tolerance = tol
 
    write(log_buf, '(A,E15.7)') &
         'Tolerance: ', this%tolerance
    call neko_log%message(log_buf)
    write(log_buf, '(A,E15.7)') &
         'Padding  : ', this%padding
    call neko_log%message(log_buf)
 
    time_start = mpi_wtime()
    call mpi_barrier(this%comm)
 
    call mpi_comm_rank(this%comm, this%pe_rank, ierr)
    call mpi_comm_size(this%comm, this%pe_size, ierr)
 
    this%gdim = gdim
    this%nelv = nelv
 
    call mpi_allreduce(nelv, this%glb_nelv, 1, mpi_integer, &
         mpi_sum, this%comm, ierr)
    lx = xh%lx
    ly = xh%ly
    lz = xh%lz
    n = nelv * lx*ly*lz
    call this%x%init(n)
    call this%y%init(n)
    call this%z%init(n)
    call copy(this%x%x, x, n)
    call this%x%copy_from(host_to_device,.false.)
    call copy(this%y%x, y, n)
    call this%y%copy_from(host_to_device,.false.)
    call copy(this%z%x, z, n)
    call this%z%copy_from(host_to_device,.false.)
    call this%Xh%init(xh%t, lx, ly, lz)
 
    ! Initialize n_dof if it was not started with a dof-based constructor
    if (this%n_dof == -1) then
       this%n_dof = n
    end if
 
 
    call get_environment_variable("NEKO_GLOBAL_INTERP_EL_FINDER", &
         mode_str, envvar_len)
 
    if (envvar_len .gt. 0) then
       if (mode_str(1:envvar_len) == 'AABB') then
          allocate(aabb_el_finder_t :: this%el_finder)
       end if
    end if
    call get_environment_variable("NEKO_GLOBAL_INTERP_PE_FINDER", &
         mode_str, envvar_len)
 
    if (envvar_len .gt. 0) then
       if (mode_str(1:envvar_len) == 'AABB') then
          allocate(aabb_pe_finder_t :: this%pe_finder)
       end if
    end if
 
    if (.not. allocated(this%el_finder)) then
       allocate(cartesian_el_finder_t :: this%el_finder)
    end if
    if (.not. allocated(this%pe_finder)) then
       allocate(cartesian_pe_finder_t :: this%pe_finder)
    end if
    select type (el_find => this%el_finder)
    type is (aabb_el_finder_t)
       call neko_log%message('Using AABB element finder')
       call el_find%init(x, y, z, nelv, xh, this%padding)
    type is (cartesian_el_finder_t)
       call neko_log%message('Using Cartesian element finder')
       boxdim = max(lx*int(real(nelv, xp)**(1.0_xp / 3.0_xp)), 2)
       boxdim = min(boxdim, 300)
       call el_find%init(x, y, z, nelv, xh, boxdim, this%padding)
    class default
       call neko_error('Unknown element finder type')
    end select
 
    select type (pe_find => this%pe_finder)
    type is (aabb_pe_finder_t)
       call neko_log%message('Using AABB PE finder')
       call pe_find%init(this%x%x, this%y%x, this%z%x, &
            nelv, xh, this%comm, this%padding)
    type is (cartesian_pe_finder_t)
       call neko_log%message('Using Cartesian PE finder')
       boxdim = lx*int(real(this%glb_nelv, xp)**(1.0_xp / 3.0_xp))
       boxdim = max(boxdim, 32)
       boxdim = min(boxdim, &
            int(8.0_xp*(30000.0_xp * this%pe_size)**(1.0_xp / 3.0_xp)))
       call pe_find%init(this%x%x, this%y%x, this%z%x, &
            nelv, xh, this%comm, boxdim, this%padding)
    class default
       call neko_error('Unknown PE finder type')
    end select
 
    call this%rst_finder%init(this%x%x, this%y%x, this%z%x, nelv, xh, &
         this%tolerance)
    if (allocated(this%n_points_pe)) deallocate(this%n_points_pe)
    if (allocated(this%n_points_pe_local)) deallocate(this%n_points_pe_local)
    if (allocated(this%n_points_offset_pe_local)) &
         deallocate(this%n_points_offset_pe_local)
    if (allocated(this%n_points_offset_pe)) deallocate(this%n_points_offset_pe)
    allocate(this%n_points_pe(0:(this%pe_size-1)))
    allocate(this%n_points_offset_pe(0:(this%pe_size-1)))
    allocate(this%n_points_pe_local(0:(this%pe_size-1)))
    allocate(this%n_points_offset_pe_local(0:(this%pe_size-1)))
    allocate(this%points_at_pe(0:(this%pe_size-1)))
    do i = 0, this%pe_size-1
       call this%points_at_pe(i)%init()
    end do
    call mpi_barrier(this%comm)
    time1 = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Global interpolation initialized (s):', time1-time_start
    call neko_log%message(log_buf)
    call neko_log%end_section()
  end subroutine global_interpolation_init_xyz
 
 
  subroutine global_interpolation_free(this)
    class(global_interpolation_t), target, intent(inout) :: this
    integer :: i
 
    call this%x%free()
    call this%y%free()
    call this%z%free()
    call this%Xh%free()
 
    this%nelv = 0
    this%gdim = 0
 
    call this%free_points()
    call this%free_points_local()
    call this%local_interp%free()
    if (allocated(this%el_finder)) then
       call this%el_finder%free()
       deallocate(this%el_finder)
    end if
    if (allocated(this%pe_finder)) then
       call this%pe_finder%free()
       deallocate(this%pe_finder)
    end if
    call this%rst_finder%free()
 
    call this%temp_local%free()
    call this%temp%free()
    call this%masked_field%free()
    if (allocated(this%points_at_pe)) then
       do i = 0, this%pe_size-1
          call this%points_at_pe(i)%free()
       end do
       deallocate(this%points_at_pe)
    end if
    if (allocated(this%n_points_pe)) deallocate(this%n_points_pe)
    if (allocated(this%n_points_pe_local)) deallocate(this%n_points_pe_local)
    if (allocated(this%n_points_offset_pe_local)) &
         deallocate(this%n_points_offset_pe_local)
    if (allocated(this%n_points_offset_pe)) deallocate(this%n_points_offset_pe)
 
 
 
  end subroutine global_interpolation_free
 
  subroutine global_interpolation_free_points(this)
    class(global_interpolation_t), target, 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%pe_owner)) deallocate(this%pe_owner)
    if (allocated(this%el_owner0)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(this%el_owner0, this%el_owner0_d)
       end if
       deallocate(this%el_owner0)
    end if
 
    call this%glb_intrp_comm%free()
 
 
  end subroutine global_interpolation_free_points
 
 
  subroutine global_interpolation_free_points_local(this)
    class(global_interpolation_t), target, intent(inout) :: this
 
    this%n_points_local = 0
    this%all_points_local = .false.
 
    if (allocated(this%xyz_local)) deallocate(this%xyz_local)
    if (allocated(this%rst_local)) deallocate(this%rst_local)
    if (allocated(this%el_owner0_local)) then
       if (neko_bcknd_device .eq. 1) then
          call device_unmap(this%el_owner0_local, this%el_owner0_local_d)
       end if
       deallocate(this%el_owner0_local)
    end if
 
  end subroutine global_interpolation_free_points_local
 
 
  subroutine global_interpolation_find_common(this)
    class(global_interpolation_t), target, intent(inout) :: this
    character(len=8000) :: log_buf
    type(vector_t) :: x_t
    type(vector_t) :: y_t
    type(vector_t) :: z_t
    type(matrix_t) :: rst_local_cand
    type(vector_t) :: resx
    type(vector_t) :: resy
    type(vector_t) :: resz
    type(c_ptr) :: el_cands_d
    type(matrix_t) :: res
    integer :: i, j, stupid_intent
    type(stack_i4_t), target :: all_el_candidates
    integer, allocatable :: n_el_cands(:)
    integer, contiguous, pointer :: el_cands(:), point_ids(:), send_recv(:)
    real(kind=rp), allocatable :: res_results(:,:)
    real(kind=rp), allocatable :: rst_results(:,:)
    integer, allocatable :: el_owner_results(:)
    integer :: ierr, ii, n_point_cand, n_glb_point_cand, point_id, rank
    real(kind=rp) :: time1, time2, time_start
    !Temp stuff for glb_intrp_comm
    type(stack_i4_t) :: send_pe, recv_pe
    type(glb_intrp_comm_t) :: glb_intrp_find, glb_intrp_find_back
    type(stack_i4_t) :: send_pe_find, recv_pe_find
 
    el_cands_d = c_null_ptr
    call neko_log%section('Global Interpolation')
    call glb_intrp_find%init_dofs(this%pe_size)
    call send_pe_find%init()
    call recv_pe_find%init()
    call mpi_barrier(this%comm)
    time_start = mpi_wtime()
    write(log_buf, '(A)') 'Global interpolation, finding points'
    call neko_log%message(log_buf)
    ! Find pe candidates that the points i want may be at
    ! Add number to n_points_pe_local
    !Working arrays
    this%n_points_pe = 0
    call this%pe_finder%find_batch(this%xyz, this%n_points, &
         this%points_at_pe, this%n_points_pe)
    call mpi_barrier(this%comm)
    time1 = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Found PE candidates time since start of findpts (s):', &
         time1 - time_start
    call neko_log%message(log_buf)
 
    !Send number of points I want to candidates
    ! n_points_local -> how many points might be at this rank
    ! n_points_pe_local -> how many points local on this rank that other pes
    !                      might want
    this%n_points_pe_local = 0
    this%n_points_local = 0
    call mpi_reduce_scatter_block(this%n_points_pe, this%n_points_local, &
         1, mpi_integer, mpi_sum, this%comm, ierr)
    call mpi_alltoall(this%n_points_pe, 1, mpi_integer,&
         this%n_points_pe_local, 1, mpi_integer, this%comm, ierr)
 
    !Set up offset arrays
    this%n_points_offset_pe_local(0) = 0
    this%n_points_offset_pe(0) = 0
    do i = 1, (this%pe_size - 1)
       this%n_points_offset_pe_local(i) = this%n_points_pe_local(i-1)&
            + this%n_points_offset_pe_local(i-1)
       this%n_points_offset_pe(i) = this%n_points_pe(i-1)&
            + this%n_points_offset_pe(i-1)
    end do
    do i = 0, (this%pe_size-1)
       if (this%n_points_pe(i) .gt. 0) then
          call send_pe_find%push(i)
          point_ids => this%points_at_pe(i)%array()
          do j = 1, this%n_points_pe(i)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+1)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+2)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j)-1)+3)
          end do
       end if
       if (this%n_points_pe_local(i) .gt. 0) then
          call recv_pe_find%push(i)
          do j = 1, this%n_points_pe_local(i)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 1)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 2)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 3)
          end do
       end if
    end do
 
 
 
    call glb_intrp_find%init(send_pe_find, recv_pe_find, this%comm)
 
    call glb_intrp_find_back%init_dofs(this%pe_size)
    ii = 0
    do i = 0, (this%pe_size-1)
       send_recv => glb_intrp_find%recv_dof(i)%array()
       do j = 1, glb_intrp_find%recv_dof(i)%size()
          call glb_intrp_find_back%send_dof(i)%push(send_recv(j))
       end do
       send_recv => glb_intrp_find%send_dof(i)%array()
       do j = 1, glb_intrp_find%send_dof(i)%size()
          ii = ii + 1
          call glb_intrp_find_back%recv_dof(i)%push(ii)
       end do
    end do
 
    call glb_intrp_find_back%init(recv_pe_find, send_pe_find, this%comm)
 
 
    if (allocated(this%xyz_local)) then
       deallocate(this%xyz_local)
    end if
    allocate(this%xyz_local(3, this%n_points_local))
    call glb_intrp_find%sendrecv(this%xyz, this%xyz_local, this%n_points*3, &
         this%n_points_local*3)
 
    call mpi_barrier(this%comm)
    time1 = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Sent to points to PE candidates, time since start of ' &
         // 'find_points (s):', time1 - time_start
    call neko_log%message(log_buf)
 
    !Okay, now we need to find the rst...
    call all_el_candidates%init()
 
    if (allocated(n_el_cands)) then
       deallocate(n_el_cands)
    end if
 
    allocate(n_el_cands(this%n_points_local))
    call this%el_finder%find_batch(this%xyz_local, this%n_points_local, &
         all_el_candidates, n_el_cands)
 
    n_point_cand = all_el_candidates%size()
    if (n_point_cand .gt. 1e8) then
       print *, 'Warning, many point candidates on rank', this%pe_rank, &
            'cands:', n_point_cand, &
            'Consider increasing number of ranks'
    end if
    call x_t%init(n_point_cand)
    call y_t%init(n_point_cand)
    call z_t%init(n_point_cand)
    ii = 0
    do i = 1 , this%n_points_local
       do j = 1, n_el_cands(i)
          ii = ii + 1
          x_t%x(ii) = this%xyz_local(1,i)
          y_t%x(ii) = this%xyz_local(2,i)
          z_t%x(ii) = this%xyz_local(3,i)
       end do
    end do
 
    call mpi_barrier(this%comm)
    time1 = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Element candidates found, now time for finding rst, time ' // &
         'since start of find_points (s):', time1 - time_start
    call neko_log%message(log_buf)
    call rst_local_cand%init(3, n_point_cand)
    call resx%init(n_point_cand)
    call resy%init(n_point_cand)
    call resz%init(n_point_cand)
 
    ! Find rst within all element candidates for target xyz (x_t, y_t, z_t)
    call mpi_barrier(this%comm)
    time1 = mpi_wtime()
    el_cands => all_el_candidates%array()
    if (neko_bcknd_device .eq. 1) then
       call x_t%copy_from(host_to_device,.false.)
 
       call y_t%copy_from(host_to_device,.false.)
 
       call z_t%copy_from(host_to_device,.false.)
       call device_map(el_cands, el_cands_d, n_point_cand)
       call device_memcpy(el_cands, el_cands_d, n_point_cand, &
            host_to_device, .true.)
    end if
 
    call this%rst_finder%find(rst_local_cand, &
         x_t, y_t, z_t, &
         el_cands, n_point_cand, &
         resx, resy, resz)
    if (neko_bcknd_device .eq. 1) then
       call rst_local_cand%copy_from(device_to_host,.false.)
       call resx%copy_from(device_to_host,.false.)
       call resy%copy_from(device_to_host,.false.)
       call resz%copy_from(device_to_host,.true.)
       call device_unmap(el_cands, el_cands_d)
    end if
    call mpi_barrier(this%comm)
 
    time2 = mpi_wtime()
    write(log_buf, '(A,E15.7)') &
         'Found rst with Newton iteration, time (s):', time2-time1
    call neko_log%message(log_buf)
 
    write(log_buf, '(A)') &
         'Checking validity of points and choosing best candidates.'
    call neko_log%message(log_buf)
    call mpi_barrier(this%comm, ierr)
 
    if (allocated(this%rst_local)) deallocate(this%rst_local)
    if (allocated(this%el_owner0_local)) deallocate(this%el_owner0_local)
    allocate(this%rst_local(3, this%n_points_local))
    allocate(this%el_owner0_local(this%n_points_local))
    ! Choose the best candidate at this rank
    ii = 0
    do i = 1 , this%n_points_local
       this%xyz_local(1,i) = 10.0
       this%xyz_local(2,i) = 10.0
       this%xyz_local(3,i) = 10.0
       this%rst_local(1,i) = 10.0
       this%rst_local(2,i) = 10.0
       this%rst_local(3,i) = 10.0
       this%el_owner0_local(i) = -1
       do j = 1, n_el_cands(i)
          ii = ii + 1
          if (rst_cmp(this%rst_local(:, i), rst_local_cand%x(:, ii), &
               this%xyz_local(:, i), [resx%x(ii), resy%x(ii), resz%x(ii)], &
               this%padding)) then
             this%rst_local(1, i) = rst_local_cand%x(1, ii)
             this%rst_local(2, i) = rst_local_cand%x(2, ii)
 
             this%rst_local(3, i) = rst_local_cand%x(3, ii)
             this%xyz_local(1,i) = resx%x(ii)
             this%xyz_local(2,i) = resy%x(ii)
             this%xyz_local(3,i) = resz%x(ii)
             this%el_owner0_local(i) = el_cands(ii)
          end if
          ! if (this%pe_rank .eq. 0) print *,i,  this%rst_local(:,i), &
          !  this%xyz_local(:,i), this%el_owner0_local(i)
       end do
    end do
    call res%init(3, this%n_points)
    n_glb_point_cand = sum(this%n_points_pe)
    if (allocated(rst_results)) deallocate(rst_results)
    if (allocated(res_results)) deallocate(res_results)
    if (allocated(el_owner_results)) deallocate(el_owner_results)
    allocate(rst_results(3, n_glb_point_cand))
    allocate(res_results(3, n_glb_point_cand))
    allocate(el_owner_results(n_glb_point_cand))
    res = 1e2_rp
    this%rst = 1e2
    this%pe_owner = -1
    this%el_owner0 = -1
    call glb_intrp_find_back%sendrecv(this%xyz_local, res_results, &
         this%n_points_local*3, n_glb_point_cand*3)
    call glb_intrp_find_back%sendrecv(this%rst_local, rst_results, &
         this%n_points_local*3, n_glb_point_cand*3)
    do i = 1, size(glb_intrp_find_back%send_pe)
       rank = glb_intrp_find_back%send_pe(i)
       call mpi_isend(this%el_owner0_local( &
            this%n_points_offset_pe_local(rank) + 1), &
            this%n_points_pe_local(rank), &
            mpi_integer, rank, 0, &
            this%comm, glb_intrp_find_back%send_buf(i)%request, ierr)
       glb_intrp_find_back%send_buf(i)%flag = .false.
    end do
    do i = 1, size(glb_intrp_find_back%recv_pe)
       rank = glb_intrp_find_back%recv_pe(i)
       call mpi_irecv(el_owner_results(this%n_points_offset_pe(rank)+1),&
            this%n_points_pe(rank), &
            mpi_integer, rank, 0, &
            this%comm, glb_intrp_find_back%recv_buf(i)%request, ierr)
       glb_intrp_find_back%recv_buf(i)%flag = .false.
    end do
    call glb_intrp_find_back%nbwait_no_op()
    ii = 0
    do i = 1, size(glb_intrp_find_back%recv_pe)
       point_ids => this%points_at_pe(glb_intrp_find_back%recv_pe(i))%array()
       do j = 1, this%n_points_pe(glb_intrp_find_back%recv_pe(i))
          point_id = point_ids(j)
          ii = ii + 1
          if (rst_cmp(this%rst(:, point_id), rst_results(:, ii), &
               res%x(:, point_id), res_results(:, ii), this%padding) .or. &
               this%pe_owner(point_ids(j)) .eq. -1 ) then
             this%rst(:, point_ids(j)) = rst_results(:, ii)
             res%x(:, point_ids(j)) = res_results(:, ii)
             this%pe_owner(point_ids(j)) = glb_intrp_find_back%recv_pe(i)
             this%el_owner0(point_ids(j)) = el_owner_results(ii)
          end if
          !  if (this%pe_rank .eq. 0) print *,point_id,  &
          !this%rst(:,point_ids(j)),res%x(:,point_ids(j)),
          ! this%el_owner0(point_ids(j))
       end do
    end do
 
    !OK, now I know the correct rst values of the points I want We now
    !send the correct rsts to the correct rank (so a point only
    !belongs to one rank)
    do i = 0, this%pe_size-1
       call this%points_at_pe(i)%clear()
       this%n_points_pe(i) = 0
    end do
 
    do i = 1, this%n_points
       stupid_intent = i
       if (this%pe_owner(i) .eq. -1 .or. this%el_owner0(i) .eq. -1) then
          print *, 'No owning rank found for',&
               ' point ', stupid_intent, ' with coords', this%xyz(:,i), &
               ' Interpolation will always yield 0.0. Try increase padding.'
       else
          call this%points_at_pe(this%pe_owner(i))%push(stupid_intent)
          this%n_points_pe(this%pe_owner(i)) = &
               this%n_points_pe(this%pe_owner(i)) + 1
       end if
    end do
    call mpi_reduce_scatter_block(this%n_points_pe, this%n_points_local, 1, &
         mpi_integer, mpi_sum, this%comm, ierr)
    call mpi_alltoall(this%n_points_pe, 1, mpi_integer, &
         this%n_points_pe_local, 1, mpi_integer, this%comm, ierr)
    this%n_points_offset_pe_local(0) = 0
    this%n_points_offset_pe(0) = 0
    do i = 1, (this%pe_size - 1)
       this%n_points_offset_pe_local(i) = this%n_points_pe_local(i-1)&
            + this%n_points_offset_pe_local(i-1)
       this%n_points_offset_pe(i) = this%n_points_pe(i-1)&
            + this%n_points_offset_pe(i-1)
    end do
    call send_pe_find%free()
    call recv_pe_find%free()
    call glb_intrp_find%free()
    call send_pe_find%init()
    call recv_pe_find%init()
    call glb_intrp_find%init_dofs(this%pe_size)
    !setup comm to send xyz and rst to chosen ranks
    do i = 0, (this%pe_size-1)
       if (this%n_points_pe(i) .gt. 0) then
          call send_pe_find%push(i)
          point_ids => this%points_at_pe(i)%array()
          do j = 1, this%n_points_pe(i)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j) - 1) + 1)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j) - 1) + 2)
             call glb_intrp_find%send_dof(i)%push(3*(point_ids(j) - 1) + 3)
          end do
       end if
       if (this%n_points_pe_local(i) .gt. 0) then
          call recv_pe_find%push(i)
          do j = 1, this%n_points_pe_local(i)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 1)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 2)
             call glb_intrp_find%recv_dof(i)%push(3*(j + &
                  this%n_points_offset_pe_local(i) - 1) + 3)
          end do
       end if
    end do
 
 
    call glb_intrp_find%init(send_pe_find, recv_pe_find, this%comm)
    call glb_intrp_find%sendrecv(this%xyz, this%xyz_local, this%n_points*3, &
         this%n_points_local*3)
    call glb_intrp_find%sendrecv(this%rst, this%rst_local, this%n_points*3, &
         this%n_points_local*3)
    ii = 0
    do i = 1, size(glb_intrp_find%send_pe)
       rank = glb_intrp_find%send_pe(i)
       point_ids => this%points_at_pe(rank)%array()
       do j = 1, this%n_points_pe(rank)
          ii = ii + 1
          el_owner_results(ii) = this%el_owner0(point_ids(j))
       end do
       call mpi_isend(el_owner_results(this%n_points_offset_pe(rank) + 1),&
            this%n_points_pe(rank), &
            mpi_integer, rank, 0, &
            this%comm, glb_intrp_find%send_buf(i)%request, ierr)
       glb_intrp_find%send_buf(i)%flag = .false.
    end do
    do i = 1, size(glb_intrp_find%recv_pe)
       rank = glb_intrp_find%recv_pe(i)
       call mpi_irecv(this%el_owner0_local( &
            this%n_points_offset_pe_local(rank) + 1), &
            this%n_points_pe_local(rank), &
            mpi_integer, rank, 0, &
            this%comm, glb_intrp_find%recv_buf(i)%request, ierr)
       glb_intrp_find%recv_buf(i)%flag = .false.
    end do
    call glb_intrp_find%nbwait_no_op()
 
    call glb_intrp_find%free()
 
    !Set up final way of doing communication
    call send_pe%init()
    call recv_pe%init()
    call this%glb_intrp_comm%init_dofs(this%pe_size)
    do i = 0, (this%pe_size-1)
       if (this%n_points_pe(i) .gt. 0) then
          call recv_pe%push(i)
          point_ids => this%points_at_pe(i)%array()
          do j = 1, this%n_points_pe(i)
             call this%glb_intrp_comm%recv_dof(i)%push(point_ids(j))
          end do
       end if
       if (this%n_points_pe_local(i) .gt. 0) then
          call send_pe%push(i)
          do j = 1, this%n_points_pe_local(i)
             call this%glb_intrp_comm%send_dof(i)%push(j + &
                  this%n_points_offset_pe_local(i))
          end do
       end if
    end do
    call this%glb_intrp_comm%init(send_pe, recv_pe, this%comm)
 
    !Initialize working arrays for evaluation
    call this%temp_local%init(this%n_points_local)
    call this%temp%init(this%n_points)
 
    !Initialize interpolator for local interpolation
    call this%local_interp%init(this%Xh, this%rst_local, &
         this%n_points_local)
 
 
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%el_owner0, this%el_owner0_d, &
            this%n_points, host_to_device, sync = .true.)
       call device_map(this%el_owner0_local, this%el_owner0_local_d, &
            this%n_points_local)
       call device_memcpy(this%el_owner0_local, this%el_owner0_local_d, &
            this%n_points_local, host_to_device, sync = .true.)
    end if
 
    call this%check_points(this%x%x, this%y%x, this%z%x)
 
    !Free stuff
    call send_pe%free()
    call recv_pe%free()
    call glb_intrp_find_back%free()
    call send_pe_find%free()
    call recv_pe_find%free()
    call x_t%free()
    call y_t%free()
    call z_t%free()
    call rst_local_cand%free()
    call resx%free()
    call resy%free()
    call resz%free()
    call res%free()
    call all_el_candidates%free()
 
    if (allocated(n_el_cands)) deallocate(n_el_cands)
    if (allocated(rst_results)) deallocate(rst_results)
    if (allocated(res_results)) deallocate(res_results)
    if (allocated(el_owner_results)) deallocate(el_owner_results)
    call mpi_barrier(this%comm, ierr)
    time2 = mpi_wtime()
    write(log_buf, '(A,E15.7)') 'Global interpolation find points ' // &
         'done, time (s):', time2-time_start
    call neko_log%message(log_buf)
    call neko_log%end_section()
    call neko_log%newline()
  end subroutine global_interpolation_find_common
 
  subroutine global_interpolation_check_points(this, x, y, z)
    class(global_interpolation_t), target, intent(inout) :: this
    real(kind=rp), intent(inout) :: x(:)
    real(kind=rp), intent(inout) :: y(:)
    real(kind=rp), intent(inout) :: z(:)
    integer :: i, j
    character(len=8000) :: log_buf
    real(kind=rp) :: xdiff, ydiff, zdiff
    logical :: isdiff
    type(vector_t) :: x_check, y_check, z_check
 
    call x_check%init(this%n_points)
    call y_check%init(this%n_points)
    call z_check%init(this%n_points)
    call this%evaluate(x_check%x, x, on_host = .true.)
    call this%evaluate(y_check%x, y, on_host = .true.)
    call this%evaluate(z_check%x, z, on_host = .true.)
    write(log_buf, '(A)') 'Checking validity of points.'
    call neko_log%message(log_buf)
    j = 0
    do i = 1 , this%n_points
       ! Check validity of points
       isdiff = .false.
       xdiff = x_check%x(i)-this%xyz(1,i)
       ydiff = y_check%x(i)-this%xyz(2,i)
       zdiff = z_check%x(i)-this%xyz(3,i)
       isdiff = norm2(real([xdiff, ydiff, zdiff], xp)) > this%tolerance
       if (isdiff) then
          write(*, *) 'Point ', i, 'at rank ', this%pe_rank, &
               'with coordinates: ', &
               this%xyz(1, i), this%xyz(2, i), this%xyz(3, i), &
               'Differ from interpolated coords: ', &
               x_check%x(i), y_check%x(i), z_check%x(i), &
               'Actual difference: ', &
               xdiff, ydiff, zdiff, norm2(real([xdiff, ydiff, zdiff], xp)), &
               'Process, element: ', &
               this%pe_owner(i), this%el_owner0(i)+1, &
               'Calculated rst: ', &
               this%rst(1,i), this%rst(2,i), this%rst(3,i)
          j = j + 1
       end if
    end do
    call x_check%free()
    call y_check%free()
    call z_check%free()
 
  end subroutine global_interpolation_check_points
 
  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()
    call this%free_points_local()
 
    this%n_points = n_points
 
    call global_interpolation_init_point_arrays(this)
 
    !Deepcopy of coordinates
    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_find_coords1d(this, x, y, z, n_points)
    class(global_interpolation_t), intent(inout) :: this
    integer :: n_points
    real(kind=rp) :: x(n_points)
    real(kind=rp) :: y(n_points)
    real(kind=rp) :: z(n_points)
    integer :: i
 
    call this%free_points()
    call this%free_points_local()
 
    this%n_points = n_points
 
    call global_interpolation_init_point_arrays(this)
 
    !Deepcopy of coordinates
    do i = 1, n_points
       this%xyz(1, i) = x(i)
       this%xyz(2, i) = y(i)
       this%xyz(3, i) = z(i)
    end do
 
    call global_interpolation_find_common(this)
 
  end subroutine global_interpolation_find_coords1d
 
 
  subroutine global_interpolation_init_point_arrays(this)
    class(global_interpolation_t) :: this
 
    if (allocated(this%xyz)) deallocate(this%xyz)
    if (allocated(this%rst)) deallocate(this%rst)
    if (allocated(this%pe_owner)) deallocate(this%pe_owner)
    if (allocated(this%el_owner0)) deallocate(this%el_owner0)
 
    allocate(this%pe_owner(this%n_points))
    allocate(this%el_owner0(this%n_points))
    allocate(this%xyz(3, this%n_points))
    allocate(this%rst(3, this%n_points))
    if (neko_bcknd_device .eq. 1) then
       call device_map(this%el_owner0, this%el_owner0_d, this%n_points)
    end if
 
  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()
    call this%free_points_local()
 
    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(:,:)
 
 
    call this%free_points()
    call this%free_points_local()
 
 
    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 this%free_points()
    this%n_points = this%n_points_local
    n_points = this%n_points_local
    call global_interpolation_init_point_arrays(this)
    if (allocated(xyz)) then
       deallocate(xyz)
    end if
    allocate(xyz(3, n_points))
 
    call copy(xyz, this%xyz_local, 3*n_points)
    call copy(this%rst, this%rst_local, 3*n_points)
    call copy(this%xyz, this%xyz_local, 3*n_points)
    this%pe_owner = this%pe_rank
    this%el_owner0 = this%el_owner0_local
    if (neko_bcknd_device .eq. 1) then
       call device_memcpy(this%el_owner0, this%el_owner0_d, &
            this%n_points, host_to_device, sync = .true.)
    end if
    this%all_points_local = .true.
 
  end subroutine global_interpolation_find_and_redist
 
  subroutine global_interpolation_init_redist_comm(this, redist_comm)
    class(global_interpolation_t), intent(inout) :: this
    type(glb_intrp_comm_t), intent(inout) :: redist_comm
    type(stack_i4_t) :: send_pe
    type(stack_i4_t) :: recv_pe
    integer, pointer :: point_ids(:) => null()
    integer :: rank
    integer :: i
 
    call send_pe%init()
    call recv_pe%init()
    call redist_comm%init_dofs(this%pe_size)
 
    do rank = 0, this%pe_size - 1
       if (this%n_points_pe(rank) .gt. 0) then
          call send_pe%push(rank)
          point_ids => this%points_at_pe(rank)%array()
          do i = 1, this%n_points_pe(rank)
             call redist_comm%send_dof(rank)%push(point_ids(i))
          end do
       end if
       if (this%n_points_pe_local(rank) .gt. 0) then
          call recv_pe%push(rank)
          do i = 1, this%n_points_pe_local(rank)
             call redist_comm%recv_dof(rank)%push( &
                  this%n_points_offset_pe_local(rank) + i)
          end do
       end if
    end do
 
    call redist_comm%init(send_pe, recv_pe, this%comm)
    nullify(point_ids)
    call send_pe%free()
    call recv_pe%free()
  end subroutine global_interpolation_init_redist_comm
 
  subroutine global_interpolation_evaluate_masked(this, interp_values, &
       field, mask, on_host)
    class(global_interpolation_t), target, intent(inout) :: this
    real(kind=rp), intent(inout), target :: interp_values(this%n_points)
    real(kind=rp), intent(inout), target :: field(this%n_dof)
    type(mask_t), intent(in) :: mask
    logical, intent(in) :: on_host
 
    call vector_masked_gather_copy(this%masked_field, field, mask, this%n_dof)
    call this%evaluate(interp_values, this%masked_field%x, on_host)
 
  end subroutine global_interpolation_evaluate_masked
 
  subroutine global_interpolation_evaluate(this, interp_values, field, on_host)
    class(global_interpolation_t), target, intent(inout) :: this
    real(kind=rp), intent(inout), target :: interp_values(this%n_points)
    real(kind=rp), intent(inout), target :: field(this%nelv*this%Xh%lxyz)
    logical, intent(in) :: on_host
    type(c_ptr) :: interp_d
 
    if (.not. this%all_points_local) then
       call this%local_interp%evaluate(this%temp_local%x, &
            this%el_owner0_local, field, this%nelv, on_host)
       if (neko_bcknd_device .eq. 1 .and. .not. on_host) then
          call device_memcpy(this%temp_local%x, this%temp_local%x_d, &
               this%n_points_local, device_to_host, .true.)
       end if
       interp_values = 0.0_rp
       call this%glb_intrp_comm%sendrecv(this%temp_local%x, interp_values, &
            this%n_points_local, this%n_points)
       if (neko_bcknd_device .eq. 1 .and. .not. on_host .and. &
            this%n_points .gt. 0) then
          interp_d = device_get_ptr(interp_values)
          call device_memcpy(interp_values, interp_d, &
               this%n_points, host_to_device, .false.)
       end if
    else
       call this%local_interp%evaluate(interp_values, this%el_owner0_local, &
            field, this%nelv, on_host)
    end if
 
  end subroutine global_interpolation_evaluate
 
 
  function rst_cmp(rst1, rst2, res1, res2, tol) result(rst2_better)
    real(kind=rp) :: rst1(3), res1(3)
    real(kind=rp) :: rst2(3), res2(3)
    real(kind=dp) :: tol
    logical :: rst2_better
    !If rst1 is invalid and rst2 is valid, take rst2
    ! If both invalidl, take smallest residual
    if (abs(rst1(1)) .gt. 1.0_xp+tol .or. &
         abs(rst1(2)) .gt. 1.0_xp+tol .or. &
         abs(rst1(3)) .gt. 1.0_xp+tol) then
       if (abs(rst2(1)) .le. 1.0_xp+tol .and. &
            abs(rst2(2)) .le. 1.0_xp+tol .and. &
            abs(rst2(3)) .le. 1.0_xp+tol) then
          rst2_better = .true.
       else
          rst2_better = (norm2(real(res2, xp)) .lt. norm2(real(res1, xp)))
       end if
    else
       rst2_better = (norm2(real(res2, xp)) .lt. norm2(real(res1, xp)) .and. &
            abs(rst2(1)) .le. 1.0_xp+tol .and. &
            abs(rst2(2)) .le. 1.0_xp+tol .and. &
            abs(rst2(3)) .le. 1.0_xp+tol)
    end if
  end function rst_cmp
 
end module global_interpolation