A Polynomial Spectral Calculus for Analysis of DG Spectral Element Methods

TL;DR

This paper introduces a polynomial spectral calculus based on Legendre-Gauss-Lobatto quadrature, simplifying the analysis of multidimensional discontinuous Galerkin spectral element methods.

Contribution

It develops a new polynomial spectral calculus leveraging summation by parts, aiding in the analysis of DG spectral element methods.

Findings

Simplifies multidimensional DG spectral element analysis

Provides a new mathematical tool based on polynomial calculus

Enhances understanding of spectral element method stability

Abstract

We introduce a polynomial spectral calculus that follows from the summation by parts property of the Legendre-Gauss-Lobatto quadrature. We use the calculus to simplify the analysis of two multidimensional discontinuous Galerkin spectral element approximations.

Equations

$${\mathbf{u}}_t+\nabla \cdot \vec{\mathbf{f}}=0,$$

$$\vec{\mathbf{f}}\left(\mathbf{u}\right)=\sum _{m=1}^3{\mathcal{A}}^{(m)}\left(\vec{x}\right)\mathbf{u}{\hat{x}}_m\equiv \vec{\mathcal{A}}\mathbf{u}$$

$${\mathcal{A}}^{(m)}={\left({\mathcal{A}}^{(m)}\right)}^T\text{and}\sum _{m=1}^3{\alpha }_m{\mathcal{A}}^{(m)}=\mathcal{R}\left(\vec{\alpha }\right)\Lambda \left(\vec{\alpha }\right){\mathcal{R}}^{-1}\left(\vec{\alpha }\right)$$

$${\Vert \nabla \cdot \vec{\mathcal{A}}\Vert }_2<\infty ,$$

$${\mathbf{u}}_t+\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u}+\vec{\mathcal{A}}\cdot \nabla \mathbf{u}=0.$$

$$\left(\mathbf{u},\mathbf{v}\right)={\int }_{\Omega }{\mathbf{u}}^T\mathbf{v}dxdydz,\Vert \mathbf{u}\Vert =\sqrt{\left(\mathbf{u},\mathbf{u}\right)}.$$

$$\left({\mathbf{u}}_t,\mathbf{\varphi }\right)+\frac{1}{2}\left(\nabla \cdot \vec{\mathbf{f}},\mathbf{\varphi }\right)+\frac{1}{2}\left\{\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u},\mathbf{\varphi }\right)+\left(\vec{\mathcal{A}}\cdot \nabla \mathbf{u},\mathbf{\varphi }\right)\right\}=0.$$

$${\int }_{\Omega }{\mathbf{u}}^T\nabla \cdot \vec{\mathbf{f}}dxdydz={\int }_{\partial \Omega }{\mathbf{u}}^T\vec{\mathbf{f}}\cdot \hat{n}dS-{\int }_{\Omega }{\left(\nabla \mathbf{u}\right)}^T\cdot \vec{\mathbf{f}}dxdydz,$$

$$\left(\mathbf{u},\nabla \cdot \vec{\mathbf{f}}\right)={\int }_{\partial \Omega }{\mathbf{u}}^T\vec{\mathbf{f}}\cdot \hat{n}dS-\left(\nabla \mathbf{u},\vec{\mathbf{f}}\right).$$

$$\left({\mathbf{u}}_t,\mathbf{\varphi }\right)+{\int }_{\partial \Omega }\vec{\mathbf{f}}\cdot \hat{n}\mathbf{\varphi }dS-\frac{1}{2}\left(\vec{\mathbf{f}},\nabla \mathbf{\varphi }\right)+\frac{1}{2}\left\{\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u},\varphi \right)-\left(\mathbf{u},\nabla \cdot \left(\vec{\mathcal{A}}\mathbf{\varphi }\right)\right)\right\}=0.$$

$$\left({\mathbf{u}}_t,\varphi \right)+\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{c},\varphi \right)=0,$$

$$\frac{d}{dt}{\int }_{\Omega }\mathbf{u}dxdydz=-{\int }_{\partial \Omega }\vec{\mathbf{f}}\cdot \hat{n}dS.$$

$$\left(\nabla \cdot \vec{\mathbf{f}}+\vec{\mathcal{A}}\cdot \nabla \mathbf{u},\mathbf{u}\right)={\int }_{\Omega }\nabla \cdot \left({\mathbf{u}}^T\vec{\mathcal{A}}\mathbf{u}\right)d\vec{x}.$$

