A Provably Stable Discontinuous Galerkin Spectral Element Approximation for Moving Hexahedral Meshes

# A Provably Stable Discontinuous Galerkin Spectral Element Approximation for Moving Hexahedral Meshes

## Abstract

We design a novel provably stable discontinuous Galerkin spectral element (DGSEM) approximation to solve systems of conservation laws on moving domains. To incorporate the motion of the domain, we use an arbitrary Lagrangian-Eulerian formulation to map the governing equations to a fixed reference domain. The approximation is made stable by a discretization of a skew-symmetric formulation of the problem. We prove that the discrete approximation is stable, conservative and, for constant coefficient problems, maintains the free-stream preservation property. We also provide details on how to add the new skew-symmetric ALE approximation to an existing discontinuous Galerkin spectral element code. Lastly, we provide numerical support of the theoretical results.

###### keywords:
discontinuous Galerkin spectral element method, summation-by-parts, moving mesh, arbitrary Lagrangian-Eulerian, energy stable, free-stream preservation

## 1 Introduction

Many applications in computational physics require the numerical approximation of a system of hyperbolic conservation laws on moving domains, for example problems in fluid dynamics etienne2009perspective (); Lomtevetal1999 (); Mavriplis:2006mb (); wang2015high () or electromagnetism censor2004 (); harfoush1989 (); ho2006 (). A common approach to approximate solutions of partial differential equations (PDEs) with moving boundaries is to use an arbitrary Lagrangian-Eulerian (ALE) formulation etienne2009perspective (); Hay2014204 (). The ALE method maps a time dependent domain, , enclosed by the moving boundaries onto a fixed reference domain, . Conveniently, systems of conservation laws in the original moving domain are transformed to a set of conservation law equations in the reference domain. In the numerical approximation on the reference domain, the new set of equations depends on the mesh velocity.

A discontinuous Galerkin spectral element method for moving domains (DGSEM-ALE) that was spectrally accurate in space, free-stream preserving, and retained the full time accuracy of the time integrator was introduced in Minoli:2010rt (). Extensions of the method were presented in ISI:000302933600007 (), and Winters:2013nx (); Winters:2013oq (); Winters:2014hl (), the latter of which addressed the approximation of problems with discontinuous and moving material interfaces and efficiency through the addition of local time stepping.

None of the papers on the DGSEM-ALE approximation directly addressed stability, however. In fact, it has been noted that even on static domains, discontinuous Galerkin spectral element approximations for variable coefficient problems Beck:2013sf (), and problems that become variable coefficient because of curved elements kopriva2015 (), can be unstable.

In this paper, we now address the additional issues of robustness and the ability of a moving mesh method to be guaranteed stable. Recent work for static meshes has focused on the use of DGSEM approximations written in a skew-symmetric form to guarantee stability, e.g. gassner2015 (); Gassner:2013uq (); kopriva2015 (). The skew-symmetric form of a problem is usually found by averaging the conservative form and a non-conservative advective form of the equation. This is problematic, because it is not obvious that discretizations of the skew-symmetric form remain conservative. Recent success has been had using diagonal norm summation-by-parts (SBP) operators to discretize the spatial derivatives in the skew-symmetric formulation skew_sbp2 (); gassner_skew_burgers (); gassner_kepdg (); gassner2015 (). There is a known link between SBP methods and the discontinuous Galerkin spectral element approximation with Legendre-Gauss-Lobatto points, e.g. gassner_skew_burgers (). The approximation developed here will also use a skew-symmetric formulation.

In addition to stability, the DGSEM approximation that we propose is conservative, high order accurate in both space and time, and ensures that for constant coefficient problems a constant solution of the conservation law remains constant, discretely, for all time, i.e, the approximation possesses free-stream preservation. Failure to satisfy free-stream preservation usually implies that the motion of the mesh can create spurious waves and may introduce discrete errors that can lead to wave misidentification, even for spectrally accurate approximations Kopriva:2006er ().

A stable deforming mesh approximation for high-order finite difference schemes with the summation by parts (SBP) property was recently proposed by Nikkar and Nordström Nikkar:2014si (). The method development here parallels theirs, and the result satisfies the same type of energy estimate, even though our use of the weak rather than strong formulation, the notation, and the approximation differ.

