scale_polynomial.F90

Go to the documentation of this file.

!> Module common / Polynomial
!!
!! @par Description
!!      A module to provide utilities for polynomials
!!
!! @par Reference
!!
!! @author Yuta Kawai, Team SCALE
!!
#include "scaleFElib.h"
module scale_polynomial
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_const, only: &
    pi => const_pi
  use scale_io
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
 
  public :: polynomial_genlegendrepoly
  public :: polynomial_genlegendrepoly_sub
  public :: polynomial_gendlegendrepoly
 
  public :: polynomial_gengausslobattopt
  public :: polynomial_gengausslobattoptintweight
 
  public :: polynomial_gengausslegendrept
  public :: polynomial_gengausslegendreptintweight
 
  public :: polynomial_genlagrangepoly
  public :: polynomial_gendlagrangepoly_lglpt
  public :: polynomial_gendlagrangepoly
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  private :: gen_jacobigaussquadraturepts
 
  !-----------------------------------------------------------------------------
 
contains
  !> A function to obtain the Lagrange basis functions related to the Gauss-Legendre-Lobatto (GLL) points 
  !!
  !! @param Nord Order of Lagrange polynomial
  !! @param x_lgl Positions of GLL points
  !! @param x Positions where the Lagrange basis functions are evaluated 
!OCL SERIAL
  function polynomial_genlagrangepoly(Nord, x_lgl, x) result(l)
    implicit none
    
    integer, intent(in) :: nord
    real(rp), intent(in) :: x_lgl(nord+1)
    real(rp), intent(in) :: x(:)
    real(rp) :: l(size(x), nord+1)
 
    integer :: n, i
    real(rp) :: p_lgl(nord+1,nord+1)    
    real(rp) :: p(size(x),nord+1)
    real(rp) :: pr(size(x),nord+1)
 
    ! real(RP) :: w(Nord+1)
    ! real(RP) :: int_w(Nord+1)
    !---------------------------------------------------------------------------
 
    p_lgl(:,:)  = polynomial_genlegendrepoly(nord, x_lgl)
    p(:,:)  = polynomial_genlegendrepoly(nord, x)
    pr(:,:) = polynomial_gendlegendrepoly(nord, x, p)
 
    do n=1, nord+1
      do i=1, size(x)
        if ( abs(x(i)-x_lgl(n)) < 1e-15_rp ) then
          l(i,n) = 1.0_rp
        else
          l(i,n) = &
            (x(i) - 1.0_rp)*(x(i) + 1.0_rp)*pr(i,nord+1)              &
            / (dble(nord*(nord+1))*p_lgl(n,nord+1)*(x(i) - x_lgl(n)))
        end if
      end do
    end do
 
    !- Calculate interpolation coefficient based on barycentric form
    ! Eq. (3.46) in Wang, Huybrechs & Vandewalle (2012): Explicit barycentric weights for polynomial interpolation in the roots or extrema of classical orthogonal polynomials 
    ! int_w(:) = Polynomial_GenGaussLobattoPtIntWeight( Nord )
    ! do n=1, Nord+1
    !   w(n) = (-1)**mod(n-1,2) * sqrt(int_w(n))
    ! end do
    ! do n=1, Nord+1
    !   do i=1, size(x)
    !     if ( abs(x(i)-x_lgl(n)) < 1E-16_RP ) then
    !       l(i,n) = 1.0_RP
    !     else
    !       l(i,n) = &
    !           ( w(n) / ( x(i) - x_lgl(n) ) ) &
    !         / sum( w(:) / ( x(i) - x_lgl(:) ) )
    !     end if
    !   end do      
    ! end do
 
    return
  end function polynomial_genlagrangepoly
 
  !> Differential values of Lagrange basis functions at the GLL points
  !!
  !! @param Nord Order of Lagrange polynomial
  !! @param x_lgl Positions of GLL points
