sdc_qnodes Subroutine

public subroutine sdc_qnodes(qnodes, flags, qtype, nnodes)

Subroutine to compute high precision quadrature nodes.

Arguments

Type IntentOptional AttributesName
real(kind=pfqp), intent(out) :: qnodes(nnodes)

The computed quadrature nodes

integer, intent(out) :: flags(nnodes)
integer, intent(in), value:: qtype

Type of nodes (see pf_dtype)

integer, intent(in), value:: nnodes

Number of nodes


Calls

proc~~sdc_qnodes~~CallsGraph proc~sdc_qnodes sdc_qnodes proc~poly_legendre poly_legendre proc~sdc_qnodes->proc~poly_legendre proc~poly_diff poly_diff proc~sdc_qnodes->proc~poly_diff proc~poly_roots poly_roots proc~sdc_qnodes->proc~poly_roots proc~qsort qsort proc~poly_roots->proc~qsort proc~poly_eval_complex poly_eval_complex proc~poly_roots->proc~poly_eval_complex

Called by

proc~~sdc_qnodes~~CalledByGraph proc~sdc_qnodes sdc_qnodes proc~pf_quadrature pf_quadrature proc~pf_quadrature->proc~sdc_qnodes proc~pf_init_sdcmats pf_init_sdcmats proc~pf_init_sdcmats->proc~pf_quadrature proc~pf_level_setup pf_level_setup proc~pf_level_setup->proc~pf_init_sdcmats proc~pf_pfasst_setup pf_pfasst_setup proc~pf_pfasst_setup->proc~pf_level_setup

Contents

Source Code


Source Code

  subroutine sdc_qnodes(qnodes, flags, qtype, nnodes) 
    integer ,       intent(in), value  :: nnodes          !!  Number of nodes
    integer ,       intent(in), value  :: qtype           !!  Type of nodes (see pf_dtype)
    real(pfqp),  intent(out)        :: qnodes(nnodes)  !!  The computed quadrature nodes
    integer ,       intent(out)        :: flags(nnodes)   !!

    integer :: j, degree
    real(pfqp), allocatable :: roots(:)
    real(pfqp), allocatable :: coeffs(:), coeffs2(:)

    real(pfqp), parameter :: pi = 3.141592653589793115997963468544185161590576171875_pfdp

    flags = 0

    select case(qtype)

    case (SDC_GAUSS_LEGENDRE)

       degree = nnodes-2
       allocate(roots(degree))
       allocate(coeffs(degree+1))

       call poly_legendre(coeffs, degree)
       call poly_roots(roots, coeffs, degree)

       qnodes(1) = 0.0_pfqp
       qnodes(2:nnodes-1) = 0.5_pfqp * (1.0_pfqp + roots)
       qnodes(nnodes) = 1.0_pfqp

       deallocate(coeffs)
       deallocate(roots)

       do j = 2, nnodes-1
          flags(j) = ibset(flags(j), 0)
       end do

    case (SDC_GAUSS_LOBATTO)
       
       degree = nnodes - 1
       allocate(roots(degree-1))
       allocate(coeffs(degree+1))

       call poly_legendre(coeffs, degree)
       call poly_diff(coeffs, degree)
       call poly_roots(roots, coeffs(:degree), degree-1)

       qnodes(1)          = 0.0_pfdp
       qnodes(2:nnodes-1) = 0.5_pfdp * (1.0_pfdp + roots)
       qnodes(nnodes)     = 1.0_pfdp

       deallocate(coeffs)
       deallocate(roots)

       do j = 1, nnodes
          flags(j) = ibset(flags(j), 0)
       end do

    case (SDC_GAUSS_RADAU)

       degree = nnodes - 1
       allocate(roots(degree))
       allocate(coeffs(degree+1))
       allocate(coeffs2(degree))


       call poly_legendre(coeffs, degree)
       call poly_legendre(coeffs2, degree-1)
       coeffs(:degree) = coeffs(:degree) + coeffs2
       call poly_roots(roots, coeffs, degree)

       qnodes(1)      = 0.0_pfdp
       do j = 2, nnodes-1
          qnodes(j) = 0.5_pfdp * (1.0_pfdp - roots(nnodes+1-j))
       end do
       qnodes(nnodes) = 1.0_pfdp

       deallocate(coeffs2)
       deallocate(coeffs)
       deallocate(roots)

       do j = 2, nnodes
          flags(j) = ibset(flags(j), 0)
       end do

    case (SDC_CLENSHAW_CURTIS)

       do j = 0, nnodes-1
          qnodes(j+1) = 0.5_pfqp * (1.0_pfdp - cos(j * pi / (nnodes-1)))
       end do

       do j = 1, nnodes
          flags(j) = ibset(flags(j), 0)
       end do

    case (SDC_UNIFORM)

       do j = 0, nnodes-1
          qnodes(j+1) = j * (1.0_pfqp / (nnodes-1))
       end do

       do j = 1, nnodes
          flags(j) = ibset(flags(j), 0)
       end do

    case default
       print *,'qtype = ',qtype
       stop "ERROR: Invalid qtype in sdc_quadrature.f90."

    end select

  end subroutine sdc_qnodes