Energy Stable Flux Reconstruction Schemes for Advection–Diffusion Problems on Tetrahedra

Abstract

The flux reconstruction (FR) methodology provides a unifying description of many high-order schemes, including a particular discontinuous Galerkin (DG) scheme and several spectral difference (SD) schemes. In addition, the FR methodology has been used to generate new classes of high-order schemes, including the recently discovered ‘energy stable’ FR schemes. These schemes, which are often referred to as VCJH (Vincent–Castonguay–Jameson–Huynh) schemes, are provably stable for linear advection–diffusion problems in 1D and on triangular elements. The VCJH schemes have been successfully applied to a wide variety of problems in 1D and 2D, ranging from linear advection–diffusion problems, to fluid mechanics problems requiring the solution of the compressible Navier–Stokes equations. Based on the results of these numerical experiments, it has been shown that certain VCJH schemes maintain the expected order of spatial accuracy and possess explicit time-step limits which rival those of the collocation-based nodal DG scheme. However, it remained to be seen whether the VCJH schemes could be extended to 3D on tetrahedral elements, enabling their convenient application to the complex geometries that arise in many real-world problems. For the first time, this article presents an extension of the VCJH schemes to tetrahedral elements. This work provides a formal proof of the stability of the new schemes and assesses their performance via numerical experiments on model problems.

Appendices

Appendix A: Constructing the Energy Stable (VCJH) Correction Fields

In this section, a procedure will be presented for constructing energy stable (VCJH) correction fields $\phi _{f,l}$ and $\psi _{f,l}$ that satisfy Eqs. (43) and (53), respectively.

1.1 Preliminaries

Recall that $\phi _{f,l} \equiv \hat{\nabla } \cdot \mathbf{h}_{f,l}$ and $\psi _{f,l} \equiv \hat{\nabla } \cdot \mathbf{g}_{f,l}$. In accordance with these definitions, it may appear natural to first construct precise expressions for the correction functions $\mathbf{h}_{f,l}$ and $\mathbf{g}_{f,l}$, and thereafter apply the divergence operator $\hat{\nabla }$ to these expressions in order to obtain $\phi _{f,l}$ and $\psi _{f,l}$. However, this is not the best approach, as $\mathbf{h}_{f,l}$ and $\mathbf{g}_{f,l}$ may not be unique. In particular, for $p > 1$ there are an unlimited number of correction functions which have the same divergence (or effectively, the same correction field). In addition, the implementation of the VCJH approach only requires the precise formulation of the correction fields $\phi _{f,l}$ and $\psi _{f,l}$, and not of the correction functions themselves. Specifically, the VCJH approach only requires definitions of the normal components of the correction functions ($\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}}$ and $\mathbf{g}_{f,l} \cdot \hat{\mathbf{n}}$) which are given by Eqs. (23) and (24). In what follows, it will be shown that Eqs. (23), (24), (43), and (53) are sufficient for constructing precise definitions of $\phi _{f,l}$ and $\psi _{f,l}$.

1.2 Constructing the Correction Fields $\phi _{f,l}$

One may arrive at a procedure for constructing the correction fields $\phi _{f,l}$ by manipulating Eq. (43) and utilizing the definition of $\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}}$ (Eqs. (23) and (24)). Recall from section 3, Lemma 3.1, that the correction functions $\mathbf{h}_{f,l}$ and correction fields $\phi _{f,l}$ are required to satisfy Eq. (43), which can be expanded as follows