!OCL SERIAL
  function polynomial_gendlagrangepoly_lglpt(Nord, x_lgl) result(lr)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x_lgl(nord+1)
    real(rp) :: lr(nord+1, nord+1) !> lr(n,k) = derivative of k-th Lagrange basis at x_lgl(n)
 
    integer :: n, k
    real(rp) :: p(nord+1,nord+1)
    real(rp) :: lr_nn
    !---------------------------------------------------------------------------
 
    p(:,:)  = polynomial_genlegendrepoly(nord, x_lgl)
 
    do n=1, nord+1
      lr_nn = 0.0_rp
      do k=1, nord+1
        if (k==1 .and. n==1) then
          lr(k,n) = - 0.25_rp*dble(nord*(nord+1))
        else if (k==nord+1 .and. n==nord+1) then
          lr(k,n) = + 0.25_rp*dble(nord*(nord+1))
        else if (k==n) then
          lr(k,n) = 0.0_rp
        else
          lr(k,n) = p(n,nord+1)/(p(k,nord+1)*(x_lgl(n) - x_lgl(k)))
        end if
 
        if ( k /= n ) then
          lr_nn = lr_nn + lr(k,n)
        end if        
      end do
      lr(n,n) = - lr_nn
    end do
    
    return
  end function polynomial_gendlagrangepoly_lglpt
 
  !> Differential values of Lagrange basis functions defined at GLL nodes, evaluated at arbitrary points.
  !!
  !! @param Nord    Polynomial order
  !! @param x_lgl   GLL interpolation nodes
  !! @param Neval   Number of evaluation points
  !! @param x_eval  Evaluation points
!OCL SERIAL
  function polynomial_gendlagrangepoly( Nord, x_lgl, x_eval ) result(lr)
    implicit none
    integer, intent(in) :: nord
    real(rp), intent(in) :: x_lgl(nord+1)
    real(rp), intent(in) :: x_eval(:)
    real(rp) :: lr(size(x_eval),nord+1) !> lr(m,k) = derivative of k-th Lagrange basis at x_eval(m)
 
    real(rp) :: p(nord+1,nord+1)
    real(rp) :: bw(nord+1)       ! Barycentric weights. Common normalization factor is unnecessary.
 
    real(rp) :: a(nord+1)
    real(rp) :: s1, s2
    real(rp) :: dx
    real(rp) :: tol
 
    integer :: k, j, m
    integer :: knode
    !-------------------------------------------------------------
 
    p(:,:) = polynomial_genlegendrepoly(nord, x_lgl)
    bw(:) = 1.0_rp / p(:,nord+1)
 
    tol = 100.0_rp * epsilon(1.0_rp)
 
    do m = 1, size(x_eval)
 
      !- Check whether evaluation point coincides with a GLL node
      knode = 0
      do k=1, nord+1
        if ( abs(x_eval(m)-x_lgl(k)) <= tol ) then
          knode = k
          exit
        end if
      end do
 
      if ( knode /= 0 ) then
        !* Evaluation exactly at a GLL node. Use the existing nodal formula.
        do k = 1, nord+1
          if ( k /= knode ) then
            lr(m,k) = p(knode,nord+1) / ( p(k,nord+1) * ( x_lgl(knode)-x_lgl(k) ) )
          else
            lr(m,k) = 0.0_rp
          end if
        end do
 
        ! Sum of derivatives of all basis functions is zero: sum_k l'_k(x) = 0
        lr(m,knode) = -sum(lr(m,:))
 
      else
        !* General evaluation point
        !  l'_k(x) = l_k(x) * [ ( sum_j w_j/(x - x_j,lgl)**2 ) / [sum_j a_j(x)] - 1 / ( x - x_k,lgl ) ]
        ! where
        !  l_k(x) = a_k(x) / sum_j a_j(x),
        !  a_k(x) = w_k / ( x - x_k,lgl )
 
        s1 = 0.0_rp; s2 = 0.0_rp
        do j = 1, nord+1
          dx = x_eval(m) - x_lgl(j)
          
          a(j) = bw(j) / dx
          s1 = s1 + a(j)
          s2 = s2 + bw(j) / ( dx * dx )
        end do
 
        !　l'_k(x)
        do k = 1, nord+1
          dx = x_eval(m) - x_lgl(k)
          lr(m,k) = ( a(k) / s1 ) * ( s2 / s1 - 1.0_rp / dx )
        end do
 
      end if
 
    end do
    return
  end function polynomial_gendlagrangepoly
 
  !> A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points. 
  !!
  !! @param Nord Order of Lagrange polynomial
  !! @param x Positions where the Legendre polynomials are evaluated
  !! @param P Values of the Legendre polynomials at x