The remainder of this paper is organized as follows: Sec. 2 reviews the ALE mapping between a reference space and a general curvilinear coordinate system. Next, in Sec. 3, we demonstrate the well-posedness of the ALE formulation, which presents the target to be approximated. In Sec. 4 we use a skew-symmetric formulation of the governing equations to develop a stable approximation of the problem on moving meshes. Sec. 5 provides proofs of the conservation, stability, and free-stream preserving properties of the skew-symmetric DGSEM-ALE formulation. We provide some details on the implementation of the newly proposed scheme in Sec. 6. Numerical results are presented in Sec. 7 to support the theoretical findings. Finally, concluding remarks are given in Sec. 8.

## 2 Arbitrary Lagrangian-Eulerian (ALE) Formulation of Conservation Laws in a Curvilinear Coordinate System

We will derive a discontinuous Galerkin spectral element method for systems of partial differential equations of the form

 qt+∇⋅→f=qt+3∑i=1∂fi∂xi=0, (2.1)

on a three dimensional domain with moving boundaries, , where . Here we denote the solution and flux vector components by bold face and spatial vectors by arrows. We assume that the system is symmetric hyperbolic, with covariant flux components

 fi=Ai(→x)q, (2.2)

where the are matrices. If the system is not symmetric to start with, then symmetrization, which is available for most systems of interest mccarthy1980 (); R.-F.-Warming:1975fk (), is applied first. We further assume that the matrices are smooth, having bounded derivatives. Since the system is hyperbolic, there exists an invertible matrix such that

 A=3∑i=1αiAi=P(A)Λ(A)P−1(A), (2.3)

for any with and some real diagonal matrix .

As a concrete example, the symmetric wave equation can be written in the form (2.1) as

 ⎡⎢ ⎢ ⎢⎣puvw⎤⎥ ⎥ ⎥⎦t+⎛⎜ ⎜ ⎜⎝⎡⎢ ⎢ ⎢⎣0c00c00000000000⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣puvw⎤⎥ ⎥ ⎥⎦⎞⎟ ⎟ ⎟⎠x+⎛⎜ ⎜ ⎜⎝⎡⎢ ⎢ ⎢⎣00c00000c0000000⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣puvw⎤⎥ ⎥ ⎥⎦⎞⎟ ⎟ ⎟⎠y+⎛⎜ ⎜ ⎜⎝⎡⎢ ⎢ ⎢⎣000c00000000c000⎤⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣puvw⎤⎥ ⎥ ⎥⎦⎞⎟ ⎟ ⎟⎠z=0. (2.4)

In the ALE formulation, one maps onto a reference domain, , by a transformation

 (→x,t)=→X(→ξ,τ), (2.5)

where is a three dimensional curvilinear coordinate on the reference domain. For convenience with the approximations later, we can take the reference domain to be the reference cube . The mapping has a set of covariant basis vectors, defined by

 →ai=∂→X∂ξii=1,2,3. (2.6)

From the covariant basis, one can formally define the contravariant basis , multiplied by the Jacobian of the transformation, , by

 J→ai=J∇ξi=→aj×→ak=∂→X∂ξj×∂→X∂ξk,(i,j,k)cyclic. (2.7)

The Jacobian itself can be written in terms of the covariant vectors as

 J=→ai⋅(→aj×→ak),(i,j,k)cyclic. (2.8)

The geometry satisfies the well-known metric identities Vinokur1974 (), (Thompsonetal1985, , Chpt. III)

 3∑i=1∂J→ai∂ξi=0, (2.9)

and the Geometric Conservation Law (GCL) Mavriplis:2006mb ()

 ∂J∂τ−3∑i=1∂J→ai⋅→xτ∂ξi=0. (2.10)

Under the transformation (2.5), see, e.g. Minoli:2010rt (), the conservation law (2.1) remains a conservation law

 ∂Jq∂τ+3∑i=1∂J→ai⋅(→F−q→xτ)∂ξi=0, (2.11)

on the reference domain. If we define the contravariant coefficient matrices,

 ~Ai=J→ai⋅3∑j=1Aj^xj,i=1,2,3, (2.12)

where denotes the unit vector in the coordinate direction, then we can re-write (2.11) in fully conservative form as

 ∂Jq∂τ+3∑i=1∂~f∂ξi=0, (2.13)

where and is the identity matrix.

Using the metric identities, the GCL and the conservative form of the equations, we can also derive a non-conservative form of the equations. From the chain rule,

 Jτq+Jqτ+3∑i=1∂~Ai∂ξiq+3∑i=1~Ai∂q∂ξi−3∑i=1∂J→ai⋅→xτ∂ξiq−3∑i=1J→ai⋅→xτ∂q∂ξi=0. (2.14)

