On the Stability of the Flux Reconstruction Schemes on Quadrilateral Elements for the Linear Advection Equation

Abstract

The flux reconstruction (FR) approach to high-order methods has proved to be a promising alternative to traditional discontinuous Galerkin (DG) schemes since they facilitate the adoption of explicit time-stepping methods suitable for parallel architectures like GPUs. The FR approach provides a parameterized family of schemes through which various classical schemes like nodal-DG and spectral difference methods can be recovered. Further, the parameters can be varied to obtain schemes with a maximum stable time-step, or minimum dispersion or dissipation errors etc., providing us a single powerful framework unifying high-order discontinuous Finite Element Methods. There have been various studies on the accuracy and the stability of these schemes and in particular, a subset of the FR schemes known as ESFR or VCJH schemes have been shown to be stable in 1D and on simplex elements in 2D and 3D for the linear advection as well as the advection–diffusion equations. However, the stability of the FR schemes on tensor product quadrilateral elements has remained an open question. Although it is the most natural extension of the 1D FR approach, it has posed a significant challenge, especially for general quadrilateral elements. In this paper, we investigate the stability of the VCJH-type FR schemes for linear advection on Cartesian quadrilateral meshes and show that the schemes could become unstable under certain conditions. However, we find that the VCJH scheme recovering the DG method is stable on all Cartesian meshes. Although we restrict ourselves to Cartesian meshes in order to circumvent the algebraic complexity posed by the variation of the Jacobian matrix inside general tensor-product quadrilateral elements, our analysis offers significant insight into the possible origins of instability in the FR approach on general quadrilaterals.

Appendix: Modified Sobolev Norm

Let us reconsider the modified partial Sobolev norm used in Theorem 3.5.

$$\begin{aligned} \Vert u^D \Vert ^2_{W^{2p,2}_{\delta }} = \sum \limits _{k = 1}^N \; \left( \int \limits _{\Omega _k} J_k \left[ (u^D_k)^2 + \frac{c}{2} \left( \left( \frac{\partial ^p u^D_k}{\partial \xi ^p}\right) ^2 + \left( \frac{\partial ^p u^D_k}{\partial \eta ^p}\right) ^2 \right) + \frac{c^2}{4} \left( \frac{\partial ^{2p} u^D_k}{\partial \xi ^p \partial \eta ^p}\right) ^2 \right] d \Omega _k \right) \end{aligned}$$

(5.1)

Notice that the norm has $u^D$ in the physical domain while the derivatives are with respect to the reference coordinates. We can use the additional notation we introduced for the Cartesian mesh geometry to rewrite this completely in the physical domain as follows

$$\begin{aligned} \Vert u^D \Vert ^2_{W^{2p,2}_{\delta }}= & {} \sum \limits _{k = 1}^N \; \left( \int \limits _{\Omega _k} J_k \bigg [ (u^D_k)^2 + \frac{c}{2} \bigg (J_x^{2p} \left( \frac{\partial ^p u^D_k}{\partial x^p}\right) ^2 + J_y^{2p}\left( \frac{\partial ^p u^D_k}{\partial y^p}\right) ^2 \bigg )\right. \nonumber \\&\left. +\,\frac{c^2}{4} J_x^{2p} J_y^{2p} \left( \frac{\partial ^{2p} u^D_k}{\partial x^p \partial y^p}\right) ^2 \bigg ] d \Omega _k \right) \end{aligned}$$

(5.2)

From our analysis of stability, we can see that we are mainly interested in $c \ge 0$, since for $c < 0$, the $\varTheta _\textit{extra}$ term contributes towards instability in the simplest case of a uniform Cartesian mesh. However, as an exercise, it is interesting to investigate the range of c for which the above is a norm.

In (5.2) we write the norm completely in the physical domain. However for algebraic manipulations, it is better to write the norm completely in the reference domain. Since the norm in the domain $\Omega$ is a sum of the norms inside each element, it is sufficient to consider the norm in a single element (kth element) as follows