!OCL SERIAL
  subroutine polynomial_genlegendrepoly_sub(Nord, x, P)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp), intent(out) :: p(size(x), nord+1)
 
    integer :: n
    !---------------------------------------------------------------------------
 
    if (nord < 0) then
       log_error("Polynomial_GenLegendrePoly",*)   "Nord must be larger than 0."
    end if
 
    p(:,1) = 1.0_rp
    if (nord==0) return
 
    p(:,2) = x(:)
    do n=2, nord
      p(:,n+1) = ( dble(2*n-1)*x(:)*p(:,n) - dble(n-1)*p(:,n-1) )/dble(n)
    end do
    
    return    
  end subroutine polynomial_genlegendrepoly_sub
 
  !> A function to obtain the values of Legendre polynomials which are evaluated at arbitrary points.   
  !!
  !! @param Nord Order of Lagrange polynomial
  !! @param x Positions where the Legendre polynomials are evaluated
  !! @param P Values of the Legendre polynomials at x
!OCL SERIAL
  function polynomial_genlegendrepoly(Nord, x) result(P)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp) :: p(size(x), nord+1)
    !---------------------------------------------------------------------------
 
    call polynomial_genlegendrepoly_sub( nord, x(:), & ! (in)
       p(:,:)                                         ) ! (out)
    return    
  end function polynomial_genlegendrepoly
 
  !> A function to obtain differential values of Legendre polynomials which are evaluated at arbitrary points. 
  !! 
  !! @param Nord Order of Lagrange polynomial
  !! @param x Positions where the Legendre polynomials are evaluated
  !! @param P Values of the Legendre polynomials at x
!OCL SERIAL
  function polynomial_gendlegendrepoly(Nord, x, P) result(GradP)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp), intent(in) :: p(:,:)
    real(rp) :: gradp(size(x), nord+1)
 
    integer :: n
    !---------------------------------------------------------------------------
 
    if (nord < 0) then
      log_error("Polynomial_GenDLegendrePoly",*)   "Nord must be larger than 0."
    end if
 
    gradp(:,1) = 0.0_rp
    if (nord == 0) return
 
    gradp(:,2) = 1.0_rp
    do n=2, nord
      gradp(:,n+1) = 2.0_rp*x(:)*gradp(:,n) - gradp(:,n-1) + p(:,n)
    end do
 
    return    
  end function polynomial_gendlegendrepoly
 
  !> A function to calculate the Legendre-Gauss-Lobatto (LGL) points.
  !!
  !! @param Nord Order of Lagrange polynomial
  !! @param pts Position of the LGL points
!OCL SERIAL
  function polynomial_gengausslobattopt(Nord) result(pts)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: pts(nord+1)
 
    integer :: n1
    !---------------------------------------------------------------------------
 
    pts(1) = -1.0_rp; pts(nord+1) = 1.0_rp
    if (nord==1) return
 
    call gen_jacobigaussquadraturepts( 1, 1, nord-2, pts(2:nord) )
    return   
  end function polynomial_gengausslobattopt
 
  !> A function to calculate the Gauss-Lobbato weights. 
  !!  
  !! @param Nord Order of Lagrange polynomial
  !! @param int_weight_lgl Gauss-Lobbato weights
