SELF_Lagrange.f90 Source File


Contents

Source Code


Source Code

! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
!
! Maintainers : support@fluidnumerics.com
! Official Repository : https://github.com/FluidNumerics/self/
!
! Copyright © 2024 Fluid Numerics LLC
!
! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
!
! 2. 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.
!
! 3. Neither the name of the copyright holder 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
! HOLDER 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 SELF_Lagrange

  use iso_fortran_env
  use iso_c_binding

  use SELF_Constants
  use SELF_Lagrange_t
  use SELF_GPU

  implicit none

  type,extends(Lagrange_t),public :: Lagrange
    character(3) :: backend = 'gpu'
    type(c_ptr) :: qWeights_gpu
    type(c_ptr) :: iMatrix_gpu
    type(c_ptr) :: dMatrix_gpu
    type(c_ptr) :: dgMatrix_gpu
    type(c_ptr) :: bMatrix_gpu

  contains

    procedure,public :: Init => Init_Lagrange
    procedure,public :: Free => Free_Lagrange

  endtype Lagrange

contains

  subroutine Init_Lagrange(this,N,controlNodeType,M,targetNodeType)
    !! Initialize an instance of the Lagrange class
    !! On output, all of the attributes for the Lagrange class are allocated and values are initialized according to the number of
    !! control points, number of target points, and the types for the control and target nodes.
    !! If a GPU is available, device pointers for the Lagrange attributes are allocated and initialized.
    implicit none
    class(Lagrange),intent(out) :: this
    !! Lagrange class instance
    integer,intent(in)          :: N
    !! The number of control points for interpolant
    integer,intent(in)          :: M
    !! The number of target points for the interpolant
    integer,intent(in)          :: controlNodeType
    !! The integer code specifying the type of control points. Parameters are defined in SELF_Constants.f90. One of GAUSS(=1),
    !! GAUSS_LOBATTO(=2), or UNIFORM(=3)
    integer,intent(in)          :: targetNodeType
    !! The integer code specifying the type of target points. Parameters are defined in SELF_Constants.f90. One of GAUSS(=1),
    !! GAUSS_LOBATTO(=2), or UNIFORM(=3)
    ! -------!
    ! Local
    real(prec) :: q(0:M)

    if(.not. GPUAvailable()) then
      print*,__FILE__,':',__LINE__,' : Error : Attempt to use GPU extension, but GPU is not available.'
      stop 1
    endif

    this%N = N
    this%M = M
    this%controlNodeType = controlNodeType
    this%targetNodeType = targetNodeType
    allocate(this%controlPoints(1:N+1), &
             this%targetPoints(1:M+1), &
             this%bWeights(1:N+1), &
             this%qWeights(1:N+1), &
             this%iMatrix(1:N+1,1:M+1), &
             this%dMatrix(1:N+1,1:N+1), &
             this%dgMatrix(1:N+1,1:N+1), &
             this%bMatrix(1:N+1,1:2))

    if(controlNodeType == GAUSS .or. controlNodeType == GAUSS_LOBATTO) then

      call LegendreQuadrature(N, &
                              this%controlPoints, &
                              this%qWeights, &
                              controlNodeType)

    elseif(controlNodeType == CHEBYSHEV_GAUSS .or. controlNodeType == CHEBYSHEV_GAUSS_LOBATTO) then

      call ChebyshevQuadrature(N, &
                               this%controlPoints, &
                               this%qWeights, &
                               controlNodeType)

    elseif(controlNodeType == UNIFORM) then

      this%controlPoints = UniformPoints(-1.0_prec,1.0_prec,0,N)
      this%qWeights = 2.0_prec/real(N,prec)

    endif

    ! Target Points
    if(targetNodeType == GAUSS .or. targetNodeType == GAUSS_LOBATTO) then

      call LegendreQuadrature(M, &
                              this%targetPoints, &
                              q, &
                              targetNodeType)

    elseif(targetNodeType == UNIFORM) then

      this%targetPoints = UniformPoints(-1.0_prec,1.0_prec,0,M)

    endif

    call this%CalculateBarycentricWeights()
    call this%CalculateInterpolationMatrix()
    call this%CalculateDerivativeMatrix()
    this%bMatrix(1:N+1,1) = this%CalculateLagrangePolynomials(-1.0_prec)
    this%bMatrix(1:N+1,2) = this%CalculateLagrangePolynomials(1.0_prec)

    print*,"Lagrange malloc"
    call gpuCheck(hipMalloc(this%iMatrix_gpu,sizeof(this%iMatrix)))
    call gpuCheck(hipMalloc(this%dMatrix_gpu,sizeof(this%dMatrix)))
    call gpuCheck(hipMalloc(this%dgMatrix_gpu,sizeof(this%dgMatrix)))
    call gpuCheck(hipMalloc(this%bMatrix_gpu,sizeof(this%bMatrix)))
    call gpuCheck(hipMalloc(this%qWeights_gpu,sizeof(this%qWeights)))

    print*,"Lagrange memcpy"

    print*,c_loc(this%iMatrix)
    call gpuCheck(hipMemcpy(this%iMatrix_gpu, &
                            c_loc(this%iMatrix), &
                            sizeof(this%iMatrix), &
                            hipMemcpyHostToDevice))
    print*,"Lagrange memcpy"

    call gpuCheck(hipMemcpy(this%dMatrix_gpu, &
                            c_loc(this%dMatrix), &
                            sizeof(this%dMatrix), &
                            hipMemcpyHostToDevice))
    print*,"Lagrange memcpy"

    call gpuCheck(hipMemcpy(this%dgMatrix_gpu, &
                            c_loc(this%dgMatrix), &
                            sizeof(this%dgMatrix), &
                            hipMemcpyHostToDevice))
    print*,"Lagrange memcpy"

    call gpuCheck(hipMemcpy(this%bMatrix_gpu, &
                            c_loc(this%bMatrix), &
                            sizeof(this%bMatrix), &
                            hipMemcpyHostToDevice))
    print*,"Lagrange memcpy"

    call gpuCheck(hipMemcpy(this%qWeights_gpu, &
                            c_loc(this%qWeights), &
                            sizeof(this%qWeights), &
                            hipMemcpyHostToDevice))

  endsubroutine Init_Lagrange

  subroutine Free_Lagrange(this)
    !! Frees all memory (host and device) associated with an instance of the Lagrange class
    implicit none
    class(Lagrange),intent(inout) :: this
    !! Lagrange class instance

    deallocate(this%controlPoints)
    deallocate(this%targetPoints)
    deallocate(this%bWeights)
    deallocate(this%qWeights)
    deallocate(this%iMatrix)
    deallocate(this%dMatrix)
    deallocate(this%dgMatrix)
    deallocate(this%bMatrix)

    call gpuCheck(hipFree(this%iMatrix_gpu))
    call gpuCheck(hipFree(this%dMatrix_gpu))
    call gpuCheck(hipFree(this%dgMatrix_gpu))
    call gpuCheck(hipFree(this%bMatrix_gpu))
    call gpuCheck(hipFree(this%qWeights_gpu))

  endsubroutine Free_Lagrange

endmodule SELF_Lagrange