Applying the GCL (2.10) we find,

 Jqτ+3∑i=1~Ai∂q∂ξi−3∑i=1J→ai⋅→xτ∂q∂ξi+3∑i=1∂~Ai∂ξiq=0, (2.15)

giving us the nonconservative system of equations

 Jqτ+3∑i=1(~Ai−J→ai⋅→xτI)∂q∂ξi+3∑i=1∂~Ai∂ξiq=0. (2.16)

Finally, let us define the matrices

 Ai(→ξ,τ)=~Ai−J→ai⋅→xτI,i=1,2,3, (2.17)

to let us write the conservative and non-conservative forms of the equations as

 ∂Jq∂τ+3∑i=1∂∂ξi(Aiq)=0, (2.18)

and

 Jqτ+3∑i=1Ai∂q∂ξi+3∑i=1∂~Ai∂ξiq=0, (2.19)

respectively.

## 3 Well-Posedness of the ALE Formulation

Since the stability proof mimics that of well-posedness of the PDE, we first show that the system of PDEs on the moving domain is well-posed when appropriate boundary conditions are applied. The derivation here follows that of Nikkar:2014si (), but uses notation that is consistent with how we will write the DGSEM and its stability proof.

Let us define the inner product on the reference domain as

 (u,v)=∫EuTvd→ξ, (3.1)

and norm . We denote the space .

We show well-posedness by bounding the time rate of change of the energy in the moving domain , which is equivalent to

 ddτ∥q∥2J≡ddτ∫EqTqJd→ξ=(q,∂Jq∂τ)+(q,Jqτ). (3.2)

The two terms on the right of (3.2) can be found from the conservative and non-conservative forms of the mapped system. The inner product of the conservative form of the equation (2.18) with the solution is

 (q,∂Jq∂τ)+3∑i=1⎛⎜ ⎜⎝q,∂(Aiq)∂ξi⎞⎟ ⎟⎠=0. (3.3)

Similarly, using the nonconservative form (2.19),

 (3.4)

Adding the conservative (3.3) and non-conservative (3.4) forms together gives

 (q,∂Jq∂τ)+(q,Jqτ)+3∑i=1⎛⎜ ⎜⎝q,∂(Aiq)∂ξi+Ai∂q∂ξi⎞⎟ ⎟⎠+3∑i=1(q,∂~Ai∂ξiq)=0. (3.5)

Also, because the matrices are symmetric, we see that

 ∂∂ξi(qTAiq)=∂∂ξi(qT)Aiq+qT∂(Aiq)∂ξi=qTAi∂∂ξi(q)+qT∂(Aiq)∂ξi, (3.6)

so

 3∑i=1(q,∂∂ξi(Aiq)+Ai∂q∂ξi)=∫E3∑i=1∂∂ξi(qTAiq)d→ξ. (3.7)

Then with (3.5) and (3.6), (3.2) can be written with the divergence of a flux as

 (3.8)

Now define

 →A=3∑i=1Ai^ξi, (3.9)

where are the cardinal bases and the outward normal on the reference box. Then Gauss’ theorem states that

 ddτ∥q∥2J+∫∂EqT→A⋅^nqdS=−3∑i=1(q,∂~Ai∂ξiq). (3.10)

Next, we determine the bound

 −3∑i=1(q,∂~Ai∂ξiq)=(q,{−1J3∑i=1∂~Ai∂ξi}Jq)⩽maxE∣∣ ∣∣1J3∑i=1∂~Ai∂ξi∣∣ ∣∣∥q∥2J. (3.11)

We then note that under the transformation rules,

 1J3∑i=1∂~Ai∂ξi=3∑i=1∂Ai∂xi, (3.12)

is assumed to be bounded. Therefore,

 ddτ∥q∥2J+∫∂EqT→A⋅^nqdS⩽2γ∥q∥2J, (3.13)

where

 γ=12maxE∣∣ ∣∣1J3∑i=1∂~Ai∂ξi∣∣ ∣∣<∞. (3.14)

Integrating both sides of (3.13) with respect to time over a time interval gives

 ∥q(T)∥2J⩽e2γT{∥q(0)∥2J−∫T0∫∂Ee−2γτqT→A⋅^nqdSdτ}. (3.15)

We now apply boundary conditions to (3.15). First we note that is diagonalizable. From the definition,

 Ai=J→ai⋅(3∑j=1Aj^xj−→xτI), (3.16)

