0
$\begingroup$

I've successfully implemented a 1D DG code with non-normalized Legendre basis and I've now moved onto developing a 2D code using tensor products. For my 2D code I've chosen to have normalized Legendre polynomials,i.e.

$$ L_{0}(x) = \sqrt{\frac{1}{2}}\\ L_{1}(x) = \sqrt{\frac{3}{2}} x. $$

My question is: If we normalize the Legendre polynomials, do I also have to normalized my Gauss-Legendre weights as well? Or do we only normalized the polynomials in terms of the basis?

subroutine GaussQuad (xq,wq,n)
use parameters 
implicit none 
integer :: n
real(kind=8),dimension(n) :: xq,wq

integer :: i,iter
real(kind=8) :: xx
real(kind=8) :: legendre,dlegendre

do i=1,n
 xx = cos(dpi*(i - 0.25d0)/(n + 0.5d0))

 do iter=1,500
    xx = xx - legendre(xx,n)/dlegendre(xx,n)
 end do

 xq(i) = xx
 wq(i) = (2.0*dble(n) + 1.0)*2.0d0/((1.0d0-xx**2.0)*dlegendre(xx,n)**2.0)
end do
end subroutine GaussQuad

In order for my 2D code to work I need to use the normalized Legendre polynomials in the Gauss quadrature routine along with the $2n + 1$ normalization on the weights (see how its been added to wq(i)). However, I only got this to work due to an ad hoc guess. I would like to avoid this as I don't personally understand why this is required for my solver to work.

EDIT 4/1/2020:

  1. In my code, I do all my operations in the reference element. In fact, double checked and my Gauss-Legendre points values are within $[-1,+1]$ and my weights also sum to 2.
  2. To map from physical space to the reference space, I use the following $$ X(\eta) = \frac{x_{i} - x_{i-1}}{2} \eta + \frac{x_{i} + x_{i-1}}{2}\\ Y(\zeta) = \frac{y_{j} - y_{j-1}}{2} \zeta + \frac{y_{j} + y_{j-1}}{2} $$
  3. To integrate the volume fluxes I use a tensor product Gauss-Legendre scheme with $M^{2}$ points

    do ix=1,nx
     do iy=1,ny
      do i=1,mx
       do j=1,my
        do inode=1,mx
          do jnode=1,my
           call Flux(un(:,ix,iy,inode,jnode),FFlux,GFlux)
    
           flux_vol1(:,ix,iy,i,j) = flux_vol1(:,ix,iy,i,j) + &
           & 0.5*FFlux(:)* &
           & dlegendre(xg(inode),i-1)*& 
           & legendre(xg(jnode),j-1)*&
           & wg(inode)*&
           & wg(jnode)
    
           flux_vol2(:,ix,iy,i,j) = flux_vol2(:,ix,iy,i,j) + &
           & 0.5*GFlux(:)* &
           & dlegendre(xg(jnode),j-1)*&
           & legendre(xg(inode),i-1)*&
           & wg(inode)* &
           & wg(jnode)
         end do
        end do
       end do
      end do
     end do
    end do
    

Above is a snippet in how I perform the tensor product quadrature. I'm a little confused on why removing the normalizing of the Gauss quadrature is causing such severe results.
Solution:

The standard way of computing tensor products is $$ u_{h}(x,y) = \sum_{k = 0}^{N} \sum_{l = 0}^{M} u_{ij}(t) \phi_{k}(x) \phi_{l}(y) (1), $$ The tensor product formula adapted in the last two papers are $$ U_{h}(x,y) = \sum_{k=0}^{\frac{1}{2}(k + 1)(k + 2)} u_{ij}^{k}(t) \phi_{k}(x,y) (2) $$ Is there a reason to use formula (1) over formula (2)? Mathematically, they seem both correct.

$\endgroup$
  • 2
    $\begingroup$ Quadrature weights/nodes do not depend on the normalisation of the basis functions. They depend on your reference interval. Usually it is [-1,1] since Legendre are defined on this interval and the weights then sum to 2. But if you use the reference interval as [0,1], then you must map the nodes to lie in this interval and weights must be scaled to sum to unity. $\endgroup$ – cfdlab Apr 1 at 6:45
  • $\begingroup$ That makes sense. I included some more info to maybe help illustrate my problem. Its quite perplexing. BTW the PDE i'm solving is 2D linear advection with a Gaussian pulse (very smooth). $\endgroup$ – NumericalKid Apr 1 at 16:02
  • $\begingroup$ My answer was only about using Legendre as basis in which case normalisation does not matter. I see that you are trying to compute the weights yourself. In this case, the weights do depend on the normalisation. The formulae e.g., here en.wikipedia.org/wiki/Gaussian_quadrature#Gauss–Legendre_quadrature are without normalisation. The weights are also given in this site, you can directly use them as they are exact. $\endgroup$ – cfdlab Apr 2 at 3:02
  • $\begingroup$ Hi CFD Lab. I managed to fix my issue. I have included what I did along with a side question. $\endgroup$ – NumericalKid Apr 8 at 16:05
  • $\begingroup$ Both types can be used in DG since there are no continuity requirements. Tensor product polynomials are normally used on quadrilateral/hexahedral elements. The other type are used on any type of elements. It is better to ask this as a new question, maybe somebody else will give more insight. $\endgroup$ – cfdlab Apr 9 at 5:00
1
$\begingroup$

Per the regulations of Stack Exchange. I have included the "answer" to my question.

So I managed to fix my problem. First, off I was missing a factor of 0.5 when calculating the surface integral. Additionally, I am using the following normalization $\phi(x) = \sqrt{2m + 1} L(\xi)$.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.