$$\frac{d}{dt}{\left|\left|\mathbf{u}\right|\right|}^2+\underset{\Omega }{\int }\nabla \cdot \left({\mathbf{u}}^T\vec{\mathcal{A}}\mathbf{u}\right)d\vec{x}+\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u},\mathbf{u}\right)=0.$$

$$\frac{d}{dt}{\left|\left|\mathbf{u}\right|\right|}^2+\underset{\partial \Omega }{\int }{\mathbf{u}}^T\vec{\mathcal{A}}\cdot \hat{n}\mathbf{u}dS=-\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u},\mathbf{u}\right).$$

$$-\left(\left(\nabla \cdot \vec{\mathcal{A}}\right)\mathbf{u},\mathbf{u}\right)⩽\underset{\Omega }{\max }{\left|\left|\nabla \cdot \vec{\mathcal{A}}\right|\right|}_2{\left|\left|\mathbf{u}\right|\right|}^2\equiv 2\gamma {\left|\left|\mathbf{u}\right|\right|}^2$$

$$\frac{d}{dt}\left(e^{-2\gamma t}{\left|\left|\mathbf{u}\right|\right|}^2\right)\le e^{-2\gamma t}\underset{\partial \Omega }{\int }{\mathbf{u}}^T\vec{\mathcal{A}}\cdot \hat{n}\mathbf{u}dS.$$

