pacman_fortran.f90

Go to the documentation of this file.

!
! This file is subject to the terms and conditions defined in
! file 'LICENSE', which is part of this source code package.
!
 
!> @brief Fortran module providing a high-level interface to the PACMAN
!>        interpolation library via ISO C bindings.
!>
!> Usage example (RBF-PUM):
!> @code{.f90}
!>   use pacman_mod
!>   real(8), allocatable :: srcPts(:,:), srcVals(:), tgtPts(:,:), tgtVals(:)
!>   integer :: ierr
!>
!>   call pacman_kokkos_initialize()
!>   allocate(tgtVals(nTarget))
!>
!>   call pacman_rbf_interpolate(3, PACMAN_OPENMP, PACMAN_WENDLANDC2, &
!>                               srcPts, srcVals, tgtPts, tgtVals, ierr)
!>
!>   call pacman_kokkos_finalize()
!> @endcode
!>
!> Point arrays must be laid out as @c pts(spaceDimension, nPoints), which
!> in Fortran column-major storage is identical to C row-major
!> @c pts[nPoints][spaceDimension].
module pacman_mod
  use iso_c_binding
  implicit none
  private
 
  ! -----------------------------------------------------------------
  ! Public procedures
  ! -----------------------------------------------------------------
  public :: pacman_kokkos_initialize
  public :: pacman_kokkos_finalize
  public :: pacman_best_execspace
  public :: pacman_rbf_interpolate
  public :: pacman_fe_interpolate
  public :: pacman_mls_interpolate
  public :: pacman_vtk_to_pacman_cell_type
 
  ! -----------------------------------------------------------------
  ! Execution space constants (match PACMAN::ExecSpaces)
  ! -----------------------------------------------------------------
  integer(c_int), parameter, public :: pacman_serial  = 0
  integer(c_int), parameter, public :: pacman_openmp  = 1
  integer(c_int), parameter, public :: pacman_threads = 2
  integer(c_int), parameter, public :: pacman_cuda    = 3
  integer(c_int), parameter, public :: pacman_hip     = 4
  integer(c_int), parameter, public :: pacman_sycl    = 5
 
  ! -----------------------------------------------------------------
  ! RBF function constants (match PACMAN::RbfFunctions)
  !   0x10 = 16, 0x11 = 17, 0x12 = 18, 0x13 = 19, 0x14 = 20
  ! -----------------------------------------------------------------
  integer(c_int), parameter, public :: pacman_wendlandc0 = 16
  integer(c_int), parameter, public :: pacman_wendlandc2 = 17
  integer(c_int), parameter, public :: pacman_wendlandc4 = 18
  integer(c_int), parameter, public :: pacman_wendlandc6 = 19
  integer(c_int), parameter, public :: pacman_wendlandc8 = 20
 
  ! -----------------------------------------------------------------
  ! FE method constants (match PACMAN::TransferMethods)
  !   0xF0 = 240, 0xF1 = 241, ..., 0xF4 = 244
  ! -----------------------------------------------------------------
  integer(c_int), parameter, public :: pacman_fe_nearest_nearest = 240
  integer(c_int), parameter, public :: pacman_fe_interp_clamp    = 241
  integer(c_int), parameter, public :: pacman_fe_interp_nearest  = 242
  integer(c_int), parameter, public :: pacman_fe_interp_zerofill = 243
  integer(c_int), parameter, public :: pacman_fe_interp_extrap   = 244
 
  ! -----------------------------------------------------------------
  ! MLS constant (match PACMAN::TransferMethods::MLS = 0xF6 = 246)
  ! -----------------------------------------------------------------
  integer(c_int), parameter, public :: pacman_mls = 246
 
  ! -----------------------------------------------------------------
  ! ISO C interfaces to the plain-C shim (fortran_interface.cpp)
  ! -----------------------------------------------------------------
  interface
    subroutine pacman_kokkos_initialize_c() &
        bind(c, name='pacman_kokkos_initialize_c')
    end subroutine
 
    subroutine pacman_kokkos_finalize_c() &
        bind(c, name='pacman_kokkos_finalize_c')
    end subroutine
 
    function pacman_best_execspace_c() result(es) &
        bind(c, name='pacman_best_execspace_c')
      import :: c_int
      integer(c_int) :: es
    end function
 
    function pacman_rbf_interpolate_c(                           &
        spaceDimension, execSpace, rbfFunction,                  &
        sourcePoints, nSourcePoints, sourceValues,               &
        targetPoints, nTargetPoints, targetValues) result(ierr)  &
        bind(c, name='pacman_rbf_interpolate_c')
      import :: c_int, c_ptr
      integer(c_int), value :: spacedimension
      integer(c_int), value :: execspace
      integer(c_int), value :: rbffunction
      type(c_ptr),    value :: sourcepoints
      integer(c_int), value :: nsourcepoints
      type(c_ptr),    value :: sourcevalues
      type(c_ptr),    value :: targetpoints
      integer(c_int), value :: ntargetpoints
      type(c_ptr),    value :: targetvalues
      integer(c_int)        :: ierr
    end function
 
    function pacman_fe_interpolate_c(                                       &
        spaceDimension, execSpace, method,                                  &
        sourcePoints, nSourcePoints, sourceValues,                          &
        connVal, connValSize, connOff, connOffSize, cellTypes,              &
        targetPoints, nTargetPoints, targetValues, targetStatus) result(ierr) &
        bind(c, name='pacman_fe_interpolate_c')
      import :: c_int, c_ptr
      integer(c_int), value :: spacedimension
      integer(c_int), value :: execspace
      integer(c_int), value :: method
      type(c_ptr),    value :: sourcepoints
      integer(c_int), value :: nsourcepoints
      type(c_ptr),    value :: sourcevalues
      type(c_ptr),    value :: connval
      integer(c_int), value :: connvalsize
      type(c_ptr),    value :: connoff
      integer(c_int), value :: connoffsize
      type(c_ptr),    value :: celltypes
      type(c_ptr),    value :: targetpoints
      integer(c_int), value :: ntargetpoints
      type(c_ptr),    value :: targetvalues
      type(c_ptr),    value :: targetstatus
      integer(c_int)        :: ierr
    end function
 
    function pacman_mls_interpolate_c(                                  &
        spaceDimension, execSpace,                                       &
        sourcePoints, nSourcePoints, sourceValues,                       &
        targetPoints, nTargetPoints, targetValues) result(ierr)          &
        bind(c, name='pacman_mls_interpolate_c')
      import :: c_int, c_ptr
      integer(c_int), value :: spacedimension
      integer(c_int), value :: execspace
      type(c_ptr),    value :: sourcepoints
      integer(c_int), value :: nsourcepoints
      type(c_ptr),    value :: sourcevalues
      type(c_ptr),    value :: targetpoints
      integer(c_int), value :: ntargetpoints
      type(c_ptr),    value :: targetvalues
      integer(c_int)        :: ierr
    end function
 
    function pacman_vtk_to_pacman_cell_type_c(vtkTypes, pacmanTypes, n) &
        result(ierr) bind(c, name='pacman_vtk_to_pacman_cell_type_c')
      import :: c_int, c_ptr
      type(c_ptr),    value :: vtktypes
      type(c_ptr),    value :: pacmantypes
      integer(c_int), value :: n
      integer(c_int)        :: ierr
    end function
  end interface
 