$$\begin{aligned}&\sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( \hat{\nabla } \cdot \hat{\mathbf{f}}^{C} \right) \hat{D}^{(p,v,w)} \left( \hat{\ell }_{j} \right) \nonumber \\&\quad = \int \limits _{\varvec{\Gamma }_S} \left( \hat{\mathbf{f}}^{C} \cdot \hat{\mathbf{n}} \right) \hat{\ell }_{j} \; d\varvec{\Gamma }_S - \int \limits _{\varvec{\Omega }_S} \left( \hat{\nabla } \cdot \hat{\mathbf{f}}^{C} \right) \hat{\ell }_{j} \; d\varvec{\Omega }_S \nonumber \\&\quad = \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( \sum _{f = 1}^{N_{fe}} \sum _{l = 1}^{N_{fp}} \left[ \left( \hat{\mathbf{f}}^{\star }_{f,l} - \hat{\mathbf{f}}^{D}_{f,l} \right) \cdot \hat{\mathbf{n}}_{f,l} \right] \, \phi _{f,l} \right) \hat{D}^{(p,v,w)} \left( \hat{\ell }_{j} \right) \nonumber \\&\quad = \int \limits _{\varvec{\Gamma }_S} \left( \sum _{f = 1}^{N_{fe}} \sum _{l = 1}^{N_{fp}} \left[ \left( \hat{\mathbf{f}}^{\star }_{f,l} - \hat{\mathbf{f}}^{D}_{f,l} \right) \cdot \hat{\mathbf{n}}_{f,l} \right] \, \mathbf{h}_{f,l} \cdot \hat{\mathbf{n}} \right) \hat{\ell }_{j} \; d\varvec{\Gamma }_S \nonumber \\&\qquad - \int \limits _{\varvec{\Omega }_S} \left( \sum _{f = 1}^{N_{fe}} \sum _{l = 1}^{N_{fp}} \left[ \left( \hat{\mathbf{f}}^{\star }_{f,l} - \hat{\mathbf{f}}^{D}_{f,l} \right) \cdot \hat{\mathbf{n}}_{f,l} \right] \, \phi _{f,l} \right) \hat{\ell }_{j} \; d\varvec{\Omega }_S. \end{aligned}$$

(144)

Upon rearranging and simplifying Eq. (144), one obtains the following

$$\begin{aligned} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( \phi _{f,l} \right) \hat{D}^{(p,v,w)} \left( \hat{\ell }_{j} \right) = \int \limits _{\varvec{\Gamma }_S} \left( \mathbf{h}_{f,l} \cdot \hat{\mathbf{n}} \right) \hat{\ell }_{j} \; d\varvec{\Gamma }_S - \int \limits _{\varvec{\Omega }_S} \left( \phi _{f,l} \right) \hat{\ell }_{j} \; d\varvec{\Omega }_S. \nonumber \\ \end{aligned}$$

(145)

Here, each basis function $\hat{\ell }_{j}$ is equal to the product of $J_k$ and $\ell _{j}$, each normal component of a correction function $(\mathbf{h}_{f,l} \cdot \hat{\mathbf{n}})$ is defined by Eq. (23) (with $\mathbf{h}_{f,l}$ in place of $\mathbf{g}_{f,l}$), and each correction field $\phi _{f,l}$ remains to be determined. In order to begin the process of computing $\phi _{f,l}$, note that it is a degree $p$ polynomial function which can be expressed as follows

$$\begin{aligned} \phi _{f,l} = \sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, L_{\jmath } \left( \hat{\mathbf{x}} \right) , \end{aligned}$$

(146)

where each $\sigma _{\jmath }$ is a constant coefficient (yet to be computed) and each basis function $L_{\jmath } \left( \hat{\mathbf{x}} \right)$ is a member of an orthonormal basis of degree $p$ on the reference element $\varvec{\Omega }_S$ defined as follows

$$\begin{aligned} L_{\jmath } \left( \hat{\mathbf{x}} \right)&= \sqrt{8} P_u \left( \hat{a} \right) P_v^{\left( 2u + 1,0\right) } \left( \hat{b} \right) \left( 1- \hat{b} \right) ^u P_w^{\left( 2u+2v+2,0\right) } \left( \hat{c} \right) \left( 1- \hat{c} \right) ^{u+v}, \\ \jmath&= 1 + \frac{\left( 11 + 12p + 3p^2 \right) u}{6} + \frac{\left( 2p + 3 \right) v}{2} + w - \frac{\left( 2 + p \right) u^2}{2} - uv - \frac{v^2}{2} + \frac{u^3}{6},\nonumber \\&0 \le u, \quad 0 \le v, \quad 0 \le w, \qquad u + v+ w \le p,\nonumber \end{aligned}$$

(147)