$$\begin{aligned} \Vert u_k^D \Vert ^2_{W^{2p,2}_{\delta }} = \int \limits _{-1}^1 \int \limits _{-1}^{1} \left[ (\hat{u}^D_k)^2 + \frac{c}{2} \left( \left( \frac{\partial ^p \hat{u}^D_k}{\partial \xi ^p}\right) ^2 + \left( \frac{\partial ^p \hat{u}^D_k}{\partial \eta ^p}\right) ^2 \right) + \frac{c^2}{4} \left( \frac{\partial ^{2p} \hat{u}^D_k}{\partial \xi ^p \partial \eta ^p}\right) ^2 \right] d \xi d \eta \end{aligned}$$

(5.3)

Till now the transformed solution $\hat{u}^D$ has been represented using the pth degree tensor-product Lagrange polynomial basis. However, we can equivalently expand our solution in a pth degree tensor product Legendre polynomial basis:

$$\begin{aligned} \hat{u}^D = \sum \limits _{i=0}^p \sum \limits _{j=0}^p L_i(\xi ) L_j(\eta ) \tilde{u}_{ij} \end{aligned}$$

(5.4)

where $\tilde{u}_{ij}$ represent the modal coefficients. This is usually referred to as the modal form while the Lagrange expansion is referred to as the nodal form of the solution. We can change from one form to the other using the corresponding Vandermonde matrix [18]. An important difference between the Lagrange and Legendre polynomials is that the nth Legendre polynomial is of degree n unlike the Lagrange polynomials which are all of degree p. Now we substitute the above expression for $\hat{u}^D$ into the norm definition (5.3). The first term can be written as follows

$$\begin{aligned} \begin{aligned} \int \limits _{-1}^1 \int \limits _{-1}^{1} (\hat{u}^D)^2 d \xi d \eta =&\sum \limits _{i=0}^p \sum \limits _{m=0}^p \sum \limits _{j=0}^p \sum \limits _{n=0}^p \tilde{u}_{ij} \tilde{u}_{mn} \int \limits _{-1}^1 \int \limits _{-1}^{1} L_i(\xi ) L_m(\xi ) L_j(\eta ) L_n(\eta ) d \xi d \eta \\ =&\sum \limits _{i=0}^p \sum \limits _{j=0}^p \tilde{u}^2_{ij} \int \limits _{-1}^1 \int \limits _{-1}^{1} L_i^2(\xi ) L_j^2(\eta ) d \xi d \eta \\&+ \sum \limits _{i=0}^p \sum \limits _{\begin{array}{c} m=0 \\ m \ne i \end{array}}^p \sum \limits _{j=0}^p \sum \limits _{n=0}^p \tilde{u}_{ij} \tilde{u}_{mn} \int \limits _{-1}^1 \int \limits _{-1}^{1} L_i(\xi ) L_m(\xi ) L_j(\eta ) L_n(\eta ) d \xi d \eta \\&+ \sum \limits _{i=0}^p \sum \limits _{j=0}^p \sum \limits _{\begin{array}{c} n=0 \\ n \ne j \end{array}}^p \tilde{u}_{ij} \tilde{u}_{in} \int \limits _{-1}^1 \int \limits _{-1}^{1} L_i^2(\xi ) L_j(\eta ) L_n(\eta ) d \xi d \eta \end{aligned} \end{aligned}$$

(5.5)

By using the orthogonality property of the Legendre polynomials, we get

$$\begin{aligned} \int \limits _{-1}^1 \int \limits _{-1}^{1} (\hat{u}^D)^2 d \xi d \eta = \sum \limits _{i=0}^p \sum \limits _{j=0}^p \bigg (\frac{2}{2i + 1}\bigg ) \bigg (\frac{2}{2j+1}\bigg ) \tilde{u}^2_{ij} \end{aligned}$$

(5.6)

Now the pth $\xi$-derivative can be written in terms of Legendre polynomials as follows

$$\begin{aligned} \frac{\partial ^p \hat{u}^D}{\partial \xi ^p} = \frac{d^p L_p(\xi )}{d \xi ^p} \sum \limits _{j=0}^p L_j(\eta ) \tilde{u}_{pj} = a_p p! \sum \limits _{j=0}^p L_j(\eta ) \tilde{u}_{pj} \end{aligned}$$

(5.7)

where one may recall that $a_p$ is the leading coefficient of $L_p$. Note that we have used the fact that the pth derivative of $L_n(\xi )$ for $n < p$ is 0 in the above expression. Therefore we have