so the normal flux matrices are given by

 3∑i=1Aini=3∑j=1Aj{3∑i=1J→ai⋅^xj^ni}−{3∑i=1J→ai⋅^xτ^ni}I=3∑j=1Ajαj−βI. (3.17)

Since is diagonalizable, so is . This means that we can split the boundary contributions into incoming and outgoing components according to the sign of the eigenvalues. So let us split the normal coefficient matrix

 3∑i=1Aini=PΛ+P−1+PΛ−P−1=A++A−, (3.18)

so that

 ∥q(T)∥2J≤e2γT{∥q(0)∥2J−∫T0∫∂Ee−2γτqTPΛ+P−1qdSdτ−∫T0∫∂Ee−2γτqTPΛ−P−1qdSdτ}. (3.19)

We then replace the incoming values (associated with the eigenvalues) with the exterior state and bound the integrals over time to get

 (3.20)

The statement (3.20) says that the initial boundary value problem is strongly well posed Lorenz:1989fk (). When we set the boundary states to zero, we see that the system of equations is well-posed,

 ∥q(T)∥J⩽eγT∥q(0)∥J. (3.21)

Finally, if the matrices are also constant, then and the energy never grows,

 ∥q(T)∥J⩽∥q(0)∥J. (3.22)

## 4 A Stable DGSEM-ALE for Moving Domains

We now derive a discontinuous Galerkin spectral element method (DGSEM) for moving elements whose stability properties mimic (3.20). A description of the standard approximation can be found in Kopriva:2009nx (). We subdivide the domain into non-overlapping, geometrically conforming hexahedral elements that cover . Since has moving boundaries, so too will the elements. We then map each element individually with a local time dependent mapping of the form (2.5) onto the reference element . Then on each element, the equations take on the conservative form (2.18) and the non-conservative form

 ∂J∂τq+J∂q∂τ+3∑i=1∂Ai∂ξiq+3∑i=1Ai∂q∂ξi=0, (4.1)

which is written without applying the metric identities and GCL.

To create the skew-symmetric form, we average the conservative and nonconservative forms Gassner:2013uq () to get

 12{∂Jq∂τ+∂J∂τq+J∂q∂τ}+123∑i=1∂(Aiq)∂ξi+12{3∑i=1∂Ai∂ξiq+3∑i=1Ai∂q∂ξi}=0. (4.2)

We construct a weak form of (4.2) by multiplying by a test function and integrating over the domain. In inner product notation, the weak form is

 (4.3)

Next we integrate terms that have derivatives of the solution by parts. Note that the coefficient matrices are symmetric so that we can write

 (4.4)

where we have introduced the shorthand notation

 (→A⋅^nq)Tϕ∣∣∣∂E≡6∑s=1{∫faces(→A⋅^nsq)TϕdSs}, (4.5)

for the boundary face contributions. Note that is the normal flux at the face . Eq. (4.4) is the weak form from which we will create our skew-symmetric approximation.

To get the approximation, (c.f. Kopriva:2009nx ()) we replace and by polynomial interpolants, the normal boundary and interface fluxes by the normal Riemann (numerical) flux, quadratic quantities by their polynomial interpolant and integrals by Legendre-Gauss-Lobatto quadrature.

We start by defining the polynomial interpolation operator. Let be the space of polynomials of degree less than or equal to . For some function defined on the reference element, the interpolant of through the tensor product of the Legendre-Gauss-Lobatto nodes is

 INv(→ξ)=N∑j,k,l=0vjklℓj(ξ)ℓk(η)ℓl(ζ)∈PN, (4.6)

where the is the Lagrange interpolating polynomial with nodes at the Legendre-Gauss-Lobatto points and is the value of at the tensor product of those points. We then approximate

 (4.7)

We also define the discrete inner product of two polynomials, as

 (U,V)N=N∑j,k,l=0UTjklVjklWjWkWl≡N∑j,k,l=0UTjklVjklWjkl, (4.8)

where the singly subscripted ’s are the one-dimensional Legendre-Gauss-Lobatto quadrature weights and the triply subscripted is the product of the three. We use a similar notation for integrals, where we add a subscript to denote quadrature.

We note two facts CHQZ:2006 () about the Legendre-Gauss-Lobatto quadrature that we will use later. The first is the exactness of the quadrature,

 (4.9)

The second is that for some function ,

 (IN(g),V)N=(g,V)N∀→V∈PN, (4.10)

which can be seen directly from the definition.

For the numerical flux (Riemann solver) we use the upwind () or central ()

 F∗(QL,QR;^n)=12{~F(QL)⋅^n+~F(QR)⋅^n}−λ∣∣IN(→A)⋅^n∣∣2{QR−QL}, (4.11)