where $P_n^{\left( \alpha , \beta \right) }$ is the normalized $n^{th}$ order Jacobi polynomial (as defined in [19]), and where

$$\begin{aligned} \hat{a} = -2 \frac{\left( 1 + \hat{x} \right) }{\hat{y} + \hat{z}} -1, \qquad \hat{b} = 2 \frac{\left( 1 + \hat{y}\right) }{1 - \hat{z}} - 1, \qquad \hat{c} = \hat{z}. \end{aligned}$$

(148)

In addition, note that each function $\hat{\ell }_{j}$ in Eq. (145) is a degree $p$ polynomial which can be expressed as follows

$$\begin{aligned} \hat{\ell }_j = \sum _{\imath = 1}^{N_p} \gamma _{\imath } \, L_{\imath } \left( \hat{\mathbf{x}} \right) , \end{aligned}$$

(149)

where each constant coefficient $\gamma _{\imath }$ is a known quantity because each function $\hat{\ell }_j$ is a known quantity.

On substituting Eqs. (146) and (149) into Eq. (145), one obtains

$$\begin{aligned}&\sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( \sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, L_{\jmath } \right) \hat{D}^{(p,v,w)} \left( \sum _{\imath = 1}^{N_p} \gamma _{\imath } \, L_{\imath } \right) \nonumber \\&\quad = \int \limits _{\varvec{\Gamma }_S} \left( \mathbf{h}_{f,l} \cdot \hat{\mathbf{n}} \right) \left( \sum _{\imath = 1}^{N_p} \gamma _{\imath } \, L_{\imath } \right) \; d\varvec{\Gamma }_S - \int \limits _{\varvec{\Omega }_S} \left( \sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, L_{\jmath } \right) \left( \sum _{\imath = 1}^{N_p} \gamma _{\imath } \, L_{\imath } \right) \; d\varvec{\Omega }_S, \end{aligned}$$

(150)

or equivalently

$$\begin{aligned}&\sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( L_{\jmath } \right) \hat{D}^{(p,v,w)} \left( L_{\imath } \right) \nonumber \\&\quad = \int \limits _{\varvec{\Gamma }_S} \left( \mathbf{h}_{f,l} \cdot \hat{\mathbf{n}} \right) L_{\imath } \; d\varvec{\Gamma }_S - \sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, \int \limits _{\varvec{\Omega }_S} L_{\jmath } \; L_{\imath } \; d\varvec{\Omega }_S. \end{aligned}$$

(151)

Upon noting that $L_{\imath }$ and $L_{\jmath }$ are orthonormal polynomials, one may derive the following expression from Eq. (151)

$$\begin{aligned} \sum _{\jmath = 1}^{N_p} \sigma _{\jmath } \, \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \, \hat{D}^{(p,v,w)} \left( L_{\jmath } \right) \hat{D}^{(p,v,w)} \left( L_{\imath } \right) = -\sigma _{\imath } + \int \limits _{\varvec{\Gamma }_S} \left( \mathbf{h}_{f,l} \cdot \hat{\mathbf{n}} \right) L_{\imath } \; d\varvec{\Gamma }_S. \qquad \end{aligned}$$

(152)

Equation (152) holds for $\imath = 1, \ldots , N_p$ and provides a system of $N_p$ equations for the $N_p$ unknown coefficients $\sigma _{\jmath }$. Together, Eqs. (152) and (146) completely define $\phi _{f,l}$.

1.3 Constructing the Correction Fields $\psi _{f,l}$

In following the approach of the previous section, one may arrive at a procedure for constructing the correction fields $\psi _{f,l}$ by manipulating Eq. (53) and utilizing the definition of $\mathbf{g}_{f,l} \cdot \hat{\mathbf{n}}$ (Eqs. (23) and (24)). Once these manipulations (which are omitted for the sake of brevity) are performed, one obtains the following formula for each field $\psi _{f,l}$

$$\begin{aligned} \psi _{f,l} = \sum _{\jmath = 1}^{N_p} \xi _{\jmath } \, L_{\jmath } \left( \hat{\mathbf{x}} \right) , \end{aligned}$$

(153)