contains
 
  ! -----------------------------------------------------------------
  ! Kokkos lifecycle wrappers
  ! -----------------------------------------------------------------
 
  !> @brief Initialize Kokkos with default settings.
  !> Must be called before any interpolation subroutine.
  subroutine pacman_kokkos_initialize()
    call pacman_kokkos_initialize_c()
  end subroutine
 
  !> @brief Finalize Kokkos.
  !> Must be called after all interpolation calls are complete.
  subroutine pacman_kokkos_finalize()
    call pacman_kokkos_finalize_c()
  end subroutine
 
  !> @brief Return the best execution space available in the current build.
  !> @return Integer constant from PACMAN_SERIAL / PACMAN_OPENMP / etc.
  function pacman_best_execspace() result(es)
    integer(c_int) :: es
    es = pacman_best_execspace_c()
  end function
 
  ! -----------------------------------------------------------------
  ! RBF-PUM interpolation
  ! -----------------------------------------------------------------
 
  !> @brief Interpolate scalar values from source to target points using
  !>        RBF-PUM (no mesh connectivity required).
  !>
  !> @param spaceDimension  Geometric dimension (1, 2, or 3).
  !> @param execSpace       Execution-space constant (e.g. PACMAN_OPENMP).
  !> @param rbfFunction     RBF basis-function constant (e.g. PACMAN_WENDLANDC2).
  !> @param sourcePoints    Source coordinates, shape (spaceDimension, nSource).
  !> @param sourceValues    Source scalar values, shape (nSource).
  !> @param targetPoints    Target coordinates, shape (spaceDimension, nTarget).
  !> @param targetValues    Interpolated output, shape (nTarget).  Must be
  !>                        allocated by the caller before the call.
  !> @param ierr            Optional error code: 0 = success.
  subroutine pacman_rbf_interpolate(spaceDimension, execSpace, rbfFunction, &
      sourcePoints, sourceValues, targetPoints, targetValues, ierr)
    integer,        intent(in)            :: spacedimension
    integer(c_int), intent(in)            :: execspace
    integer(c_int), intent(in)            :: rbffunction
    real(c_double), intent(in),  contiguous, target :: sourcepoints(:,:)
    real(c_double), intent(in),  contiguous, target :: sourcevalues(:)
    real(c_double), intent(in),  contiguous, target :: targetpoints(:,:)
    real(c_double), intent(out), contiguous, target :: targetvalues(:)
    integer(c_int), intent(out), optional :: ierr
 
    integer(c_int) :: ret, nsrc, ntgt
 
    nsrc = int(size(sourcevalues), c_int)
    ntgt = int(size(targetvalues), c_int)
    ret  = pacman_rbf_interpolate_c(                                        &
               int(spacedimension, c_int), execspace, rbffunction,          &
               c_loc(sourcepoints), nsrc, c_loc(sourcevalues),              &
               c_loc(targetpoints), ntgt, c_loc(targetvalues))
    if (present(ierr)) ierr = ret
  end subroutine
 
  ! -----------------------------------------------------------------
  ! FE interpolation
  ! -----------------------------------------------------------------
 
  !> @brief Interpolate scalar values from a source mesh to target points
  !>        using finite-element transfer methods.
  !>
  !> @param spaceDimension  Geometric dimension (1, 2, or 3).
  !> @param execSpace       Execution-space constant.
  !> @param method          FE method constant (e.g. PACMAN_FE_INTERP_CLAMP).
  !> @param sourcePoints    Source coordinates, shape (spaceDimension, nSource).
  !> @param sourceValues    Source scalar values, shape (nSource).
  !> @param connVal         CSR connectivity values, shape (connValSize).
  !> @param connOff         CSR connectivity offsets, shape (nElems + 1).
  !> @param cellTypes       PACMAN cell-type codes, shape (nElems).
  !> @param targetPoints    Target coordinates, shape (spaceDimension, nTarget).
  !> @param targetValues    Interpolated output, shape (nTarget).
  !> @param targetStatus    Transfer-status codes per target point, shape (nTarget).
  !> @param ierr            Optional error code: 0 = success.
  subroutine pacman_fe_interpolate(spaceDimension, execSpace, method,          &
      sourcePoints, sourceValues, connVal, connOff, cellTypes,                 &
      targetPoints, targetValues, targetStatus, ierr)
    integer,        intent(in)            :: spacedimension
    integer(c_int), intent(in)            :: execspace
    integer(c_int), intent(in)            :: method
    real(c_double), intent(in),  contiguous, target :: sourcepoints(:,:)
    real(c_double), intent(in),  contiguous, target :: sourcevalues(:)
    integer(c_int), intent(in),  contiguous, target :: connval(:)
    integer(c_int), intent(in),  contiguous, target :: connoff(:)
    integer(c_int), intent(in),  contiguous, target :: celltypes(:)
    real(c_double), intent(in),  contiguous, target :: targetpoints(:,:)
    real(c_double), intent(out), contiguous, target :: targetvalues(:)
    integer(c_int), intent(out), contiguous, target :: targetstatus(:)
    integer(c_int), intent(out), optional :: ierr
 
    integer(c_int) :: ret, nsrc, ntgt, nconnval, nconnoff
    integer :: i
 
    nsrc     = int(size(sourcevalues), c_int)
    ntgt     = int(size(targetvalues), c_int)
    nconnval = int(size(connval),      c_int)
    nconnoff = int(size(connoff),      c_int)
    write(*,*) 'DEBUG: nsrc=', nsrc, ' ntgt=', ntgt, ' nconnval=', nconnval, ' nconnoff=', nconnoff, &
     ' size(connVal)=', size(connval), ' size(connOff)=', size(connoff)
    do i = 1, size(connval)
      if (connval(i) < 0) then
        ! LCOV_EXCL_START
        write(*,'(A,I0,A,I0)') 'ERROR: connVal(', i, ') = ', connval(i)
        if (present(ierr)) ierr = -1
        return
        ! LCOV_EXCL_STOP
      end if
    end do
    ret = pacman_fe_interpolate_c(                                             &
              int(spacedimension, c_int), execspace, method,                   &
              c_loc(sourcepoints), nsrc,  c_loc(sourcevalues),                 &
              c_loc(connval), nconnval, c_loc(connoff), nconnoff, c_loc(celltypes), &
              c_loc(targetpoints), ntgt, c_loc(targetvalues), c_loc(targetstatus))
    if (present(ierr)) ierr = ret
  end subroutine
 
  ! -----------------------------------------------------------------
  ! MLS interpolation
  ! -----------------------------------------------------------------
 
  !> @brief Interpolate scalar values from source to target points using
  !>        Moving Least Squares (no mesh connectivity required).
  !>
  !> @param spaceDimension  Geometric dimension (1, 2, or 3).
  !> @param execSpace       Execution-space constant (e.g. PACMAN_OPENMP).
  !> @param sourcePoints    Source coordinates, shape (spaceDimension, nSource).
  !> @param sourceValues    Source scalar values, shape (nSource).
  !> @param targetPoints    Target coordinates, shape (spaceDimension, nTarget).
  !> @param targetValues    Interpolated output, shape (nTarget).  Must be
  !>                        allocated by the caller before the call.
  !> @param ierr            Optional error code: 0 = success.
  subroutine pacman_mls_interpolate(spaceDimension, execSpace, &
      sourcePoints, sourceValues, targetPoints, targetValues, ierr)
    integer,        intent(in)            :: spacedimension
    integer(c_int), intent(in)            :: execspace
    real(c_double), intent(in),  contiguous, target :: sourcepoints(:,:)
    real(c_double), intent(in),  contiguous, target :: sourcevalues(:)
    real(c_double), intent(in),  contiguous, target :: targetpoints(:,:)
    real(c_double), intent(out), contiguous, target :: targetvalues(:)
    integer(c_int), intent(out), optional :: ierr
 
    integer(c_int) :: ret, nsrc, ntgt
 
    nsrc = int(size(sourcevalues), c_int)
    ntgt = int(size(targetvalues), c_int)
    ret  = pacman_mls_interpolate_c(                                        &
               int(spacedimension, c_int), execspace,                       &
               c_loc(sourcepoints), nsrc, c_loc(sourcevalues),              &
               c_loc(targetpoints), ntgt, c_loc(targetvalues))
    if (present(ierr)) ierr = ret
  end subroutine
 
  ! -----------------------------------------------------------------
  ! VTK cell-type helpers
  ! -----------------------------------------------------------------
 
  !> @brief Convert VTK cell-type IDs to PACMAN CellType codes.
  !>
  !> @param vtkTypes    Input VTK cell-type IDs, shape (n).
  !> @param pacmanTypes Output PACMAN cell-type codes, shape (n).
  !> @param ierr        Optional error code: 0 = success.
  subroutine pacman_vtk_to_pacman_cell_type(vtkTypes, pacmanTypes, ierr)
    integer(c_int), intent(in),  contiguous, target :: vtktypes(:)
    integer(c_int), intent(out), contiguous, target :: pacmantypes(:)
    integer(c_int), intent(out), optional           :: ierr
    integer(c_int) :: ret, n
    n   = int(size(vtktypes), c_int)
    ret = pacman_vtk_to_pacman_cell_type_c(c_loc(vtktypes), c_loc(pacmantypes), n)
    if (present(ierr)) ierr = ret
  end subroutine
 
end module pacman_mod