$$\begin{aligned} \int \limits _{-1}^1 \int \limits _{-1}^{1} \left( \frac{\partial ^p \hat{u}^D}{\partial \xi ^p}\right) ^2 d \xi d \eta = 2 (a_p p!)^2 \sum \limits _{j=0}^p \int \limits _{-1}^{1} L^2_j(\eta ) \tilde{u}^2_{pj} d \eta = 2 (a_p p!)^2 \sum \limits _{j=0}^p \bigg (\frac{2}{2j+1}\bigg ) \tilde{u}^2_{pj} \end{aligned}$$

(5.8)

Similarly we have the pth $\eta$-derivative

$$\begin{aligned} \int \limits _{-1}^1 \int \limits _{-1}^{1} \left( \frac{\partial ^p \hat{u}^D}{\partial \eta ^p}\right) ^2 d \xi d \eta = 2 (a_p p!)^2 \sum \limits _{i=0}^p \bigg (\frac{2}{2i+1}\bigg ) \tilde{u}^2_{ip} \end{aligned}$$

(5.9)

Now we consider the last term of the norm

$$\begin{aligned} \frac{\partial ^{2p} \hat{u}^D}{\partial \xi ^p \partial \eta ^p} = \frac{d^p L_p(\xi )}{d \xi ^p} \frac{d^p L_p(\eta )}{d \eta ^p} \tilde{u}_{pp} = (a_p p!)^2 \tilde{u}_{pp} \end{aligned}$$

(5.10)

Therefore

$$\begin{aligned} \int \limits _{-1}^1 \int \limits _{-1}^{1} \bigg (\frac{\partial ^{2p} \hat{u}^D}{\partial \xi ^p \partial \eta ^p} \bigg )^2 d\xi d\eta = 4 (a_p p!)^4 \tilde{u}^2_{pp} \end{aligned}$$

(5.11)

From (5.6), (5.8), (5.9) and (5.11), we can see that the norm inside the kth element can be written as

$$\begin{aligned} \begin{aligned} \Vert u_k^D \Vert ^2_{W^{2p,2}_{\delta }} =&\sum \limits _{i=0}^{p-1} \sum \limits _{j=0}^{p-1} \bigg (\frac{2}{2i + 1}\bigg ) \bigg (\frac{2}{2j+1}\bigg ) \tilde{u}^2_{ij} \\&+ \sum \limits _{j=0}^{p-1} \bigg [ \bigg (\frac{2}{2p + 1}\bigg ) \bigg (\frac{2}{2j+1}\bigg ) + \frac{c}{2}\; 2 (a_p p!)^2 \bigg (\frac{2}{2j+1}\bigg ) \bigg ] \tilde{u}^2_{pj} \\&+ \sum \limits _{i=0}^{p-1} \bigg [ \bigg (\frac{2}{2p + 1}\bigg ) \bigg (\frac{2}{2i+1}\bigg ) + \frac{c}{2}\; 2 (a_p p!)^2 \bigg (\frac{2}{2i+1}\bigg ) \bigg ] \tilde{u}^2_{ip} \\&+ \bigg [ \frac{4}{(2p+1)^2} + \frac{2 (a_p p!)^2}{2p+1} c + (a_p p!)^4 c^2\bigg ] \tilde{u}^2_{pp} \end{aligned} \end{aligned}$$

(5.12)

In order for this to be a norm, we need the co-efficients of each $\tilde{u}_{ij}$ have to be non negative. Therefore we have the following 2 conditions,

$$\begin{aligned}&\displaystyle \frac{1}{2p+1} + \frac{c}{2} (a_p p!)^2 \ge 0 \Longrightarrow c \ge \frac{-2}{(2p+1) (a_p p!)^2} \end{aligned}$$

(5.13)

$$\begin{aligned}&\displaystyle \frac{4}{(2p+1)^2} + \frac{2 (a_p p!)^2}{2p+1} c + (a_p p!)^4 c^2 \ge 0 \end{aligned}$$

(5.14)

The first condition is the same as the one obtained for the 1D case in [6]. The LHS of the second condition is a convex quadratic with a negative discriminant, therefore condition 2, i.e., (5.14) is always satisfied. Therefore, the condition on c for (5.3) to be a norm is the same as obtained in 1D.