which provides a unique flux given a left and right state relative to the (normal) vector .

We then make the substitutions of the approximations into the continuous weak form to get the discrete weak form

 12(∂IN(JQ)∂τ+∂J∂τQ+IN{J∂Q∂τ},ϕ)N+(F∗)Tϕ∣∣∂E,N−12{3∑i=1(~Fi,∂ϕ∂ξi)N}+12⎧⎪ ⎪⎨⎪ ⎪⎩3∑i=1⎛⎜ ⎜⎝∂IN(Ai)∂ξiQ,ϕ⎞⎟ ⎟⎠N−3∑i=1⎛⎜ ⎜⎝Q,∂IN(IN(Ai)ϕ)∂ξi⎞⎟ ⎟⎠N⎫⎪ ⎪⎬⎪ ⎪⎭=0. (4.12)

Finally, we separate out the parts that contribute to the GCL to get a formal statement of the approximate weak form

 (4.13)

Reduction to a Static Domain

Before moving on, we note that the approximation (4.13) is identical to the approximation in kopriva2015 () when the elements do not move. In the static case, time derivatives of the mesh and Jacobian vanish. Next, if we expand the quadratures,

 12(∂IN(JQ)∂τ+J∂Q∂τ,ϕ)N=12N∑j,k,l=0Wjklϕ(→ξjkl){∂JjklQjkl∂τ+Jjkl∂Qjkl∂τ}=(J∂Q∂τ,ϕ)N. (4.14)

Finally, since , . Then for a static domain,

 (4.15)

### 4.1 Approximation of the GCL

The Geometric Conservation Law (2.10) can be written as

 ˙J+∇ξ⋅→ψ=0, (4.16)

where

 →ψ=−3∑i=1J→ai⋅→xτ^ξi. (4.17)

For convenience, we approximate this with a DGSEM approximation simultaneously with the solution, so we write a weak form of the GCL as

 (Jτ,ϕ)+(∇ξ⋅→ψ,ϕ)=0. (4.18)

Integrating by parts (to put it into the same equation form as the solution),

 (Jτ,ϕ)+→ψ⋅^nϕ∣∣∂E−(→ψ,∇ξϕ)=0. (4.19)

We now approximate (4.19). This means that we approximate and . Note that by (2.7), the flux function, , is actually independent of the Jacobian and is dependent only on the current mesh and its velocity. As a result, (4.19) doesn’t describe a PDE for the Jacobian, but rather is an ODE. Following the recipe above, we replace inner products by Legendre-Gauss-Lobatto quadrature. Formally, we would also replace the boundary term by a Riemann solver. However, the normals (for a conforming mesh) and the mesh velocity are continuous at the faces, so we can simply use the computed values there. The approximation of the Jacobian is therefore

 (˙J,ϕ)N+→Ψ⋅^nϕ∣∣∂E,N−(→Ψ,∇ϕ)N=0. (4.20)

Furthermore, the discrete inner product satisfies the summation-by-parts (SBP) property gassner2010 (). For some and ,

 (∇⋅~F,ϕ)N=ϕ~F⋅^n∣∣∂E,N−(~F,∇ϕ)N, (4.21)

so with continuity at the boundaries for , the approximation (4.20) is algebraically equivalent to

 (˙J,ϕ)N+(∇⋅→Ψ,ϕ)N=0. (4.22)

Finally, we can combine the two inner products to get the equivalent statement

 (˙J+∇⋅→Ψ,ϕ)N=0. (4.23)

We now show that the last term in (4.13), which contains the GCL, vanishes if we compute the Jacobian using the DG approximation. To see this, it is important to note that the weak form is satisfied pointwise by using the quadrature. That is, we take the test function to be the tensor product basis, i.e. . Using the last form of the approximation, (4.23), we see that satisfies

 ˙Jjkl+(∇⋅→Ψ)jkl=0,j,k,l=0,1,…,N, (4.24)

or

 ˙Jjkl−⎛⎜ ⎜⎝3∑i=1∂IN(J→ai⋅→xτ)∂ξi⎞⎟ ⎟⎠jkl=0,j,k,l=0,1,…,N, (4.25)

when we expand the definition of . Eq. (4.25) holds at each Legendre-Gauss-Lobatto node. Therefore, multiplying by the solution vector, quadrature weight and test function at each node, and summing over all nodes,

 ⎛⎜ ⎜⎝⎡⎢ ⎢⎣∂J∂τ−3∑i=1∂IN(J→ai⋅→xτ)∂ξi⎤⎥ ⎥⎦Q,ϕ⎞⎟ ⎟⎠N=0. (4.26)

