Subroutine to compute high precision quadrature nodes.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
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 |
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