$${\left|\left|\mathbf{u}(T)\right|\right|}^2\le e^{2\gamma T}{\left|\left|\mathbf{u}(0)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{e}^{2\gamma (T-t)}{\mathbf{u}}^T\vec{\mathcal{A}}\cdot \hat{n}\mathbf{u}dSdt.$$

$$\vec{\mathcal{A}}\cdot \hat{n}=\sum _{m=1}^3{\mathcal{A}}^{(m)}{\hat{n}}_m=\mathcal{R}\Lambda {\mathcal{R}}^{-1}=\mathcal{P}{\Lambda }^{+}{\mathcal{R}}^{-1}+\mathcal{R}{\Lambda }^{-}{\mathcal{R}}^{-1}\equiv {\mathcal{A}}^{+}+{\mathcal{A}}^{-},$$

$$\begin{array}{rl}{\left|\left|\mathbf{u}(T)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{\mathbf{u}}^T{\mathcal{A}}^{+}\mathbf{u}dSdt& \le e^{2\gamma T}{\left|\left|\mathbf{u}(0)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{e}^{2\gamma (T-t)}{\mathbf{g}}^T\left|{\mathcal{A}}^{-}\right|\mathbf{g}dSdt\\ & \le e^{2\gamma T}\left\{{\left|\left|\mathbf{u}(0)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{\mathbf{g}}^T\left|{\mathcal{A}}^{-}\right|\mathbf{g}dSdt\right\}.\end{array}$$

$${\left|\left|\mathbf{u}(T)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{\mathbf{u}}^T{\mathcal{A}}^{+}\mathbf{u}dSdt\le {\left|\left|\mathbf{u}(0)\right|\right|}^2+{\int }_0^T{\int }_{\partial \Omega }{\mathbf{g}}^T\left|{\mathcal{A}}^{-}\right|\mathbf{g}dSdt.$$

$$U=\sum _{i,j,k=0}^N{U}_{ijk}{\ell }_i(\xi ){\ell }_j(\eta ){\ell }_k(\zeta ),$$

$${\ell }_l\left(s\right)=\prod _{i=0;i\ne l}^N\frac{s-s_i}{s_l-s_i}$$

$${\mathbb{I}}^N\left(u\right)=\sum _{i,j,k=0}^N{u}_{ijk}{\ell }_i(\xi ){\ell }_j(\eta ){\ell }_k(\zeta ).$$

$${\frac{\partial U}{\partial \xi }|}_{nml}=\sum _{i,j,k=0}^N{U}_{ijk}{\ell }_i^{\prime }({\xi }_n){\ell }_j({\eta }_m){\ell }_k({\zeta }_l)=\sum _{i=0}^N{U}_{iml}{{\ell }^{\prime }}_i({\xi }_n)\equiv \sum _{i=0}^N{U}_{iml}{\mathcal{D}}_{ni},$$

$$\begin{array}{c}{\nabla U|}_{nml}=\sum _{i=0}^N{U}_{iml}{\mathcal{D}}_{ni}\hat{\xi }+\sum _{j=0}^N{U}_{njl}{\mathcal{D}}_{mj}\hat{\eta }+\sum _{k=0}^N{U}_{nmk}{\mathcal{D}}_{lk}\hat{\zeta },\text{hfill}\\ {\nabla \cdot \vec{F}|}_{nml}=\sum _{i=0}^N{F}_{iml}^{(1)}{\mathcal{D}}_{ni}+\sum _{j=0}^N{F}_{njl}^{(2)}{\mathcal{D}}_{mj}+\sum _{k=0}^N{F}_{nmk}^{(3)}{\mathcal{D}}_{lk}.\text{hfill}\end{array}$$

$${\int }_{-1}^{1}gds\approx \sum _{i=0}^{N}g\left(s_i\right){\omega }_i\equiv {\int }_{N}gds,$$

$${\int }_{E,N}gd\xi d\eta d\zeta \equiv \sum _{i,j,k=0}^N{g}_{ijk}{\omega }_{ijk},$$

$$\begin{array}{rl}{\int }_{\partial E,N}\vec{g}\cdot \hat{n}dS& =\sum _{i,j=0}^N{{\omega }_{ij}{g}^{(1)}\left(\xi ,{\eta }_i,{\zeta }_j\right)|}_{\xi =-1}^1+\sum _{i,j=0}^N{{\omega }_{ij}{g}^{(2)}\left({\xi }_i,\eta ,{\zeta }_j\right)|}_{\eta =-1}^1+\sum _{i,j=0}^N{{\omega }_{ij}{g}^{(3)}\left({\xi }_i,{\eta }_j,\zeta \right)|}_{\zeta =-1}^1\\ & \equiv {\int }_N{g^{(1)}d\eta d\zeta |}_{\xi =-1}^1+{\int }_N{g^{(2)}d\xi d\zeta |}_{\eta =-1}^1+{\int }_N{g^{(3)}d\xi d\eta |}_{\zeta =-1}^1.\end{array}$$

$${\left(f,g\right)}_{E,N}={\int }_{E,N}fgd\xi d\eta d\zeta \equiv \sum _{i,j,k=0}^N{f}_{ijk}{g}_{ijk}{\omega }_{ijk},{\Vert f\Vert }_{E,N}=\sqrt{{\left(f,f\right)}_{E,N}}.$$

Full text

1 Introduction

The discontinuous Galerkin Spectral Element Method (DGSEM) introduced by Black [4],[5] has the desired properties of spectral accuracy, geometric flexibility, and excellent phase and dissipation properties [22],[10]. Spectral accuracy comes from the use of high order polynomial approximations to the solutions and fluxes, and high order Gauss quadratures for the inner products, e.g. [20]. Geometric flexibility comes from the multi-element subdivision of the domain. The DGSEM is now developed to the point of being efficient for large scale engineering level computations, e.g. [1],[3],[8], among others.

Robustness, however, has been an issue with the DGSEM at high order. It usually works, but it can go unstable even when the solutions are smooth. For nonlinear problems, this is probably not surprising. Examples are demonstrated in the computation of the Taylor-Green vortex problem, where instability at high orders is seen [11]. But instability arises even in linear problems when the coefficients are variable, which can come from inherent variability [2] or from variability introduced by curved elements [15]. The instability, we will show, comes from aliasing errors associated with the products of polynomials and insufficient Gauss quadrature precision.

Robust (provably stable) versions of the DGSEM that start from a split form of the partial differential equation (PDE) have recently been developed for linear hyperbolic systems for static [15] and moving domains [16]. In addition to stability, the approximations match the additional conservative and constant state preserving properties of the PDE [19]. The approach is applicable to nonlinear problems, where, depending on the equations and split form, the methods are energy or entropy stable [9],[13],[12].

In this paper, we introduce a polynomial spectral calculus that allows us to mirror the continuous PDE analysis to show stability of Black’s and the split-form approximations with a simple, compact notation applicable to any number of space dimensions. For the split form method, we also show how to use the calculus to demonstrate conservation and constant state preservation. The key starting point of the calculus is the summation by parts property satisfied by the Gauss-Lobatto quadrature [14], which allows us to write discrete versions of the Gauss law and its variants. Those discrete Gauss laws, in turn, allow us to write algebraically equivalent forms of the approximations, with which we can easily analyze their properties.

2 Linear Hyperbolic Problems on Bounded Domains

As examples of the use of the discrete calculus, we will analyze two discontinuous Galerkin spectral element approximations to the linear system of conservation laws

[TABLE]

where $\mathbf{u}\left(\vec{x},t\right)=[u_{1}\;u_{2}\;\ldots\;u_{p}]^{T}$ is the state vector and

[TABLE]

is the linear flux space-state vector. For simplicity we will assume that the system has been symmetrized and is hyperbolic so that

[TABLE]

for any $\left\|\vec{\alpha}\right\|^{2}_{2}=\sum\limits_{m=1}^{3}{\alpha_{m}^{2}}\neq 0$ and some real diagonal matrix $\Lambda$. We will also assume that the matrices $\mathcal{A}^{(m)}$ have bounded derivatives in the sense that

[TABLE]

where $\left\|\cdot\right\|_{2}$ is the matrix 2-norm. Additional constraints on the coefficient matrices need to be added later to ensure that the derivatives of their interpolants converge in the maximum norm. The product rule applied to (1 leads to the nonconservative form of the system

[TABLE]

With appropriate initial and characteristic boundary conditions on a bounded domain $\Omega\in\mathbb{R}^{3}$ the problem is (i) well posed, (ii) conservative, and, under conditions on $\vec{\mathcal{A}}$, (iii) preserves a constant state. These properties are demonstrated from a weak form of the average of the conservative, (1), and nonconservative, (5), forms of the equation, the so-called “split-form”. To write the weak form, we define the $\mathbb{L}^{2}$ inner product and norm

[TABLE]

Then for any state vector $\boldsymbol{\phi}\in\mathbb{L}^{2}\left(\Omega\right)$,

[TABLE]

From vector calculus, we have the extended Gauss law,

[TABLE]

where $\hat{n}$ is the outward unit normal. We write (8) in inner product form as

[TABLE]

We can apply the extended Gauss law to the inner products in the braces in (7) and use the fact that $\vec{\mathcal{A}}$ is symmetric to get an equivalent form that separates the boundary and volume contributions

[TABLE]

Constant state preservation, conservation and well-posedness are shown with judicious choices of $\mathbf{u}$ and $\boldsymbol{\phi}$. To find under what conditions a constant state is preserved, set $\mathbf{u}=\mathbf{c}=\text{constant}$ in (7) to see that

[TABLE]

from which it follows that $\mathbf{u}_{t}=0$ if ${\nabla\cdot\vec{\mathcal{A}}}=0$.

Global conservation is shown by selectively choosing each component of the state vector $\boldsymbol{\phi}$ in (10) to be unity and again noting that the coefficient matrices are symmetric to see that the terms in the braces cancel to leave

[TABLE]

To find conditions under which the initial boundary value problem is well-posed, we choose $\boldsymbol{\phi}=\mathbf{u}$ in (7) and note that

[TABLE]

Replacing those terms in (7) and multiplying the equation by two gives

[TABLE]

Gauss’ theorem allows us to replace the second term by a surface integral so

[TABLE]

We bound the right hand side by

[TABLE]

so

[TABLE]

Integrating over the time interval $[0,T]$ we write the energy in terms of the initial value and a boundary integral

[TABLE]

To properly pose the problem we must impose appropriate boundary conditions. From (3), we separate the waves traveling to the left and right of the boundary relative to $\hat{n}$ as

[TABLE]

where $\Lambda^{\pm}=\Lambda\pm\left|\Lambda\right|$ and we have left off the explicit dependence on $\hat{n}$. When we replace the values of $\mathbf{u}$ along the boundary associated with the incoming $\Lambda^{-}$ waves with a boundary state, $\mathbf{g}$, the solution can be bounded in terms of the initial and boundary data,

[TABLE]

Furthermore, if $\nabla\cdot\vec{\mathcal{A}}=0$, $\gamma=0$ and the energy does not grow in time except for energy introduced at the boundaries,

[TABLE]

3 A Polynomial Spectral Calculus

To follow the continuous problem analysis as closely as possible, we introduce a discrete calculus that looks and behaves like the continuous one as much as possible. We define the calculus for the reference domain $E=[-1,1]^{3}$ with coordinates $\vec{\xi}=\left(\xi,\eta,\zeta\right)=\xi\hat{\xi}+\eta\hat{\eta}+\zeta\hat{\zeta}=\sum\limits_{m=1}^{3}{{\xi^{(m)}}{{\hat{\xi}}^{m}}}$. Corresponding forms hold for two dimensional problems.

We represent functions of the reference domain coordinates by polynomials of degree $N$ or less, i.e. as elements of $\mathbb{P}^{N}(E)\subset\mathbb{L}^{2}(E)$. A basis for the polynomials on $E$ is the tensor product of the one dimensional Lagrange basis. Using that basis, we write a polynomial, $U$, in terms of nodal values ${U}_{ijk}={U}\left(\xi_{i},\eta_{j},\zeta_{k}\right)$ as an upper case letter, which for three space dimensions is

[TABLE]

where

[TABLE]

is the one-dimensional Lagrange interpolating polynomial with the property $\ell_{l}\left(s_{m}\right)=\delta_{lm}$, $l,m=0,1,2,\ldots,N$. The points $s_{i},\;i=0,1,2,\ldots,N$ are the interpolation points, whose locations are chosen below. We also write the interpolation operator, $\mathbb{I}^{N}:\mathbb{L}^{2}\rightarrow\mathbb{P}^{N}$, which projects square integrable functions on $E$ onto polynomials, as

[TABLE]

The use of the tensor product means that one and two dimensions are special cases of three dimensions, which is why we concentrate on three dimensional geometries here.

Derivatives of polynomials on $E$ evaluated at the nodes can be represented by matrix-vector multiplication. For instance,

[TABLE]

where $\mathcal{D}$ is the derivative matrix. The gradient and divergence of a polynomial in three space dimensions evaluated at a point $\left(\xi_{n},\eta_{m},\zeta_{l}\right)$ are therefore

[TABLE]

The use of the calculus that we develop depends on the choice that the interpolation nodes, $s_{i}$, are the nodes of the Legendre-Gauss-Lobatto (LGL) quadrature. We represent the one dimensional LGL quadrature of a function $g(s)$ as

[TABLE]

where the $\omega_{i}$ are the LGL quadrature weights. The quadrature is exact if $g\in\mathbb{P}^{2N-1}$. By tensor product extension, we write three dimensional volume integral approximations as

[TABLE]

where $\omega_{ijk}=\omega_{i}\omega_{j}\omega_{k}$. Two-dimensional surface integral approximations are

[TABLE]

Two space dimensional areas and edge integrals are defined similarly.

We define the discrete inner product of two functions $f$ and $g$ and the discrete norm of $f$ from the quadrature

[TABLE]

The definition is extended for vector arguments like

[TABLE]

for a state vector $\mathbf{f}^{(m)}=[f^{(m)}_{1}\;f^{(m)}_{2}\;\ldots\;f^{(m)}_{p}]^{T}$ as

[TABLE]

and similarly for other arguments.

The Lagrange basis functions are orthogonal with respect to the discrete inner product [6]. In one space dimension, for instance, $\left(\ell_{i},\ell_{j}\right)_{E,N}=\omega_{j}\delta_{ij}$. Also, from the definitions of the interpolation operator and the discrete inner product,

[TABLE]

Finally, the discrete norm is equivalent to the continuous norm [7] in that for $U\in\mathbb{P}^{N}$,

[TABLE]

where $C$ is a constant.

The crucial property for the analysis of the discrete approximation is the summation by parts (SBP) property satisfied by the LGL quadrature. Let $U,V\in\mathbb{P}^{N}$. Then exactness of the LGL quadrature implies that

[TABLE]

The result extends to all space dimensions [14] with

[TABLE]

We can use (35) and (36) to formulate a discrete integral calculus. If we replace $U$ in (36) by the components of a vector $\vec{F}$, and sum, we get the Discrete Extended Gauss Law (DXGL) originally derived in [14]: For any vector of polynomials $\vec{F}\in\mathbb{P}^{N}$ and any polynomial $V\in\mathbb{P}^{N}$,

[TABLE]

where $\hat{n}$ is the unit outward normal at the faces of $E$. Carrying this further, if we set $V=1$ we get the Discrete Gauss Law (DGL)

[TABLE]

The DGL is exact for polynomial arguments. By using the appropriate definitions for the inner products, both discrete Gauss laws extend to hold for state vectors $\vec{\mathbf{F}}$ and $\mathbf{V}$.

Next, we see that if we replace the vector flux $\vec{F}$ in (37) with $\nabla\Phi\in\mathbb{P}^{N}$, then we get the discrete version of Green’s first identity,

[TABLE]

Swapping the variables $\Phi$ and $V$ and subtracting from the original gives Green’s second identity

[TABLE]

The discrete Green’s identities would be useful to prove stability of continuous Galerkin spectral element methods of second order problems.

Other identities that do not involve quadratic products of polynomial arguments hold discretely through exactness of the LGL quadrature. For instance,

[TABLE]

and

[TABLE]

What we see, then, is that the well-known integral identities hold due to either integration or summation by parts.

Whereas integration rules hold discretely, product differentiation rules do not usually hold because differentiation and interpolation do not always commute. For instance, the product rule does not generally hold. That is, for polynomials $U,V$,

[TABLE]

unless the product $UV\in\mathbb{P}^{N}$. Differentiation and interpolation do not commute because of aliasing errors that arise from projecting the product onto a polynomial of degree less than or equal to $N$ [6].

4 Discontinuous Galerkin Spectral Element Approximations

We now use the polynomial calculus introduced in (3) to formulate and analyze discontinuous Galerkin spectral element approximations in three space dimensions. The steps to derive two dimensional approximations are identical. The domain $\Omega$ is subdivided into $N_{el}$ nonoverlapping hexahedral elements, $e^{r},r=1,2,\ldots,N_{el}$. We assume here that the subdivision is conforming. Each element is mapped from the reference element $E$ by a transformation $\vec{x}=\vec{X}\left(\vec{\xi}\right)$. From the transformation, we define the three covariant basis vectors

[TABLE]

and (volume weighted) contravariant vectors, formally written as

[TABLE]

where

[TABLE]

is the Jacobian of the transformation.

Under the mapping, the divergence of a spatial vector flux can be written compactly in terms of the reference space variables as

[TABLE]

The vector $\tilde{{\mathbf{f}}}$ is the volume weighted contravariant flux whose components are $\tilde{{\mathbf{f}}}^{i}={\mathcal{J}{{\vec{a}}^{i}}\cdot\vec{{\mathbf{f}}}}$.

The conservation law is then represented on the reference domain by another conservation law

[TABLE]

where we have defined the (volume weighted) contravariant coefficient matrices

[TABLE]

and

[TABLE]

We can also construct the nonconservative form of the system on the reference domain using the chain rule,

[TABLE]

We construct weak forms of the conservative and nonconservative equations by taking the inner product of the equations with a test function $\boldsymbol{\phi}\in\mathbb{L}^{2}(E)$ and applying extended Gauss Law to the space derivative terms,

[TABLE]

and

[TABLE]

When we average the two equations (52) and (53) we get the split weak form

[TABLE]

4.1 The DGSEM

The original DG spectral element method introduced by Black [4] starts from the conservative weak form (52). We use the calculus now to show that it is stable if the coefficient matrices $\tilde{\mathcal{A}}$ are constant. If, in addition, characteristic boundary conditions are used at physical boundaries, the approximation is optimally stable in the sense that the global energy discretely matches (21).

To construct the approximation, one approximates the solutions, fluxes, coefficient matrices and Jacobian with polynomial interpolants on element $e^{r}\rightarrow E$ by

[TABLE]

From this point, we leave off the superscripts $r$ and subscripts $\xi$ on $\nabla_{\xi}$ unless necessary.

To continue the contruction, one replaces the continuous inner products by the discrete inner products, here being Gauss-Lobatto quadratures. The normal boundary flux is replaced by a consistent numerical flux, $\tilde{\mathbf{f}}\leftarrow\tilde{\mathbf{F}}^{*}\left(\mathbf{U}^{L},\mathbf{U}^{R};\hat{n}\right)$ where $\mathbf{U}^{L,R}$ are the left and right states at the element boundary, measured with respect to the outward normal, $\hat{n}$. The numerical flux ensures continuity of the normal flux at element faces. Finally, $\boldsymbol{\phi}$ is restricted to elements of $\mathbb{P}^{N}$. The result of the approximations is the formal statement of the method

[TABLE]

Details for going from the formal statement to the form to implement can be found in [17].

Alternate, yet algebraically equivalent forms of the DGSEM can be derived by applying the DXGL. For instance, if we apply the DXGL to the last inner product in (56) we get the algebraically equivalent form

[TABLE]

If the contravariant coefficient matrices are constant, implying that the original problem is constant coefficient and the elements are rectangular in shape, then the DGSEM approximation is strongly stable. To show stability, we set $\boldsymbol{\phi}=\mathbf{U}$ in (57) and define the volume weighted norm

[TABLE]

Then

[TABLE]

With constant coefficients, the volume term in (59) can be converted to a surface quadrature. The coefficient matrices being constant and symmetric implies that

[TABLE]

The key step is the second marked with the “*”, where the product rule applies because $\tilde{\mathcal{A}}\mathbf{U}\in\mathbb{P}^{N}$ when $\tilde{\mathcal{A}}$ is constant. We then substitute the equivalence into the DXGL to see that

[TABLE]

Therefore, the local energy changes according to

[TABLE]

and stability depends solely on what happens on the element faces.

The change in the total energy is found by summing over all the elements. Although the numerical flux is continuous at element interfaces, the solution and flux are discontinuous. If we define the jump in a quantity with the usual notation \left\llbracket V\right\rrbracket=V^{R}-V^{L}, then

[TABLE]

Stability is determined, therefore, only by the influence of the jumps at the element boundaries and the physical boundary approximations through the numerical flux. For linear problems, it is natural to choose an upwinded or central flux,

[TABLE]

where $\sigma=0$ is the central flux and $\sigma=1$ is the fully upwind flux. With this flux [18],

[TABLE]

so that the interior face terms in (63) are dissipative. To match the PDE energy bound, (21), the fully upwind flux must be used at the physical boundaries. With exterior values $\mathbf{g}$ set along incoming characteristics [21] and when $\sigma=1$,

[TABLE]

If we define the total energy by

[TABLE]

and integrate (63)in time, the total energy satisfies (c.f. (20))

[TABLE]

Finally, if the interpolant of the Jacobian is bounded from below, $J>0$, then for some positive constants $c$ and $C$ [18],

[TABLE]

which says that, like the continuous solution, the energy approximate solution is bounded by the data in the continuous norm over the entire domain.

4.2 Stabilization by Split Form

If the contravariant coefficient matrices are not constant, then the key step in (60) does not hold because interpolation and differentiation do not commute. We show now that stability hangs on whether or not the dissipation introduced by the numerical flux at the element interfaces and by the characteristic boundary conditions is sufficient to counterbalance the aliasing errors associated with the volume term that remains. That balance shows why the approximation [DGSEM] can be, but does not have to be, stable for variable coefficient problems or curved elements.

The use of the polynomial calculus allows us to quickly and compactly construct four algebraically equivalent representations of a split form approximation [18] that is strongly stable, constant state preserving and globally conservative for non-constant coefficient problems where the coefficient variation is due to inherent variability in the PDE and/or due to variability introduced by the coefficient mappings from curved elements to the reference element. It also allows us to simplify the analysis done, for example, in [19].

The result of applying the approximations (55) and LGL quadrature to (54) is the first split form of the DGSEM used in [18]. In accordance to common terminology, this is the “weak” form

[TABLE]

We get alternative, yet algebraically equivalent forms by applying the DXGL (37) to selected terms in (70). When we apply the DXGL to the first two inner products in the braces and use the fact that the coefficient matrices are symmetric we get the “strong” form

[TABLE]

If we rearrange the terms in $[S]$ to “strong+correction” form

[TABLE]

we see that the split form approximation is the strong form of the original DGSEM (57) plus a correction term in the amount by which the product rule (43) does not hold. When the product rule does hold, such as when the contravariant coefficient matrices are constant, the correction term vanishes and we are back to the original scheme of Black, [DGSEM].

We get a fourth algebraically equivalent “directly stable” form by applying the DXGL to only the first inner product in the braces of the weak form $[W]$,

[TABLE]

Any of the four equivalent forms $[W]\Leftrightarrow[S]\Leftrightarrow[SC]\Leftrightarrow[DS]$ can be used as is convenient for computation or theory. For instance, to show conservation, choose the form $[W]$ and selectively set each component of $\boldsymbol{\phi}$ to one. Then the first inner product in the braces vanishes and the second and third cancel leaving

[TABLE]

Summing over all elements, the interior face contributions cancel leaving the global conservation statement

[TABLE]

To find conditions under which the approximation is constant state preserving, use the form $[S]$ with $\mathbf{U}=\mathbf{c}=\text{const}$ in all elements. The first and third inner products in the braces vanish provided that $\nabla\cdot\left(\mathbb{I}^{N}\left(\tilde{\mathcal{A}}\right)\right)=0$, and the second is explicitly zero. Consistency of the numerical flux implies that $\tilde{\mathbf{F}}^{*}\left(\mathbf{c},\mathbf{c};\hat{n}\right)=\tilde{\mathbf{F}}\cdot\hat{n}$. Therefore, $\left(J\mathbf{U}_{t},\boldsymbol{\phi}\right)_{N}=0$ for all $\boldsymbol{\phi}\in\mathbb{P}^{N}$, which implies that at each node $nml$ in each element $r$, $d\mathbf{U}^{r}_{nml}/dt=0$.

Finally, the split form approximation is optimally stable in the sense that with the numerical flux (64), the norm of the approximate solution satisfies an energy statement like (20). We show stability using $[DS]$ and $\boldsymbol{\phi}=\mathbf{U}$. With the substitution, the volume terms represented by the first two inner products in the braces immediately cancel. The third inner product in the braces can be bounded

[TABLE]

and under assumptions on the smoothness of $\vec{\mathcal{A}}$ [23] and positivity of the Jacobian [18] the coefficient $\hat{\gamma}$ will converge spectrally to $\gamma$. (If the divergence of the interpolant vanishes, then $\hat{\gamma}=0$.) With the bound on the divergence of the coefficient matrices,

[TABLE]

The change in the total energy is again found by summing over all the elements. If we introduce the integrating factor $\hat{\gamma}=\mathop{\max}\limits_{r}\hat{\gamma}^{r}$,

[TABLE]

The interface and boundary terms on the right hand side of (78) are identical to what appeared in the original DGSEM, (63). Therefore, the total energy satisfies

[TABLE]

As with the continuous solution, (21), if $\nabla\cdot\mathbb{I}^{N}\left(\tilde{\mathcal{A}}\right)=0$, $\hat{\gamma}=0$, and

[TABLE]

Applying the norm equivalence (69), we see that the split form approximation is strongly stable for variable coefficient problems and/or curved elements.

Finally, we can use the stability analysis of the split form approximation to write the conditions needed for the original DGSEM to be stable when the coefficients are variable. For simplicity, let us suppose that $\nabla\cdot\vec{\mathcal{A}}=0$ and external boundary states $\mathbf{g}=0$ so that the global energy should not increase and instability is not masked by natural growth. Let us also assume that $\nabla\cdot\mathbb{I}^{N}\left(\tilde{\mathcal{A}}\right)=0$ so that $\hat{\gamma}$ also vanishes. Then by (79), the energy of the split form approximation does not grow. With $[DS]\Leftrightarrow[SC]$ and $[DGSEM]\Leftrightarrow[SC]-[C]$, where $[C]$ is the correction term

[TABLE]

the elemental energy for the DGSEM satisfies

[TABLE]

The term in the braces of the volume term is non-zero unless the product rule holds. Therefore, for the DGSEM to be stable when the coefficients are variable, the surface terms (including the dissipation arising from the physical boundaries seen in (66) must be sufficiently large to counteract any destabilizing influence of the volume term, which might require trying more dissipative numerical fluxes than the characteristic upwind flux. Practice has shown that at least at low order one can often find numerical fluxes for which the influence of the surface terms is sufficiently dissipative. But (82) shows that the approximation can be unstable if the aliasing growth contribution is larger than the dissipation contribution from the element faces.

5 Summary

In this paper, we described a discrete integral spectral calculus for polynomial spectral methods using Legendre-Gauss-Lobatto quadrature. This calculus allowed us to write and analyze discontinuous Galerkin spectral element approximations in a compact notation consistent with the continuous version. In particular, it is possible to easily derive four algebraically equivalent forms of a split form approximation for linear hyperbolic systems. These four equivalent forms can then be used to show global conservation, constant state preservation (when applicable) and, most importantly, strong stability of the split form approximation for variable coefficient problems on curved elements.