where the $N_p$ unknown coefficients $\xi _{\jmath }$ can be obtained from the following system of $\imath = 1, \ldots , N_p$ equations

$$\begin{aligned} \sum _{\jmath = 1}^{N_p} \xi _{\jmath } \, \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} \kappa _{vw} \, \hat{D}^{(p,v,w)} \left( L_{\jmath } \right) \hat{D}^{(p,v,w)} \left( L_{\imath } \right) = -\xi _{\imath } + \int \limits _{\varvec{\Gamma }_S} \left( \mathbf{g}_{f,l} \cdot \hat{\mathbf{n}} \right) L_{\imath } \; d\varvec{\Gamma }_S. \qquad \end{aligned}$$

(154)

Appendix B: Energy Stable (VCJH) Filter Matrices

This section discusses the procedure for forming the energy stable (VCJH) filter matrices and examines the resulting structure of the filter matrices.

1.1 Procedure for Forming the Filter Matrices

The filters $\mathbf{F}_1$ and $\mathbf{F}_2$ are obtained from the matrices $\mathbf{M}^{k}, \mathcal M ^{k}, \mathbf{K}^{k}$, and $\mathcal K ^{k}$ via Eq. (142). The mass matrices $\mathbf{M}^{k}$ and $\mathcal M ^{k}$ can be constructed from inner products of the nodal basis functions $\ell _j \left( \hat{\mathbf{x}} \right)$ in accordance with Eq. (70). However, it is often difficult to compute the nodal basis functions $\ell _j \left( \hat{\mathbf{x}} \right)$ directly, and one is often better served by using the orthonormal basis functions $L_{j} \left( \hat{\mathbf{x}} \right)$ in their place. In particular, one may use the orthonormal basis functions to define a ‘Vandermonde matrix’ $\mathbf{V}^k$ which has the following entries

$$\begin{aligned}{}[\mathbf{V}^k]_{ij} = L_j \left( \hat{\mathbf{x}}_i \right) , \end{aligned}$$

(155)

and next, one may use $\mathbf{V}^k$ to compute the mass matrices as follows

$$\begin{aligned} \mathbf{M}^{k}&= J_k \left( \mathbf{V}^k \left( \mathbf{V}^k \right) ^T \right) ^{-1}, \end{aligned}$$

(156)

$$\begin{aligned} \mathcal M ^k&= J_k \begin{bmatrix} \left( \mathbf{V}^k \left( \mathbf{V}^k \right) ^T \right) ^{-1}&\mathbf 0&\mathbf 0 \\ \mathbf{0}&\left( \mathbf{V}^k \left( \mathbf{V}^k \right) ^T \right) ^{-1}&\mathbf 0 \\ \mathbf{0}&\mathbf{0}&\left( \mathbf{V}^k \left( \mathbf{V}^k \right) ^T \right) ^{-1} \end{bmatrix}. \end{aligned}$$

(157)

Note that the derivation of Eqs. (156) and (157) is discussed in detail in [19].

In a similar fashion, one may construct expressions for $\mathbf{K}^{k}$ and $\mathcal K ^{k}$ in terms of the orthonormal basis functions. Recall that Eqs. (76) and (77) provide expressions for $\mathbf{K}^{k}$ and $\mathcal K ^{k}$ in terms of the matrices $\left( \hat{\mathbf{D}}^{(p,v,w)} \right) ^{T} \hat{\mathbf{D}}^{(p,v,w)}$, where each $\hat{\mathbf{D}}^{(p,v,w)}$ is a matrix which has the following entries

$$\begin{aligned}{}[\hat{\mathbf{D}}^{(p,v,w)}]_{ij} = \hat{D}^{(p,v,w)} \left( \ell _j \left( \hat{\mathbf{x}} \right) \right) \bigg |_{\hat{\mathbf{x}}_i}. \end{aligned}$$

(158)

One may use the orthonormal basis functions $L_{j} \left( \hat{\mathbf{x}} \right)$ to construct similar matrices $\hat{\hat{\mathbf{D}}}^{(p,v,w)}$ which have the following entries

