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 LagrangianEulerian formulation to map the governing equations to a fixed reference domain. The approximation is made stable by a discretization of a skewsymmetric formulation of the problem. We prove that the discrete approximation is stable, conservative and, for constant coefficient problems, maintains the freestream preservation property. We also provide details on how to add the new skewsymmetric 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, summationbyparts, moving mesh, arbitrary LagrangianEulerian, energy stable, freestream preservation1 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 LagrangianEulerian (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 (DGSEMALE) that was spectrally accurate in space, freestream 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 DGSEMALE 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 skewsymmetric form to guarantee stability, e.g. gassner2015 (); Gassner:2013uq (); kopriva2015 (). The skewsymmetric form of a problem is usually found by averaging the conservative form and a nonconservative advective form of the equation. This is problematic, because it is not obvious that discretizations of the skewsymmetric form remain conservative. Recent success has been had using diagonal norm summationbyparts (SBP) operators to discretize the spatial derivatives in the skewsymmetric 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 LegendreGaussLobatto points, e.g. gassner_skew_burgers (). The approximation developed here will also use a skewsymmetric 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 freestream preservation. Failure to satisfy freestream 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 highorder 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 wellposedness of the ALE formulation, which presents the target to be approximated. In Sec. 4 we use a skewsymmetric 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 freestream preserving properties of the skewsymmetric DGSEMALE 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 LagrangianEulerian (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
(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
(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
(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
(2.4) 
In the ALE formulation, one maps onto a reference domain, , by a transformation
(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
(2.6) 
From the covariant basis, one can formally define the contravariant basis , multiplied by the Jacobian of the transformation, , by
(2.7) 
The Jacobian itself can be written in terms of the covariant vectors as
(2.8) 
The geometry satisfies the wellknown metric identities Vinokur1974 (), (Thompsonetal1985, , Chpt. III)
(2.9) 
and the Geometric Conservation Law (GCL) Mavriplis:2006mb ()
(2.10) 
Under the transformation (2.5), see, e.g. Minoli:2010rt (), the conservation law (2.1) remains a conservation law
(2.11) 
on the reference domain. If we define the contravariant coefficient matrices,
(2.12) 
where denotes the unit vector in the coordinate direction, then we can rewrite (2.11) in fully conservative form as
(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 nonconservative form of the equations. From the chain rule,
(2.14) 
Applying the GCL (2.10) we find,
(2.15) 
giving us the nonconservative system of equations
(2.16) 
Finally, let us define the matrices
(2.17) 
to let us write the conservative and nonconservative forms of the equations as
(2.18) 
and
(2.19) 
respectively.
3 WellPosedness of the ALE Formulation
Since the stability proof mimics that of wellposedness of the PDE, we first show that the system of PDEs on the moving domain is wellposed 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
(3.1) 
and norm . We denote the space .
We show wellposedness by bounding the time rate of change of the energy in the moving domain , which is equivalent to
(3.2) 
The two terms on the right of (3.2) can be found from the conservative and nonconservative forms of the mapped system. The inner product of the conservative form of the equation (2.18) with the solution is
(3.3) 
Similarly, using the nonconservative form (2.19),
(3.4) 
Adding the conservative (3.3) and nonconservative (3.4) forms together gives
(3.5) 
Also, because the matrices are symmetric, we see that
(3.6) 
so
(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
(3.9) 
where are the cardinal bases and the outward normal on the reference box. Then Gauss’ theorem states that
(3.10) 
Next, we determine the bound
(3.11) 
We then note that under the transformation rules,
(3.12) 
is assumed to be bounded. Therefore,
(3.13) 
where
(3.14) 
Integrating both sides of (3.13) with respect to time over a time interval gives
(3.15) 
We now apply boundary conditions to (3.15). First we note that is diagonalizable. From the definition,
(3.16) 
so the normal flux matrices are given by
(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.18) 
so that
(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 wellposed,
(3.21) 
Finally, if the matrices are also constant, then and the energy never grows,
(3.22) 
4 A Stable DGSEMALE 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 nonoverlapping, 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 nonconservative form
(4.1) 
which is written without applying the metric identities and GCL.
To create the skewsymmetric form, we average the conservative and nonconservative forms Gassner:2013uq () to get
(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
(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 skewsymmetric 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 LegendreGaussLobatto 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 LegendreGaussLobatto nodes is
(4.6) 
where the is the Lagrange interpolating polynomial with nodes at the LegendreGaussLobatto 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
(4.8) 
where the singly subscripted ’s are the onedimensional LegendreGaussLobatto 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 LegendreGaussLobatto quadrature that we will use later. The first is the exactness of the quadrature,
(4.9) 
The second is that for some function ,
(4.10) 
which can be seen directly from the definition.
For the numerical flux (Riemann solver) we use the upwind () or central ()
(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
(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,
(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
(4.16) 
where
(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
(4.18) 
Integrating by parts (to put it into the same equation form as the solution),
(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 LegendreGaussLobatto 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
(4.20) 
Furthermore, the discrete inner product satisfies the summationbyparts (SBP) property gassner2010 (). For some and ,
(4.21) 
so with continuity at the boundaries for , the approximation (4.20) is algebraically equivalent to
(4.22) 
Finally, we can combine the two inner products to get the equivalent statement
(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
(4.24) 
or
(4.25) 
when we expand the definition of . Eq. (4.25) holds at each LegendreGaussLobatto node. Therefore, multiplying by the solution vector, quadrature weight and test function at each node, and summing over all nodes,
(4.26) 
We will call (4.26) the weak discrete geometric conservation law or WDGCL. It is also equivalent to the other forms above.
4.2 The SkewSymmetric Approximation on Moving Meshes
Formally and compactly, provided that the discrete metric identities are satisfied Kopriva:2006er (), the skewsymmetric approximation on moving meshes for the Jacobian and solution is the geometric conservation law
(4.27) 
and solution approximation
(4.28) 
where

.

.

.
5 Properties of the SkewSymmetric Approximation
We now show that the skewsymmetric DGSEMALE 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 skewsymmetric 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
(5.1) 
Therefore,
(5.2) 
We then apply summationbyparts (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
(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
(5.5) 
where
(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
(5.7) 
We now compute the boundary contributions. The external boundary contributions are
(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 ,