!OCL SERIAL
  function polynomial_gengausslobattoptintweight(Nord) result(int_weight_lgl)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: int_weight_lgl(nord+1)
 
    real(rp) :: lglpts1d(nord+1)
    real(rp) :: p1d_ori(nord+1, nord+1)
    !---------------------------------------------------------------------------
 
    lglpts1d(:)      = polynomial_gengausslobattopt( nord )
    p1d_ori(:,:)     = polynomial_genlegendrepoly( nord, lglpts1d )
 
    int_weight_lgl(:) = 2.0_rp/(dble(nord*(nord+1))*p1d_ori(:,nord+1)**2)
 
    return
  end function polynomial_gengausslobattoptintweight
 
  !> A function to calculate the Gauss-Legendre (GL) points.
  !!
  !! @param Nord Order of the Legendre polynomial
  !! @param pts Position of the GL points
!OCL SERIAL
  function polynomial_gengausslegendrept(Nord) result(pts)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: pts(nord)
    !---------------------------------------------------------------------------
 
    call gen_jacobigaussquadraturepts( 0, 0, nord-1, pts(:) )
    return   
  end function polynomial_gengausslegendrept
 
  !> A function to calculate the Gauss-Legendre (GL) weights. 
  !! 
  !! @param Nord Order of the Legendre polynomial
  !! @param int_weight_gl Gauss-Legendre weights
!OCL SERIAL
  function polynomial_gengausslegendreptintweight(Nord) result(int_weight_gl)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: int_weight_gl(nord)
 
    real(rp) :: glpts1d(nord)
    real(rp) :: p1d_ori(nord, nord+1)
    real(rp) :: dp1d_ori(nord, nord+1)
    !---------------------------------------------------------------------------
 
    glpts1d(:)    = polynomial_gengausslegendrept( nord )
    p1d_ori(:,:)    = polynomial_genlegendrepoly( nord, glpts1d )   
    dp1d_ori(:,:) = polynomial_gendlegendrepoly( nord, glpts1d, p1d_ori )
 
    int_weight_gl(:) = 2.0_rp / ( (1.0_rp - glpts1d(:)**2) * dp1d_ori(:,nord+1)**2 )
 
    return
  end function polynomial_gengausslegendreptintweight
 
  !- private -------------------------------
 
  !> Calculate the N'th-order Gauss quadrature points and weights associated the Jacobi polynomial of type (alpha,beta).
!OCL SERIAL
  subroutine gen_jacobigaussquadraturepts( alpha, beta, N, &
      x )
 
    implicit none
    integer, intent(in) :: alpha
    integer, intent(in) :: beta
    integer, intent(in) :: n
    real(rp), intent(out) :: x(n+1)
 
    integer :: i
    real(dp) :: d(n+1), e(n)
    real(dp) :: work(2*(n+1)-2), z(n+1,n+1)
    real(dp) :: h1(n+1)
    integer :: info
    !--------------------------------------------------------------
 
    if (n==0) then
      x(1) = - dble(alpha - beta) / dble(alpha + beta + 2)
      return
    end if
 
    do i=0, n
      h1(i+1) = dble( 2*i + alpha + beta )
    end do
 
    do i=1, n+1
      d(i) = - dble(alpha**2 - beta**2) / (h1(i) * (h1(i) + 2d0))
    end do
    do i=1, n
      e(i) = 2d0 / (h1(i) + 2d0) &
           * sqrt( dble(i * (i + alpha + beta) * (i + alpha) * (i + beta)) &
                   / ((h1(i) + 1d0) * (h1(i) + 3d0))                       )
    end do
    if ( dble(alpha + beta) < 1d-16) d(1) = 0.0_rp
    
    call dstev( 'Vectors', n+1, d, e, z, n+1, work, info )
    x(:) = d(:) 
 
    return
  end subroutine gen_jacobigaussquadraturepts
 
end module scale_polynomial