$$\begin{aligned}{}[\hat{\hat{\mathbf{D}}}^{(p,v,w)}]_{ij} = \hat{D}^{(p,v,w)} \left( L_j \left( \hat{\mathbf{x}} \right) \right) \bigg |_{\hat{\mathbf{x}}_i}. \end{aligned}$$

(159)

One may then relate each $\hat{\hat{\mathbf{D}}}^{(p,v,w)}$ to each $\hat{\mathbf{D}}^{(p,v,w)}$ by using the following identity

$$\begin{aligned} \hat{\mathbf{D}}^{(p,v,w)} = \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k \right) ^{-1}. \end{aligned}$$

(160)

Note that the derivation of Eq. (160) is discussed in detail in [19].

Upon substituting Eq. (160) into Eqs. (76) and (77), one obtains the following expressions for $\mathbf{K}^{k}$ and $\mathcal K ^{k}$

$$\begin{aligned} \mathbf{K}^k&= \frac{J_k}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \left( \mathbf{V}^k \right) ^{-T} \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k \right) ^{-1}, \end{aligned}$$

(161)

$$\begin{aligned} \mathcal K ^k&= \frac{J_k}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} \kappa _{vw} \begin{bmatrix} \left( \mathbf{V}^k \right) ^{-T} \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \left( \mathbf{V}^k\right) ^{-T} \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \left( \mathbf{V}^k\right) ^{-T} \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \end{bmatrix}. \end{aligned}$$

(162)

Next, on substituting Eqs. (161), (162), (156), and (157) into Eq. (142), one obtains the following expressions for $\mathbf{F}_1$ and $\mathbf{F}_2$

$$\begin{aligned} \mathbf{F}_1&= \left( \mathbf{I}+ \frac{1}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \; \mathbf{V}^k \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k \right) ^{-1} \right) ^{-1} \nonumber \\&= \mathbf{V}^k \, \hat{\hat{ \mathbf{F}}}_1 \left( \mathbf{V}^k \right) ^{-1}, \end{aligned}$$

(163)

where

$$\begin{aligned} \hat{\hat{ \mathbf{F}}}_1 = \left( \mathbf{I}+ \frac{1}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} c_{vw} \; \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \right) ^{-1}, \end{aligned}$$

(164)

and

$$\begin{aligned} \nonumber \mathbf{F}_{2}&= \left( \mathcal I + \frac{1}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} \kappa _{vw} \begin{bmatrix} \mathbf{V}^k \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \mathbf{V}^k \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \mathbf{V}^k \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \, \left( \mathbf{V}^k\right) ^{-1} \end{bmatrix} \right) ^{-1} \\&= \mathcal V ^k \, \hat{\hat{\mathbf{F}}}_2 \left( \mathcal V ^k \right) ^{-1}, \end{aligned}$$

(165)

where

$$\begin{aligned} \hat{\hat{\mathbf{F}}}_2 = \left( \mathcal I + \frac{1}{N_p} \sum _{v = 1}^{p+1} \sum _{w = 1}^{v} \kappa _{vw} \begin{bmatrix} \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \qquad \qquad \qquad \mathbf 0 \\ \mathbf 0 \qquad \qquad \qquad \mathbf 0 \qquad \qquad \qquad \left( \hat{\hat{\mathbf{D}}}^{(p,v,w)} \right) ^{T} \hat{\hat{\mathbf{D}}}^{(p,v,w)} \end{bmatrix} \right) ^{-1}, \end{aligned}$$

(166)

and

$$\begin{aligned} \mathcal V ^k = \begin{bmatrix} \mathbf{V}^k&\mathbf 0&\mathbf 0 \\ \mathbf 0&\mathbf{V}^k&\mathbf 0 \\ \mathbf 0&\mathbf 0&\mathbf{V}^k \end{bmatrix}. \end{aligned}$$

(167)

Note that Eqs. (164) and (166) define new filtering matrices $\hat{\hat{\mathbf{F}}}_1$ and $\hat{\hat{\mathbf{F}}}_2$. These matrices can be viewed as filters which act on the orthonormal basis (in contrast to $\mathbf{F}_1$ and $\mathbf{F}_2$ which can be viewed as filters which act on the nodal basis). Conveniently, the two sets of filters are related via left and right multiplication by the Vandermonde matrix and its inverse (as shown in Eqs. (163) and (165)).