We will call (4.26) the weak discrete geometric conservation law or WDGCL. It is also equivalent to the other forms above.

###### Remark 0.

The approximations (4.20) (4.22) (4.23) (4.26). The equivalences will be important later when we want derivatives on the test function or not. They say we can use any of the discrete forms of the GCL as is convenient for theory or computation.

### 4.2 The Skew-Symmetric Approximation on Moving Meshes

Formally and compactly, provided that the discrete metric identities are satisfied Kopriva:2006er (), the skew-symmetric approximation on moving meshes for the Jacobian and solution is the geometric conservation law

 (Jτ,ϕ)N+→Ψ⋅^nϕ∣∣∂E,N−(→Ψ,∇ϕ)N=0, (4.27)

and solution approximation

 (4.28)

where

• .

• .

• .

## 5 Properties of the Skew-Symmetric Approximation

We now show that the skew-symmetric DGSEM-ALE approximations (4.27) and (4.28) are stable, conservative, and preserve a constant state when the ’s are constant.

### 5.1 Stability

The key feature of the skew-symmetric approximation is that it is stable, which follows if the WDGCL is satisfied as described in Sec. 4.1.

We first derive a bound on the contribution to the energy of a single element. When we set , the time derivative term in (4.28) is

 12ddτ∥Q∥2J,N=12ddτ(JQ,Q)N=12ddτN∑j,k,l=0WjklJjklQTjklQjkl=12N∑j,k,l=0Wjkl{ddτ(JjklQTjkl)Qjkl+JjkldQTjkldτQjkl}=12(∂IN(JQ)∂τ+J∂Q∂τ,Q)N. (5.1)

Therefore,

 12ddτ∥Q∥2J,N+(F∗)TQ∣∣∂E,N−123∑i=1(~Fi,∂∂ξiQ)N−123∑i=1⎛⎜ ⎜⎝Q,∂IN(IN(Ai)Q)∂ξi⎞⎟ ⎟⎠N+123∑i=1⎛⎜ ⎜⎝∂IN(~Ai)∂ξiQ,Q⎞⎟ ⎟⎠N=0. (5.2)

We then apply summation-by-parts (4.21) to the first sum in (5.2) to move the derivative onto the flux

 (5.3)

Since , the two volume terms cancel, leaving

 12ddτ∥Q∥2J,N+(F∗−12~F⋅^n)TQ∣∣∣∂E,N=−123∑i=1⎛⎜ ⎜⎝∂IN(~Ai)∂ξiQ,Q⎞⎟ ⎟⎠N. (5.4)

If the contravariant coefficient matrices are sufficiently smooth so that the derivatives of the interpolants can be bounded, and if the interpolant of the Jacobian is bounded away from zero Gassner:2013uq (), then

 −3∑i=1⎛⎜ ⎜⎝∂IN(~Ai)∂ξiQ,Q⎞⎟ ⎟⎠N⩽2^γ∥Q∥2J,N, (5.5)

where

 ^γ=12max→ξ∈E∣∣ ∣ ∣∣1J∂IN(~Ai)∂ξi∣∣ ∣ ∣∣. (5.6)

The total energy change is found by summing over all elements. When summed over all elements, the boundary terms along internal faces combine, whereas the boundary terms along physical boundaries do not. Let be the solution on the -th element, . If we call the interior face contributions and the boundary contributions , then

 12ddτNel∑m=1∥Qm∥2J,N+ΣI+ΣB≤2^γ2Nel∑m=1∥Qm∥2J,N. (5.7)

We now compute the boundary contributions. The external boundary contributions are

 ΣB=∑boundaryfaces{N∑r,s=0WrWs(F∗−12~F⋅^n)T→m(r,s)Q→m(r,s)}, (5.8)

where we use the subscript “” to represent the appropriate nodal value on that face. For instance, if the face is the right face of the reference hexahedron then . At each nodal face point, the left state is the computed solution value, and the right state is taken from the exterior of the domain, i.e., from the boundary condition, which we denote as in Sec. 3 by . At each nodal point along a boundary surface,

 (5.9)

To guarantee the right kind of bound and therefore stability, we use the upwind solver () at the physical boundaries. For convenience, let us define the intermediate matrix value . Then the contribution from each boundary point in (5.8) is

 (5.10)

where . Since ,

 (F∗−12