Now, having established a method for constructing $\mathbf{F}_1$ and $\mathbf{F}_2$ using the orthornormal basis (via $\hat{\hat{\mathbf{F}}}_1$ and $\hat{\hat{\mathbf{F}}}_2$), one may obtain insights into the overall effects of the filtering process by examining the sparsity patterns of $\hat{\hat{\mathbf{F}}}_1$ and $\hat{\hat{\mathbf{F}}}_2$.

1.2 Sparsity Patterns of the Filter Matrices

On evaluating Eq. (164), one finds that $\hat{\hat{\mathbf{F}}}_1$ has the following block structure

$$\begin{aligned} \hat{\hat{\mathbf{F}_1}} = \begin{bmatrix} \mathbf{I}_{1}^{B}&\mathbf 0 \\ \mathbf 0&\mathbf{F}_{1}^{B} \end{bmatrix}, \end{aligned}$$

(168)

where $\mathbf{I}_{1}^{B} \in \mathbb R ^{N_p^{\ell } \times N_p^{\ell }}$ is an identity matrix, $\mathbf{F}_{1}^{B} \in \mathbb R ^{N_p^{u} \times N_p^{u}}$ is a dense matrix of filtering coefficients, $N_p^{\ell } = N_p - N_p^{u}$ is the number of orthonormal basis functions of degree $\le \left( p - 1\right)$, and $N_p^{u} = \frac{1}{2} \left( p + 1\right) \left( p + 2 \right)$ is the number of orthonormal basis functions of degree $p$. The structure of $\hat{\hat{\mathbf{F}_1}}$ ensures that only the degree $p$ orthonormal basis functions are effected by the filtering matrix. All basis functions of degree $\le \left( p-1\right)$ are multiplied by the identity matrix and remain unaffected.

Similarly, on evaluating Eq. (166), one finds that $\hat{\hat{\mathbf{F}}}_2$ has the following structure

$$\begin{aligned} \hat{\hat{\mathbf{F}_2}} = \left[ \begin{array}{llllll} \mathbf{I}_{2}^B &{} \mathbf 0 &{} \cdots &{} \cdots &{} \cdots &{} \mathbf 0 \\ \mathbf 0 &{} \mathbf{F}_{2}^B &{} \ddots &{} &{} &{} \vdots \\ \vdots &{} \ddots &{} \mathbf{I}_{2}^B &{} \mathbf 0 &{} &{} \vdots \\ \vdots &{} &{} \mathbf 0 &{} \mathbf{F}_{2}^B &{} \ddots &{} \vdots \\ \vdots &{} &{} &{} \ddots &{} \mathbf{I}_{2}^B &{} \mathbf 0 \\ \mathbf 0 &{} \cdots &{} \cdots &{} \cdots &{} \mathbf 0 &{} \mathbf{F}_{2}^B\\ \end{array}\right] , \end{aligned}$$

(169)

where $\mathbf{I}_{2}^{B} \in \mathbb R ^{N_p^{\ell } \times N_p^{\ell }}$ is an identity matrix and $\mathbf{F}_{2}^{B} \in \mathbb R ^{N_p^{u} \times N_p^{u}}$ is a dense matrix of filtering coefficients. The structure of $\hat{\hat{\mathbf{F}_2}}$ is similar to the structure of $\hat{\hat{\mathbf{F}_1}}$, as the structure of $\hat{\hat{\mathbf{F}}}_2$ also ensures that only the degree $p$ orthonormal basis functions are effected by the filtering matrix.

In summary, the filters $\mathbf{F}_1$ and $\mathbf{F}_2$ are related (via the Vandermonde matrix) to the filters $\hat{\hat{\mathbf{F}_1}}$ and $\hat{\hat{\mathbf{F}_2}}$ which act on the orthonormal basis functions $L_j (\hat{\mathbf{x}})$, and effect only the highest (degree $p$) modes of the residual and the auxiliary variable, respectively.