- Email: [email protected]
Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
Applied Mathematics and Computation 000 (2015) 1–17
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable David A. Kopriva∗, Gregor J. Gassner Department of Mathematics, The Florida State University, Tallahassee, FL, United States and Mathematical Institute, University of Cologne, Cologne, Germany
a r t i c l e
i n f o
Article history: Available online xxx Keywords: Discontinuous Galerkin spectral element method Gauss–Lobatto Legendre Curved elements Skew-symmetry Metric-identities Energy stability
a b s t r a c t We investigate three effects of the variable geometric terms that arise when approximating linear conservation laws on curved elements with a provably stable skew-symmetric variant of the discontinuous Galerkin spectral element method (DGSEM). We show for a constant coefficient system that the non-constant coefficient problem generated by mapping a curved element to the reference element is stable and has energy bounded by the initial value as long as the discrete metric identities are satisfied. Under those same conditions, the skewsymmetric approximation is also constant state preserving and discretely conservative, just like the original DGSEM. © 2015 Elsevier Inc. All rights reserved.
1. Introduction In this paper, we investigate geometrical effects of curved elements on the stability of the discontinuous Galerkin collocation spectral element method (DGSEM) for linear hyperbolic advection problems. Curved elements introduce spatially dependent metric terms in the equations, which in turn create a variable coefficient problem. Depending on the polynomial degree of the curved element mapping, the integration precision of the Gauss–Lobatto points is not enough to guarantee exact L2 -projections. The result is that stability is not guaranteed when using a standard DGSEM approximation on curved elements. Standard DGSEM approximations are widely in use for linear advection problems, such as, e.g. electromagnetics and optics [1–5] and aeroacoustics [6–9]. The DGSEM is an attractive high order approximation scheme, as it offers low dispersion and dissipation errors [10], which are determined by the local polynomial degree. Besides being a high order method, the discretization is applicable to unstructured grids and only needs direct von Neumann neighbors to compute the spatial residuals, leading to highly efficient parallel computing [11]. The extension of the method to multi-dimensions is straightforward via tensor products on quads or hexahedra, yielding in some sense a semi-structured-unstructured method: the discretization is structured inside the element allowing highly efficient implementations, e.g. [12,13], while the elements themselves are managed in an unstructured way allowing for geometric flexibility. When simulating problems in complex geometries, the individual elements (at least the ones at the domain boundary, possibly more) are curved. To simplify the implementation and to increase overall efficiency with respect to CPU time and memory, elements are interpolated from a polynomial mapping, which is then used to transform the equations onto a unit reference element. It is this transformation and its metric terms [14] that generate spatially dependent coefficients in a constant coefficient advection problem. ∗
Corresponding author. Tel.: +1 850 645 0185. E-mail address: [email protected] (D.A. Kopriva).
http://dx.doi.org/10.1016/j.amc.2015.08.047 0096-3003/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 2
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Linear advection problems with variable coefficients introduce a new issue for nodal approximations, for instance documented in [15], where the authors investigated the aliasing errors associated with the polynomial interpolation of a variable coefficient flux. These aliasing errors can dominate the discretization and consequently lead to an unstable approximation. In [16], we proposed an alternative way to construct provably stable approximations by using ideas from summation-byparts finite difference discretizations, e.g. [17,18]. In this paper we focus on the case of linear advection problems on curved elements and specifically study the influence of the variable geometric terms on such approximations. We show how to add skew-symmetric terms to the standard discretization to get a provably stable DGSEM in three space dimensions. With this formulation we investigate the role of the metric terms and show an interesting connection between stability and the metric identities needed to guarantee free stream preservation. The rest of the paper is organized as follows: As a start we present the continuous analysis of linear hyperbolic systems in skew-symmetric form where we derive the bounds on the energy. Section 3 introduces the skew-symmetric DG formulation and its stability analysis. The application to constant coefficient advection problems on curved elements in Section 4 presents the analysis and the connection to the metric terms arising in the curved element mapping. The conclusion of this work is drawn in the last section. 2. Linear hyperbolic systems in skew-symmetric form We will derive a skew-symmetric DGSEM approximation to a linear hyperbolic system of conservation laws
t + ∇ · f = 0 q
(1)
∈ Rr , and the flux f = to be solved on a fixed three dimensional domain x = (x, y, z) = (x1 , x2 , x3 ) ∈ . In (1), the solution q 3 3 r 3 3 ∂ fi . We i=1 f i xˆi = i=1 Ai (x)qxˆi ∈ R × R is linear. With this notation, the divergence of the flux is defined as ∇ · f = i=1 ∂ x i
assume that the r × r coefficient matrices, Ai , have bounded derivatives and can be simultaneously symmetrized. Finally, we assume that appropriate dissipative boundary conditions are applied at physical boundaries, the exact form of which we will specify later. The spectral element approximation [19,12] begins by subdividing the domain into non-overlapping and conforming hexahedral elements, em , such that the union of the elements is . Each element is mapped individually to a reference element ξ , where ξ = (ξ , η, ζ ) = ξ 1 , ξ 2 , ξ 3 is the reference space coordinate. This map is typically E = [−1, 1]3 by a mapping x = X an isoparametric transfinite linear blending of the surfaces of the element [20], where the faces are represented by polynomials i , i = 1, 2, 3 of the same degree as the solution [21]. Under this mapping, we can define a set of contravariant basis vectors, a through
i = J ∇ξ i = Ja
∂ X ∂ X × , i = 1, 2, 3; i, j, k cyclic, j ∂ξ ∂ξ k
(2)
where J ξ > 0 is the Jacobian of the transformation. The conservation law (1) mapped onto the reference element E then becomes
t + q
3 1 J i=1
∂ J ai · f = 0, ∂ξ i
(3)
where
i · f = Ja
3
i · xˆn fn . Ja
(4)
n=1
Under the transformation, therefore, the system of conservation laws becomes another conservation law,
t + q
1 ∇ · f˜ = 0, J ξ
(5)
, the divergence can be written as Since we are assuming a linear flux here, with covariant flux components fi = Ai (x)q
∇ξ · f˜ =
3 i=1
∂ A˜ i q , ∂ξ i
(6)
where
i · A˜ i = J a
3
A j xˆ j
(7)
j=1
are the contravariant coefficient matrices. Note that even if the Aj are constant, the contravariant coefficient matrices will not be unless the volume weighted contravariant basis vectors, J ai , are constant. Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
3
The energy method proofs for well-posedness and stability rely on the coefficient matrices being self-adjoint [22, Sec. 5.1], so we will symmetrize the system (5). We assume that there exists a sufficiently smooth matrix B with ||B|| < ∞ and B−1 < ∞ such that
T
A˜ is = B−1 A˜ i B = A˜ is
,
i = 1, 2, 3.
(8)
The divergence of the flux can be re-written in terms of the symmetric contravariant coefficient matrices through 3 i=1
3 3 3 ∂ A˜ i q ∂ BA˜ is B−1 q ∂ A˜ is B−1 q ∂ B ˜ i −1 A = = B . q + B ∂ξ i ∂ξ i ∂ξ i s ∂ξ i i=1 i=1 i=1
(9)
and note that with the Ai autonomous, B is not a function of time, to write (5) as We now define q¯ = B−1 q
i 3 1 ∂ A˜ s q¯ ∂ B ˜i A q¯ + B = 0. J ∂ξ i s ∂ξ i i=1
3 1 Bq¯t + J i=1
(10)
If we pre-multiply by B−1 then
q¯t +
3 ∂ A˜ is q¯ 1 1 ∂ B ˜i ¯ A q¯ ≡ Sq, = − B−1 i i s J J ∂ξ ∂ξ i=1
3 1 J i=1
(11)
where 3
S = −B−1
i=1
∂ B ˜i A. ∂ξ i s
(12)
(If B does depend on time, it just adds another term to S.) The coefficient matrices are now symmetric, and the right hand side of (11) represents lower order terms that contain the ¯ we arrive at the final form of derivatives of the symmetrization matrix. If we define the symmetric contravariant flux f˜si = A˜ is q, the symmetric conservation law that we want to approximate on the reference element
q¯t +
3 1 J
1 ∂ f˜si ¯ = Sq. ∂ξ i J
i=1
(13)
The skew-symmetric approximation is derived from a combination of the conservation form (13) and a non-conservation form of the PDE [16]. Applying the chain rule, 3 i=1
3 3 ∂ A˜ i ∂ A˜ is q¯ ∂ q¯ s ¯ A˜ is i . = q + i ∂ξ i ∂ξ ∂ξ i=1 i=1
(14)
Then we can write the system of Eq. (13) also in the non-conservative form
1 q¯t + J
3 ∂ A˜ is i=1
∂ξ i
q¯ +
3 1 1 ˜ i ∂ q¯ ¯ As i = Sq. J J ∂ξ
(15)
i=1
The approximation will start with a weak form of the original equation, so we write weak forms of (13) and (15). Let
(u , v)E =
E
T vdξ u
(16)
E = (u , u )E be the associated norm. When we multiply Eqs. (13) and be the L2 inner product on the reference element and u ∈ L2 (E ) and integrate over the volume, we construct the two weak forms of the original system (15) each by a test function φ
∂ A˜ is q¯ + ¯ φ J q¯t , φ ,φ = Sq, i E E ∂ξ i=1 E i 3 3 ∂ A˜ s i ∂ q¯ , Sq¯ . ˜ φ , J q¯t E + φ , As i + φ , = φ q¯ i E ∂ξ ∂ξ E E i=1 i=1
3
(17)
If we average the two equations in (17), see [16], and use the symmetry of the inner product, we get the skew-symmetric weak form
J q¯t , φ
E
i ˜ ∂ A˜ is q¯ ¯ ∂ q , A˜ i , ∂ As q¯ . ¯ φ ,φ + φ + φ = Sq, s i i i E ∂ξ ∂ξ ∂ξ E E i=1
3 1 + 2
(18)
E
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 4
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
We then integrate the first term in the braces of (18) by parts
3 i T 1 i ∂ A˜ is q¯ ∂ φ dS − A˜ i q, ˜ s q¯ φ ¯ A , φ = , s ∂ξ i ∂ξ i E −1 ∂ Ei i=1 i=1
3
(19)
E
where ∂ Ei is the boundary of the reference element. The second term in the braces is also integrated by parts,
1 3 3 ∂ A˜ is φ T i ∂ q¯ i ∂ q¯ i i ˜ ˜ ˜ ¯ As φ , A φ , As i = = φ dS − , q¯ . q s ∂ξ E i=1 ∂ξ i E i=1 ∂ E i ∂ξ i −1 i=1
3
(20)
E
Replacing these terms, and recalling that A˜ is q¯ = f˜si is the symmetric contravariant flux, yields
J q¯t , φ
E
3 i T 1 i ∂ φ dS − 1 f˜s φ f˜si , 2 ∂ξ i E −1 ∂ Ei i=1 i=1 i 3 ˜ is A˜ s φ ∂ A 1 ∂ . ¯ φ + φ , − , q¯ = Sq, q¯ i i E 2 ∂ξ ∂ξ E i=1 +
3
(21)
E
or
J q¯t , φ
where
E
+
1 ¯ φ = 0, BE q, 2
(22)
3
3 i T 1 i ∂(A˜ is φ) i ∂φ dS − ˜f φ ˜ ¯ As i , q¯ + q, s ∂ξ i ∂ξ E −1 ∂ Ei E i=1 i=1 3 ∂ A˜ is . ¯ φ + φ , − 2 Sq, q¯ i E ∂ξ E i=1
=2 ¯ φ BE q,
(23)
2.1. Energy bounds on the skew-symmetric equations To ensure later that the energy growth of the numerical approximation matches that of the analytical solution, we find bounds on the latter. We are only interested in the form of the growth factor in the energy bound, so in this section we assume that the physical domain, , can be mapped by a single map X: E → . The idea is to show that the two norm (energy) of the exact solution can be bounded by the initial energy times a factor that grows at most exponentially fast. The growth factor in the exponent will depend only on the derivatives of the variable coefficient matrices and the source term S. The energy bounds rely on the fact that the coefficient matrices have been symmetrized so that the volume terms containing spatial derivatives of the solution vanish, leaving only an ODE for the energy to be integrated. We first get bounds on the spatial derivative and boundary terms in (22), i.e., BE . Theorem 1. For a simply connected domain mapped onto a reference domain E, if the boundary values are specified so that the boundary conditions are dissipative, i.e.
∂
Ei
i T 1 i f˜s q¯ dS = −1
∂
Ei
1 q¯ T A˜ is q¯ dSi 0, −1
i = 1, 2, 3
(24)
then
BE (v, v) −γ v2 , where
(25)
3 1 ∂ A˜ is γ = max − 2 S . i x∈ J ∂ξ i=1
(26)
2
in (23) by v we get Proof. If we replace q¯ and φ
BE (v, v) = 2
3 i=1
∂ Ei
1 vT A˜ isv dSi − −1
3 i=1
v ∂(A˜ isv) i ∂ ˜ ,v + v , As i ∂ξ i ∂ξ E E
i 3 ∂ A˜ s v, v + − (2 Sv, v)E . ∂ξ i E i=1
(27)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Integration by parts for the first term in the second sum of (27) yields 3
BE (v, v) =
∂ Ei
i=1
1
vT A˜ isv
−1
dSi +
3
A˜ is v,
i=1
i ∂ A˜ s v, v + − (2 Sv, v)E . ∂ξ i E i=1
∂v ∂ξ i
− v, A˜ is E
∂v ∂ξ i
5
E
3
(28)
The symmetry of the coefficient matrices A˜ is allows the second sum in (28) to vanish, so
BE (v, v) =
3
∂ Ei
i=1
T
v
v
1 A˜ is dSi −1
+ v,
3 ∂ A˜ is
∂ξ i
i=1
− 2 S v
,
(29)
E
or
BE (v, v) =
3
∂ Ei
i=1
where
M=
3 ∂ A˜ is i=1
∂ξ i
1 vT A˜ isv dSi + (v, Mv)E ,
(30)
−1
− 2S .
(31)
Since M can be positive or negative, we can bound the inner product with it from below by
M (v, Mv)E − max (v, J v)E = −γ v2 ,
(32)
M γ = max .
(33)
x∈
where
x∈
J
J
2
2
With the assumption that the boundary conditions are dissipative, the required result follows. Corollary 1. If the conditions of Theorem 1 hold and the original problem is constant coefficient, then
B(v, v) 0
(34)
Proof. When the original problem is constant coefficient, M ≡ 0, and the result follows. A consequence to Theorem 1 is that showing well-posedness [22, Section 7.3] is a short operation. Theorem 2. For the same simply connected domain mapped onto a reference domain E as in Theorem 1, the skew-symmetric hyperbolic system of PDEs with dissipative boundary conditions is well-posed, i.e.,
q Ceαt q0 . = q¯ so that Proof. In (22) we set φ
d ¯ q¯ ) = 0. q¯ 2 + BE (q, dt
(35)
From Theorem 1 and the fact that J q¯t , φ
E
= q¯t , φ
,
d ¯ q¯ ) γ q¯ 2 . q¯ 2 = −BE (q, dt
(36)
q¯ eαt q¯ 0 ,
(37)
Then
where α = γ /2. The matrix B and its inverse have bounded norms so that [16]
q Ceαt q0 .
(38)
Remark 1. If the original problem is constant coefficient, then from Corollary 1 and the energy bound in Theorem 2, (38) becomes
q C q0 .
(39)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 6
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
3. A skew-symmetric DG approximation Eq. (21) is the form of the equation that we now approximate on each element. Let PN represent the space of polynomials of degree less than or equal to N in each space direction on an element and PNpw be the space of piecewise polynomials over the domain. For simplicity of exposition, we assume the same order polynomial in each direction, but this is not necessary in practice. We approximate the volume weighted solution by a polynomial J Q¯ ∈ PN and write it in the Lagrange form
JQ¯ =
N
(JQ¯ )i, j,k i (ξ ) j (η) k (ζ ),
(40)
i, j,k=0
where the i , j , k are the one-dimensional Lagrange interpolating polynomials and J ∈ PN is the polynomial interpolation of the Jacobian of the transformation. Similarly, we approximate all other quantities by their Lagrange form interpolants. For instance, the flux is approximated by N
f˜si ≈ F˜si =
f˜si Q¯ i, j,k , ξi, j,k i (ξ ) j (η) k (ζ ).
(41)
i, j,k=0
is also restricted to the piecewise polynomial space. We take φ ∈ PN , for which the basis on each element is The test function φ pw the product i j k . Next we approximate all integrals by Gauss–Lobatto quadrature, which we write as
N
N
udξ ≡
ui, j,k wi w j wk ,
(42)
i, j,k=0
where the w’s are the Gauss–Lobatto quadrature weights and ui, j,k = u ξi , η j , ζk with the nodes being the Gauss–Lobatto nodes. When approximating the inner products with Gauss–Lobatto quadrature we write ( ·, ·)N . As a discontinuous Galerkin method, we replace boundary fluxes by a numerical flux to couple the elements. We consider here either the exact upwind flux (λ = 1) or a centered flux (λ = 0),
F˜si,∗
L
R
Q ,Q
A˜ s R 1 i L R . L − Q = +λ F˜s Q + F˜si Q Q 2 2 i
(43)
Finally, we have options for how to compute the derivatives of the coefficient matrices expressed in (21). We will approximate the derivatives in the reference space by representing the coefficient matrices as we do other quantities as polynomial interpolants
∂ A˜ is ∂ IN A˜ is ≈ , ∂ξ i ∂ξ i
(44)
where IN : L2 (E ) → PN is the polynomial interpolation operator. An alternative is to compute the derivatives analytically in physical space, if that is practical. The elemental skew-symmetric spectral approximation is therefore derived from (21) as
JQ¯t , φ
N
3 T 1 i 1 ∂ φ φ dS − F˜si , 2 ∂ξ i N −1 ∂ E i ,N i=1 i=1 3 ∂ IN A˜ is ∂ IN A˜ is φ 1 . ¯ ¯ φ + φ , − , Q¯ = IN (S)Q, Q i i N 2 ∂ξ ∂ξ i=1 +
3
F˜si,∗
N
(45)
N
To re-write (45) in a form suitable for computation, see [12].
3.1. Stability of the skew-symmetric discontinuous Galerkin approximation The skew-symmetric approximation is stable, and has energy bounds that approximate those of the analytical solution. The procedure mirrors that of finding energy bounds on the analytic solution. Therefore we re-write the approximation (45) on element em as
JQ¯t , φ
N
+
1 m ¯ φ = 0, B Q, 2 N
(46)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
where
Bm N
=2 ¯ φ Q,
3 ∂ E i ,N
i=1 3
+
i=1
T F˜si,∗
1
φ
−1
dS
i
3
−
∂ φ ∂ξ i N
7
F˜si ,
i=1
N i ∂ IN A˜ is ∂ I A˜ s φ φ , − , Q¯ Q¯ i ∂ξ ∂ξ i
¯ φ − 2 IN (S)Q,
N
N
.
(47)
N
We then start with the discrete equivalent to Theorem 1: Theorem 3. If the boundary values are specified so that the approximation at the boundaries is dissipative, i.e.
3
∂ E i ,N
i=1
F˜si,∗
T
1 i T 1 − F˜ V dSi 2 s −1
0
(48)
∈ PN along physical boundaries, then for V pw
,V = BN V
m m −γ V 2 Bm N V ,V
where
(49)
J,N
m
m
3 N ˜i Mm A I ∂ s , Mm = γ = max − 2IN (S) m,x J m ∂ξ i 2
(50)
i=1
and Nel 2 m 2 V = V J,N J,N
(51)
m=1
is the discrete norm whose elemental contributions are N m 2 m 2 V = wi w j wk . Ji,mj,k V i, j,k J,N
(52)
i, j,k=0
Remark 2. Boundary approximations that are dissipative with the Riemann solver (43) are found by setting trivial external state conditions. Strong stability conditions when setting nontrivial external states can be found in [16]. Proof. Replacing the arguments,
Bm N
m = 2 m, V V
3 i=1
+
3
∂ E i ,N
m, V
i=1
i
−1
∂ I A˜ s m V ∂ξ i N
1
T m F˜si,∗ V
dS
i
−
3
i=1
− N
∂I
N
F˜si ,
A˜ isV
∂ V m ∂ξ i m
N
m m, V − 2IN (S)V
m ,V
∂ξ i
N
.
(53)
N
m = IN A˜ i V m . After summation by parts, the volume terms of the We then use summation by parts [23,24] and note that F˜si V s second sum cancel the second term of the third sum leaving
m
m, V BN V
m
=2
3
∂ E i ,N
i=1
i,∗ T
F˜s
When we sum over all elements,
,V = BN V
Nel
Bm N
m
m
V ,V
= 2{ I + B } +
m=1
where
I =
interior faces
1 i T m 1 − F˜ V dSi 2 s −1
Nel
m
V ,
m=1
∂ E i ,N
F˜si,∗
T
−1 V 2
+
m, V
∂ξ i
T 1 i i ˜ Fs V dS
∂ξ i
− 2IN (S)
N
− 2I
m V
.
(54)
N
m
3 ∂ IN A˜ is
m
3 ∂ IN A˜ is i=1
i=1
−1
(S)
m
,
V
(55)
N
(56)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 8
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
represents the contribution of the interior interfaces in terms of the jumps in quantities represented by [ · ] and
B =
boundary faces
∂ E i ,N
F˜si,∗
T
−
1 i T 1 F˜ V dSi 2 s −1
(57)
contains the contributions for the flux at physical boundaries. From [16],
1 λ (F ∗ )T [V ] = [F T V ] + [V ]2 .
(58)
2
2
1 2
λ 2 [V ] ≥ 0
Therefore,
(F ∗ )T [V ] − [F T V ] =
(59)
2
and I ≥ 0 for either the upwind or central Riemann solvers. So
,V 2 B + BN V
Nel m m m ,M V V
(60)
N
m=1
where
m
M =
m
3 ∂ IN A˜ is
∂ξ i
i=1
N
− 2I
(S)
(61)
is the matrix element in element em . We then bound the inner product in (60) from below by
m , M mV m V
where
N
2
m , −γ m V
(62)
J,N
Mm γ m = max m . J
ξ ∈E
(63)
2
Next, we can bound globally as
,V 2 B − γˆ V 2 , BN V
(64)
γˆ = max γ m .
(65)
J,N
where m
With dissipative boundary approximations, B ≥ 0, we are left with
,V −γˆ V 2 . BN V
(66)
J,N
Remark 3. As in the analytical case, note that if
BN
Mm
= 0 then γˆ = 0 and
,V 0. V
(67)
In parallel with the well-posedness result, with Theorem 3, we can show the stability result Theorem 4. With dissipative boundary conditions, the skew-symmetric discontinuous Galerkin spectral element approximation is stable, i.e.,
Q Ceαˆ t Q0 ,
where C and αˆ are constants independent of N and t. = Q¯ so that on element em Proof. Starting from (46), we set φ
m
J Q¯tm , Q¯ m
N
+
1 m m m = 0. B Q¯ , Q¯ 2 N
(68)
Then by Theorem 3 Nel m 2 2 d Q¯ = d Q¯ 2 = −BN Q, ¯ Q¯ γˆ Q¯ . J,N J,N J,N dt dt
(69)
m=1
Therefore,
2 Q¯ eγˆ t Q¯ 0 2 . J,N J,N
(70)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
9
The discrete and continuous norms are equivalent [25], which leads to the bound
Q¯ C1 eαˆ t Q¯ 0 ,
(71)
where αˆ = γˆ /2 and C1 is a constant independent of N and t. Finally, we use the fact that the symmetrization matrix, B, and its inverse are bounded [16] again to see that
Q C2 eαˆ t Q0 ,
(72)
where C2 is also a constant independent of N and t.
Remark 4. The growth rate, αˆ = γˆ /2, for the approximation is the same as the growth rate in the energy estimate for the exact solution (33) to within interpolation errors. If analytical derivatives are used to compute the terms that go into M, then the growth rate of energy of the approximate solution is identical to that of the analytical. Otherwise, we see how close the two are depends on the convergence of an interpolant and its derivatives in the maximum norm, which puts strong conditions on the smoothness of the coefficient and symmetrization matrices. From Remark 3, it follows that Corollary 2. If Mm = 0, then
d Q 0 J,N dt
(73)
and the approximation is stable and the energy is bounded by the initial energy,
Q C Q0 .
(74)
4. Approximation of constant coefficient systems on curved elements As seen in (7), having curved elements in a mesh leads to a non-constant coefficient problem even if the original system of equations has constant coefficients when J ai ξ = const. It is interesting, then, to see what affects the mapping alone has on the approximate solution, and conversely, what properties the transformations from the reference element must have to retain stability. For constant coefficient problems the energy of the skew-symmetric approximation should be bounded by the initial energy, as it is for the analytic solution (Remark 1). Constant coefficient problems also preserve constant states – “free stream preservation” in fluid dynamics parlance. We show that the skew symmetric approximation preserves constant states under the same conditions as for stability. The influence of the element transformation on the differential equations is seen in (6) and (7). If the coefficient matrices of the system and the state = c are constant, then q
∇
3 1 · f = ∇ξ · f˜ = J i=1
3 ∂ Jai · A cxˆ . ∂ξ i j=1 j j
(75)
Analytically, the metric terms satisfy three metric identities [26, pg. 158], in component form written as 3 i=1
∂ Jain = 0 n = 1, 2, 3, ∂ξ i
(76)
= c, q t = −∇ f˜ = 0 and the where Jain = J ai · xˆn is the component of J ai in the nth Cartesian coordinate direction. Then for q constant state is preserved by the PDE on the reference domain as well as in physical space. It was shown in [14] that constant states are preserved in the standard discontinuous Galerkin spectral element approximation iff the interpolant of the metric terms is divergence free, that is,
3 ∂ IN Jain i=1
∂ξ i
=0
n = 1, 2, 3.
(77)
We will show that the skew-symmetric approximation is also constant state preserving if this condition holds. 4.1. Stability of the skew-symmetric approximation for constant coefficient systems We start by showing that stability is ensured, and the energy norm does not grow for the skew symmetric approximation, when the discrete metric identities (77) are satisfied. Theorem 5. If the boundary conditions are dissipative, the skew-symmetric approximation (45) of the constant coefficient system is stable and
Q C Q0
(78)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 10
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
if the discrete metric identities,
3 ∂ IN Jain
=0
∂ξ i
i=1
n = 1, 2, 3
(79)
are satisfied. Proof. By Corollary 2, we have the result (78) if we can show that
3
Mm =
i=1
m ∂ IN A˜ is − 2IN S = 0. ∂ξ i
(80)
When the matrices Ai are constant, the symmetrization matrix B is also constant, so the matrix
S = B−1
3 i=1
∂ B ˜i A = 0. ∂ξ i s
(81)
Also when B is constant, 3
i=1
Now, 3
i=1
3 ∂ IN A˜ is ∂ IN A˜ i −1 = B B. ∂ξ i ∂ξ i i=1
(82)
3 3 ∂ IN A˜ i ∂ = IN Jain An . ∂ξ i ∂ξ i n=1 i=1
(83)
Using the fact that the An are constants, 3
i=1
3 3 ∂ IN Jain ∂ IN A˜ is −1 =B An ∂ξ i ∂ξ i n=1 i=1
(84)
Therefore, if the sums,
3 ∂ IN Jain i=1
∂ξ i
=0
n = 1, 2, 3,
(85)
i.e., the discrete metric identities hold, then Mm = 0.
4.1.1. Example As an example of the stabilization that the use of the skew-symmetric approximation provides, we approximate a constant solution to the symmetric wave equation,
⎡ ⎤ p
⎡
0
⎣u⎦ + ⎣ c v t 0
⎤⎡ ⎤ p
⎡
⎤⎡ ⎤
c
0
0
0
c
0
0⎦⎣u⎦ + ⎣0
0
0⎦⎣u⎦ = 0,
0
0
0
0
v
x
c
p
v
(86)
y
where c = 1 is the constant wave speed, on the mesh shown in Fig. 1. The boundaries are approximated as polynomials of degree five. We choose the initial and external boundary states so that the solution is the Gaussian plane wave
⎡ ⎤
⎡
1
⎤
⎢ ⎥ p ⎢ kx ⎥ (kx x2 +ky y2 −ct )2 ⎣u⎦ = ⎢ c ⎥e− w2 ⎢ ⎥ ⎣ ⎦ v ky
(87)
c
√ with kx = ky = 1/ 2 and w = 1.6. t ) as a function of time step number for the central Riemann solver, λ = 0, and N = 5. Fig. 2 shows the maximum residual (Q The original, non-skew-symmetric, DGSEM approximation is weakly unstable in that the maximum residual grows slowly with time, whereas the skew-symmetric approximation formulation stabilizes the approximation. Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
[m3Gsc;September 3, 2015;10:37] 11
Fig. 1. Grid with curved elements for testing.
Log10(Max. Residual)
0
1
2 Skew-Symmetric DGSEM 3
4 0
20000 Time Step
40000
Fig. 2. Residual growth in the DGSEM shows weak instability for the central Riemann solver whereas the skew-symmetric approximation is stable as a plane wave propagates through and then leaves the domain.
4.2. Constant state preservation Next we show that a constant state is preserved by the skew-symmetric approximation of a constant coefficient system if the metric identities hold. We start with the following result that is needed to show that the flux contributions at element interfaces cancel when the solution is constant. Lemma 1. The numerical flux is consistent if the mesh is conforming, i.e.
Q = F˜ i Q . F˜si,∗ Q, s
(88)
Proof. Recall that
ai · F˜si Q¯ = A˜ is Q¯ = B−1 A˜ i BQ¯ = B−1 J
3
¯ A j xˆ j BQ.
(89)
j=1
The numerical flux is
A˜ s L R 1 i L R = R . ,Q + F˜ i Q L − Q +λ F˜si,∗ Q F˜ Q Q s 2 s 2 i
(90)
Then for either the upwind or central numerical flux
Q = 1 F˜ i Q + F˜ i Q . F˜si,∗ Q, s 2 s
(91)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 12
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Expanding out,
i,∗
F˜s
L i
J a
Q = B−1 Q,
R
+ J ai 2
·
3
A j xˆ j BQ.
(92)
j=1
But J ai = ∇ξ i is continuous across element boundaries if the mesh is conforming, because
∇ξ
i
1 = J
∂ X ∂ X × ∂ξ j ∂ξ k
i, j, k cyclic
(93)
contains only derivatives tangent to the surface. For conforming boundaries, these tangential derivatives match exactly and the result follows. Constant state preservation is then guaranteed by Theorem 6. For constant coefficient systems, if the discrete metric identities (79) hold, then the skew-symmetric approximation (45) preserves a constant state Q¯ = c, that is, the time derivative at each nodal point (i, j, k) in any element m, Qtm = 0. i, j,k
Proof. We write the approximation (45) as
JQ¯t , φ
N
where
1 1 1 Term2 + Term3 − Term4 + Term5 = Term1 2 2 2
−
¯ φ Term1 = IN (S)Q, Term2 =
3
∂ φ
F˜si ,
N
∂ξ i N 3 ∂ IN A˜ is Term3 = φ , Q¯ ∂ξ i i=1 N 3 ∂ IN A˜ is φ Term4 = , Q¯ ∂ξ i i=1 N 3 i,∗ T 1 i Term5 = φ dS F˜s i i=1
(94)
∂ E ,N
i=1
(95)
−1
and study each term individually. • Term1 When the coefficients are constant, the matrix B is constant so that
S = B−1
3 i=1
Therefore,
∂ B ˜i A = 0. ∂ξ i s
¯ φ Term1 = IN (S)Q, • Term2 Next, 3
∂ φ ∂ξ i
F˜si ,
i=1
Now,
= N
N
(96)
= 0.
3 ∂ E i ,N
i=1
(97)
3 i ∂ F˜ i T 1 i s ˜ ,φ . Fs φ dS − ∂ξ i −1 N i=1
F˜si = IN A˜ is B−1 c = IN B−1 A˜ i BB−1 c = IN B−1 A˜ i c, so 3 ∂ F˜si i=1
∂ξ i
=
3 i=1
3 ∂ c IN B−1 J ai · A j xˆ j ∂ξ i j=1
(98)
(99)
(100)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
13
Then for constant coefficients, 3 i=1
3 3 3 i ∂ IN Jai ∂ −1 N −1 B I J a · A j xˆ j c = B · A j xˆ j c, i i ∂ξ ∂ξ j=1 i=1 j=1
(101)
and the term in the parentheses vanishes by the assumption that the metric identities are satisfied discretely. Therefore, only the boundary integral term in (98) remains and
Term2 = I1, where 3
I1 =
∂ E i ,N
i=1
(102)
i T 1 i dS . F˜s φ
(103)
−1
• Term3 From (84) and the assumption that the discrete metric identities hold, 3 ∂ IN A˜ is
∂ξ i
i=1
= 0.
(104)
Therefore,
Term3 = 0.
(105)
• Term4 If we sum by parts,
i T 1 i ¯ ∂ IN A˜ is φ dS − A˜ i φ , ∂Q ¯ ˜s φ , Q = . F s ∂ξ i ∂ξ i N −1 ∂ E i ,N
(106)
N
But by hypothesis, Q¯ = Bc so the volume term vanishes, leaving
Term4 = I1.
(107)
At this point, we have
JQ¯t , φ
N
+ Term5 − I1 = 0,
(108)
so the time derivative depends only on the contribution of the boundary terms:
Term5 − I1 =
3 i=1
By Lemma 1,
∂ E i ,N
F˜si,∗
T
T 1 i φ dS .
− F˜si
−1
¯ Q¯ = F˜si Q¯ F˜si,∗ Q,
(109)
(110)
so Term5 − I1 = 0, and
JQ¯t , φ
N
= 0.
(111)
is arbitrary, we can take it to be the Lagrange basis function φ = (ξ ) (η) (ζ )eˆn , where eˆn is the unit vector Finally, since φ i k k whose nth component is one and the rest are zero. The discrete orthogonality then implies that
Ji, j,k Q¯t i.e.,
Q¯t
i, j,k
i, j,k
wi w j wk = 0,
=0
(112)
(113)
so the constant state is preserved and the result follows since B is constant. 4.2.1. Example As an example of what happens when we do or do not satisfy the metric identities, we approximate a constant solution to the symmetric wave equation, (86), on the mesh shown in Fig. 1. A constant is a solution to the wave equation and that constant should remain constant for all time. For this example, we compute the metric terms by
J a1 = Yη xˆ − Xη yˆ ˆ J a2 = −Yξ xˆ + Xξ y.
(114)
= X xˆ + Y y, ˆ is the linear transfinite interpolation [12] of the element edges to the reference element. The mapping, X Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 14
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
2 N=3 2 4
Log10(Max. Error)
Log10(Max Error)
4 6
8
10
N=4
6
8
10 N=5-8
12
12 3
4
5
6
7
8
0
N
20000 Time Step
40000
Fig. 3. Maximum error for the constant state solved on the mesh of Fig. 1 (left) and time history of the maximum error for the initially constant state (right). For N ≥ 5 the discrete metric identities are satisfied.
∈ PN , where N is the order at which It was shown in [14] that the discrete metric identities are satisfied if and only if X the solution is approximated. It was noted there that with isoparametric approximation of the boundaries by polynomials the condition is satisfied, but if the boundaries are approximated by polynomials of higher order than the solution then it is not. Therefore, if the boundaries are represented at a fixed order then the constant state is preserved if and only if the solution is approximated at that order or higher. For the example, we approximate the curved boundaries of the domain in Fig. 1 by polynomials of degree five and set the external boundary states and initial conditions to be a constant. The metric identities are satisfied therefore only for N ≥ 5 and we would not expect the constant solution to satisfy the discrete equations for solution approximation orders less than five. Fig. 3 shows that the constant state is preserved when the discrete metric identities are satisfied, but otherwise not. Fig. 3 shows the maximum error as a function of approximation order for the constant state solution. We see large errors for N = 3 and N = 4, i.e., when the metric identities are not satisfied, but as soon as they are at N = 5, the errors drop almost eight orders of magnitude to the level of rounding. The time histories of the errors computed with the skew-symmetric approximation are also shown in Fig. 3. When the discrete metric identities are not satisfied (N = 3, 4) the constant state is not a solution of the discrete equations and the solution adjusts. For N ≥ 5 the metric identities are satisfied and the time derivatives remain near rounding error for all time. 5. Conservation Finally, we show that the addition of the non-conservative form of the equation to the DGSEM does not change the conservation properties of the approximation. Theorem 7. If the symmetric equation (13) is strongly conservative, i.e., S = 0, the approximation (45) is discretely conservative in that Nel
Jm Q¯tm dξ =
m=1E,N
∂ E i ,N boundary faces
1
F˜si,∗
−1
dSi .
(115)
Remark 5. For constant coefficient PDEs, S = 0 since it depends on derivatives of the covariant coefficient matrices and the derivatives of the symmetrization matrix. = eˆn , then each component of the system satisfies Proof. If we set φ
JQ¯t , eˆn
3
N
3 T 1 i ˆn 1 i ∂e ˜ + eˆn dS − Fs , 2 ∂ξ i N −1 ∂ E i ,N i=1 i=1 3 ∂ IN A˜ is eˆn 1 ∂ IN A˜ is ¯ ¯ eˆn + − , Q¯ = IN (S)Q, eˆn , Q i i N 2 ∂ξ ∂ξ N i=1
F˜si,∗
(116)
N
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
First, when S = 0,
¯ eˆn IN (S)Q,
N
15
= 0.
(117)
Next, it is immediate that 3
∂ eˆn ∂ξ i
F˜si ,
i=1
= 0.
(118)
N
The most important fact is that the inclusion of the “nonconservative” terms have no effect since eˆn is constant and A˜ is is symmetric so that
N i ∂ IN A˜ is ∂ I A˜ s eˆn ∂ IN A˜ is ∂ IN A˜ is − , Q¯ = eˆn , − = 0. eˆn , eˆn , Q¯ Q¯ Q¯ ∂ξ i ∂ξ i ∂ξ i ∂ξ i N
N
N
(119)
N
Therefore, conservation depends only on the surface integral of the fluxes, as required:
JQ¯t , eˆn
N
+
3
∂ E i ,N
i=1
F˜si,∗
T 1 i eˆn dS = 0. −1
(120)
This says that on each element, the solutions satisfy a discrete conservation property element-wise and that N
JQ¯
i jk
wi w j wk
(121)
i jk
is conserved. Finally, global conservation is ensured by the continuity of the numerical fluxes, which cancel leaving only the integral of the fluxes at physical boundaries. 5.1. Example As an example showing the conservation properties of the skew-symmetric approximation, we integrate the wave equation (86) with the initial condition
p(x, y, 0) = exp − ln 2
x2 +y2 2.32
(122)
u(x, y, 0) = v(x, y, 0) = 0 which produces an expanding p wave
p(x, y, t ) = −
∞ 0
e−ω /4b ωJ0 (rω) cos (ωt )dω, 2b 2
(123)
where J0 is the Bessel function of the first kind of order zero, and b = ln 2/w2 with w = 2.3. We compute the solution on the circular mesh shown in Fig. 4 on the left. A snapshot of p is shown in Fig. 4 on the right.
Fig. 4. Mesh and contours of the expanding p wave to test conservation of the skew-symmetric approximation.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 16
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
1.0
Normalized Total
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
Time Fig. 5. Normalized total p for the expanding wave as a function of time showing conservation until the wave is lost through the outer boundary
Fig. 5 shows the normalized total amount of p over the entire domain. That is, we plot Nel N
Ji,mj p¯ m i, j wi w j
m=1 i, j=0
normalized to its initial value as a function of time. We see that the total remains constant until it is lost through the outer mesh boundary, as expected. Fig. 5 supports the result of Theorem 7 that the skew-symmetric approximation is conservative. 6. Conclusions In this work, we analyze the stability of a skew symmetric DGSEM for constant coefficient advection problems on curved domains. The continuous problem is discretized with the novel DGSEM introduced in [16]. We show the importance of the so-called metric identities for the skew-symmetric DGSEM. Satisfying the metric identities is important to guarantee the socalled ‘free-stream preservation’ in a standard DGSEM. We prove that the same condition guarantees free-stream preservation for the new scheme and furthermore that the skew-symmetric formulation being the average of the conservative and the nonconservative form of the equation is discretely conservative when the metric identities hold and a conforming mesh is used. Thus, in conclusion, the new DGSEM formulation satisfies the same properties as the original: free-stream preservation and element-wise conservation but is provably stable on curved elements. Additionally, positive properties of a high order DGSEM, such as the good parallel efficiency, are still valid since the additional skew-symmetric terms depend on local element data only. References [1] G. Cohen, X. Ferrieres, S. Pernet, A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell’s equations in time domain, J. Comput. Phys. 217 (2) (2006) 340–363. [2] D.A. Kopriva, S.L. Woodruff, M.Y. Hussaini, Discontinuous spectral element approximation of Maxwell’s Equations, in: B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Proceedings of the International Symposium on Discontinuous Galerkin Methods, Springer-Verlag, New York, 2000, pp. 355–361. [3] D.A. Kopriva, S.L. Woodruff, M.Y. Hussaini, Computation of electromagnetic scattering with a non-conforming discontinuous spectral element method, Int. J. Numer. Methods Eng. 53 (1) (2002) 105–122. [4] S. Deng, Numerical simulation of optical coupling and light propagation in coupled optical resonators with size disorder, Appl. Numer. Math. 57 (5–7) (2007) 475–485, doi:10.1016/j.apnum.2006.07.001. [5] S. Deng, W. Cai, V. Astratov, Numerical study of light propagation via whispering gallery modes in microcylinder coupled resonator optical waveguides, Opt. Express 12 (26) (2004) 6468–6480. [6] P. Rasetarinera, D. Kopriva, M. Hussaini, Discontinuous spectral element solution of acoustic radiation from thin airfoils, AIAA J. 39 (11) (2001) 2070–2075. [7] D. Stanescu, F. Farassat, M. Hussaini, Aircraft engine noise scattering – parallel discontinuous Galerkin spectral element method, Paper 2002-0800, AIAA, 2002. [8] D. Stanescu, J. Xu, F. Farassat, M. Hussaini, Computation of engine noise propagation and scattering off an aircraft, Aeroacoustics 1 (4) (2002) 403–420. [9] N. Castel, G. Cohen, M. Durufle, Application of discontinuous Galerkin spectral method on hexahedral elements for aeroacoustic, J. Comput. Acoust. 17 (2) (2009) 175–196. [10] G. Gassner, D.A. Kopriva, A comparison of the dispersion and dissipation errors of Gauss and Gauss–Lobatto discontinuous Galerkin spectral element methods, SIAM J. Sci. Comput. 33 (5) (2010) 2560–2579. [11] C. Altmann, A.D. Beck, F. Hindenlang, M. Staudenmaier, G.J. Gassner, C.-D. Munz, An efficient high performance parallelization of a discontinuous Galerkin spectral element method, in: R. Keller, D. Kramer, J.-P. Weiss (Eds.), Facing the Multicore-Challenge III, Lecture Notes in Computer Science, vol. 7686, Springer, Berlin, Heidelberg, 2013, pp. 37–47.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
JID: AMC
ARTICLE IN PRESS D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
[m3Gsc;September 3, 2015;10:37] 17
[12] D.A. Kopriva, Implementing Spectral Methods for Partial Differential Equations, Scientific Computation, Springer, 2009. [13] F. Hindenlang, G.J. Gassner, C. Altmann, A. Beck, M. Staudenmaier, C.-D. Munz, Explicit discontinuous Galerkin methods for unsteady problems, Comput. Fluids 61 (2012) 86–93, doi:10.1016/j.compfluid.2012.03.006. [14] D.A. Kopriva, Metric identities and the discontinuous spectral element method on curvilinear meshes, J. Sci. Comput. 26 (3) (2006) 301–327. [15] J. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer, 2008. [16] D.A. Kopriva, G.J. Gassner, An energy stable discontinuous Galerkin spectral element discretization for variable coefficient advection problems, SIAM J. Sci. Comput. 36 (4) (2014) A2076–A2099. [17] J. Nordström, Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation, J. Sci. Comput. 29 (3) (2006) 375–404, doi:10.1007/s10915-005-9013-4. [18] J. Kozdon, E. Dunham, J. Nordström, Simulation of dynamic earthquake ruptures in complex geometries using high-order finite difference methods, J. Sci. Comput. 55 (1) (2013) 92–124, doi:10.1007/s10915-012-9624-5. [19] M. Deville, P. Fischer, E. Mund, High Order Methods for Incompressible Fluid Flow, Cambridge University Press, 2002. [20] W. Gordon, C. Hall, Construction of curvilinear co-ordinate systems and their applications to mesh generation, Int. J. Numer. Methods Eng. 7 (1973) 461–477. [21] K. Korczak, A. Patera, An isoparametric spectral element method for solution of the Navier–Stokes equations in complex-geometry, J. Comput. Phys. 62 (2) (1986) 361–382. [22] H.-O. Kreiss, J. Lorenz, Initial-Boundary Value Problems and the Navier–Stokes Equations, Pure and Applied Mathematics, 136, Academic Press, San Diego, CA, USA, 1989. [23] D.A. Kopriva, G. Gassner, On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods, J. Sci. Comput. 44 (2) (2010) 136–155. [24] G. Gassner, A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods, SIAM J. Sci. Comput. 35 (3) (2013) A1233–A1253. [25] C. Canuto, A. Quarteroni, Approximation results for orthogonal polynomials in Sobolev spaces, Math. Comput. 38 (157) (1982) 67–86, doi:10.2307/2007465. [26] J. Thompson, Z. Warsi, C. Mastin, Numerical Grid Generation, Elsevier Science Publishing, Amsterdam, 1985.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
JID: AMC
[m3Gsc;September 3, 2015;10:37]
Applied Mathematics and Computation 000 (2015) 1–17
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable David A. Kopriva∗, Gregor J. Gassner Department of Mathematics, The Florida State University, Tallahassee, FL, United States and Mathematical Institute, University of Cologne, Cologne, Germany
a r t i c l e
i n f o
Article history: Available online xxx Keywords: Discontinuous Galerkin spectral element method Gauss–Lobatto Legendre Curved elements Skew-symmetry Metric-identities Energy stability
a b s t r a c t We investigate three effects of the variable geometric terms that arise when approximating linear conservation laws on curved elements with a provably stable skew-symmetric variant of the discontinuous Galerkin spectral element method (DGSEM). We show for a constant coefficient system that the non-constant coefficient problem generated by mapping a curved element to the reference element is stable and has energy bounded by the initial value as long as the discrete metric identities are satisfied. Under those same conditions, the skewsymmetric approximation is also constant state preserving and discretely conservative, just like the original DGSEM. © 2015 Elsevier Inc. All rights reserved.
1. Introduction In this paper, we investigate geometrical effects of curved elements on the stability of the discontinuous Galerkin collocation spectral element method (DGSEM) for linear hyperbolic advection problems. Curved elements introduce spatially dependent metric terms in the equations, which in turn create a variable coefficient problem. Depending on the polynomial degree of the curved element mapping, the integration precision of the Gauss–Lobatto points is not enough to guarantee exact L2 -projections. The result is that stability is not guaranteed when using a standard DGSEM approximation on curved elements. Standard DGSEM approximations are widely in use for linear advection problems, such as, e.g. electromagnetics and optics [1–5] and aeroacoustics [6–9]. The DGSEM is an attractive high order approximation scheme, as it offers low dispersion and dissipation errors [10], which are determined by the local polynomial degree. Besides being a high order method, the discretization is applicable to unstructured grids and only needs direct von Neumann neighbors to compute the spatial residuals, leading to highly efficient parallel computing [11]. The extension of the method to multi-dimensions is straightforward via tensor products on quads or hexahedra, yielding in some sense a semi-structured-unstructured method: the discretization is structured inside the element allowing highly efficient implementations, e.g. [12,13], while the elements themselves are managed in an unstructured way allowing for geometric flexibility. When simulating problems in complex geometries, the individual elements (at least the ones at the domain boundary, possibly more) are curved. To simplify the implementation and to increase overall efficiency with respect to CPU time and memory, elements are interpolated from a polynomial mapping, which is then used to transform the equations onto a unit reference element. It is this transformation and its metric terms [14] that generate spatially dependent coefficients in a constant coefficient advection problem. ∗
Corresponding author. Tel.: +1 850 645 0185. E-mail address: [email protected] (D.A. Kopriva).
http://dx.doi.org/10.1016/j.amc.2015.08.047 0096-3003/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 2
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Linear advection problems with variable coefficients introduce a new issue for nodal approximations, for instance documented in [15], where the authors investigated the aliasing errors associated with the polynomial interpolation of a variable coefficient flux. These aliasing errors can dominate the discretization and consequently lead to an unstable approximation. In [16], we proposed an alternative way to construct provably stable approximations by using ideas from summation-byparts finite difference discretizations, e.g. [17,18]. In this paper we focus on the case of linear advection problems on curved elements and specifically study the influence of the variable geometric terms on such approximations. We show how to add skew-symmetric terms to the standard discretization to get a provably stable DGSEM in three space dimensions. With this formulation we investigate the role of the metric terms and show an interesting connection between stability and the metric identities needed to guarantee free stream preservation. The rest of the paper is organized as follows: As a start we present the continuous analysis of linear hyperbolic systems in skew-symmetric form where we derive the bounds on the energy. Section 3 introduces the skew-symmetric DG formulation and its stability analysis. The application to constant coefficient advection problems on curved elements in Section 4 presents the analysis and the connection to the metric terms arising in the curved element mapping. The conclusion of this work is drawn in the last section. 2. Linear hyperbolic systems in skew-symmetric form We will derive a skew-symmetric DGSEM approximation to a linear hyperbolic system of conservation laws
t + ∇ · f = 0 q
(1)
∈ Rr , and the flux f = to be solved on a fixed three dimensional domain x = (x, y, z) = (x1 , x2 , x3 ) ∈ . In (1), the solution q 3 3 r 3 3 ∂ fi . We i=1 f i xˆi = i=1 Ai (x)qxˆi ∈ R × R is linear. With this notation, the divergence of the flux is defined as ∇ · f = i=1 ∂ x i
assume that the r × r coefficient matrices, Ai , have bounded derivatives and can be simultaneously symmetrized. Finally, we assume that appropriate dissipative boundary conditions are applied at physical boundaries, the exact form of which we will specify later. The spectral element approximation [19,12] begins by subdividing the domain into non-overlapping and conforming hexahedral elements, em , such that the union of the elements is . Each element is mapped individually to a reference element ξ , where ξ = (ξ , η, ζ ) = ξ 1 , ξ 2 , ξ 3 is the reference space coordinate. This map is typically E = [−1, 1]3 by a mapping x = X an isoparametric transfinite linear blending of the surfaces of the element [20], where the faces are represented by polynomials i , i = 1, 2, 3 of the same degree as the solution [21]. Under this mapping, we can define a set of contravariant basis vectors, a through
i = J ∇ξ i = Ja
∂ X ∂ X × , i = 1, 2, 3; i, j, k cyclic, j ∂ξ ∂ξ k
(2)
where J ξ > 0 is the Jacobian of the transformation. The conservation law (1) mapped onto the reference element E then becomes
t + q
3 1 J i=1
∂ J ai · f = 0, ∂ξ i
(3)
where
i · f = Ja
3
i · xˆn fn . Ja
(4)
n=1
Under the transformation, therefore, the system of conservation laws becomes another conservation law,
t + q
1 ∇ · f˜ = 0, J ξ
(5)
, the divergence can be written as Since we are assuming a linear flux here, with covariant flux components fi = Ai (x)q
∇ξ · f˜ =
3 i=1
∂ A˜ i q , ∂ξ i
(6)
where
i · A˜ i = J a
3
A j xˆ j
(7)
j=1
are the contravariant coefficient matrices. Note that even if the Aj are constant, the contravariant coefficient matrices will not be unless the volume weighted contravariant basis vectors, J ai , are constant. Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
3
The energy method proofs for well-posedness and stability rely on the coefficient matrices being self-adjoint [22, Sec. 5.1], so we will symmetrize the system (5). We assume that there exists a sufficiently smooth matrix B with ||B|| < ∞ and B−1 < ∞ such that
T
A˜ is = B−1 A˜ i B = A˜ is
,
i = 1, 2, 3.
(8)
The divergence of the flux can be re-written in terms of the symmetric contravariant coefficient matrices through 3 i=1
3 3 3 ∂ A˜ i q ∂ BA˜ is B−1 q ∂ A˜ is B−1 q ∂ B ˜ i −1 A = = B . q + B ∂ξ i ∂ξ i ∂ξ i s ∂ξ i i=1 i=1 i=1
(9)
and note that with the Ai autonomous, B is not a function of time, to write (5) as We now define q¯ = B−1 q
i 3 1 ∂ A˜ s q¯ ∂ B ˜i A q¯ + B = 0. J ∂ξ i s ∂ξ i i=1
3 1 Bq¯t + J i=1
(10)
If we pre-multiply by B−1 then
q¯t +
3 ∂ A˜ is q¯ 1 1 ∂ B ˜i ¯ A q¯ ≡ Sq, = − B−1 i i s J J ∂ξ ∂ξ i=1
3 1 J i=1
(11)
where 3
S = −B−1
i=1
∂ B ˜i A. ∂ξ i s
(12)
(If B does depend on time, it just adds another term to S.) The coefficient matrices are now symmetric, and the right hand side of (11) represents lower order terms that contain the ¯ we arrive at the final form of derivatives of the symmetrization matrix. If we define the symmetric contravariant flux f˜si = A˜ is q, the symmetric conservation law that we want to approximate on the reference element
q¯t +
3 1 J
1 ∂ f˜si ¯ = Sq. ∂ξ i J
i=1
(13)
The skew-symmetric approximation is derived from a combination of the conservation form (13) and a non-conservation form of the PDE [16]. Applying the chain rule, 3 i=1
3 3 ∂ A˜ i ∂ A˜ is q¯ ∂ q¯ s ¯ A˜ is i . = q + i ∂ξ i ∂ξ ∂ξ i=1 i=1
(14)
Then we can write the system of Eq. (13) also in the non-conservative form
1 q¯t + J
3 ∂ A˜ is i=1
∂ξ i
q¯ +
3 1 1 ˜ i ∂ q¯ ¯ As i = Sq. J J ∂ξ
(15)
i=1
The approximation will start with a weak form of the original equation, so we write weak forms of (13) and (15). Let
(u , v)E =
E
T vdξ u
(16)
E = (u , u )E be the associated norm. When we multiply Eqs. (13) and be the L2 inner product on the reference element and u ∈ L2 (E ) and integrate over the volume, we construct the two weak forms of the original system (15) each by a test function φ
∂ A˜ is q¯ + ¯ φ J q¯t , φ ,φ = Sq, i E E ∂ξ i=1 E i 3 3 ∂ A˜ s i ∂ q¯ , Sq¯ . ˜ φ , J q¯t E + φ , As i + φ , = φ q¯ i E ∂ξ ∂ξ E E i=1 i=1
3
(17)
If we average the two equations in (17), see [16], and use the symmetry of the inner product, we get the skew-symmetric weak form
J q¯t , φ
E
i ˜ ∂ A˜ is q¯ ¯ ∂ q , A˜ i , ∂ As q¯ . ¯ φ ,φ + φ + φ = Sq, s i i i E ∂ξ ∂ξ ∂ξ E E i=1
3 1 + 2
(18)
E
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 4
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
We then integrate the first term in the braces of (18) by parts
3 i T 1 i ∂ A˜ is q¯ ∂ φ dS − A˜ i q, ˜ s q¯ φ ¯ A , φ = , s ∂ξ i ∂ξ i E −1 ∂ Ei i=1 i=1
3
(19)
E
where ∂ Ei is the boundary of the reference element. The second term in the braces is also integrated by parts,
1 3 3 ∂ A˜ is φ T i ∂ q¯ i ∂ q¯ i i ˜ ˜ ˜ ¯ As φ , A φ , As i = = φ dS − , q¯ . q s ∂ξ E i=1 ∂ξ i E i=1 ∂ E i ∂ξ i −1 i=1
3
(20)
E
Replacing these terms, and recalling that A˜ is q¯ = f˜si is the symmetric contravariant flux, yields
J q¯t , φ
E
3 i T 1 i ∂ φ dS − 1 f˜s φ f˜si , 2 ∂ξ i E −1 ∂ Ei i=1 i=1 i 3 ˜ is A˜ s φ ∂ A 1 ∂ . ¯ φ + φ , − , q¯ = Sq, q¯ i i E 2 ∂ξ ∂ξ E i=1 +
3
(21)
E
or
J q¯t , φ
where
E
+
1 ¯ φ = 0, BE q, 2
(22)
3
3 i T 1 i ∂(A˜ is φ) i ∂φ dS − ˜f φ ˜ ¯ As i , q¯ + q, s ∂ξ i ∂ξ E −1 ∂ Ei E i=1 i=1 3 ∂ A˜ is . ¯ φ + φ , − 2 Sq, q¯ i E ∂ξ E i=1
=2 ¯ φ BE q,
(23)
2.1. Energy bounds on the skew-symmetric equations To ensure later that the energy growth of the numerical approximation matches that of the analytical solution, we find bounds on the latter. We are only interested in the form of the growth factor in the energy bound, so in this section we assume that the physical domain, , can be mapped by a single map X: E → . The idea is to show that the two norm (energy) of the exact solution can be bounded by the initial energy times a factor that grows at most exponentially fast. The growth factor in the exponent will depend only on the derivatives of the variable coefficient matrices and the source term S. The energy bounds rely on the fact that the coefficient matrices have been symmetrized so that the volume terms containing spatial derivatives of the solution vanish, leaving only an ODE for the energy to be integrated. We first get bounds on the spatial derivative and boundary terms in (22), i.e., BE . Theorem 1. For a simply connected domain mapped onto a reference domain E, if the boundary values are specified so that the boundary conditions are dissipative, i.e.
∂
Ei
i T 1 i f˜s q¯ dS = −1
∂
Ei
1 q¯ T A˜ is q¯ dSi 0, −1
i = 1, 2, 3
(24)
then
BE (v, v) −γ v2 , where
(25)
3 1 ∂ A˜ is γ = max − 2 S . i x∈ J ∂ξ i=1
(26)
2
in (23) by v we get Proof. If we replace q¯ and φ
BE (v, v) = 2
3 i=1
∂ Ei
1 vT A˜ isv dSi − −1
3 i=1
v ∂(A˜ isv) i ∂ ˜ ,v + v , As i ∂ξ i ∂ξ E E
i 3 ∂ A˜ s v, v + − (2 Sv, v)E . ∂ξ i E i=1
(27)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Integration by parts for the first term in the second sum of (27) yields 3
BE (v, v) =
∂ Ei
i=1
1
vT A˜ isv
−1
dSi +
3
A˜ is v,
i=1
i ∂ A˜ s v, v + − (2 Sv, v)E . ∂ξ i E i=1
∂v ∂ξ i
− v, A˜ is E
∂v ∂ξ i
5
E
3
(28)
The symmetry of the coefficient matrices A˜ is allows the second sum in (28) to vanish, so
BE (v, v) =
3
∂ Ei
i=1
T
v
v
1 A˜ is dSi −1
+ v,
3 ∂ A˜ is
∂ξ i
i=1
− 2 S v
,
(29)
E
or
BE (v, v) =
3
∂ Ei
i=1
where
M=
3 ∂ A˜ is i=1
∂ξ i
1 vT A˜ isv dSi + (v, Mv)E ,
(30)
−1
− 2S .
(31)
Since M can be positive or negative, we can bound the inner product with it from below by
M (v, Mv)E − max (v, J v)E = −γ v2 ,
(32)
M γ = max .
(33)
x∈
where
x∈
J
J
2
2
With the assumption that the boundary conditions are dissipative, the required result follows. Corollary 1. If the conditions of Theorem 1 hold and the original problem is constant coefficient, then
B(v, v) 0
(34)
Proof. When the original problem is constant coefficient, M ≡ 0, and the result follows. A consequence to Theorem 1 is that showing well-posedness [22, Section 7.3] is a short operation. Theorem 2. For the same simply connected domain mapped onto a reference domain E as in Theorem 1, the skew-symmetric hyperbolic system of PDEs with dissipative boundary conditions is well-posed, i.e.,
q Ceαt q0 . = q¯ so that Proof. In (22) we set φ
d ¯ q¯ ) = 0. q¯ 2 + BE (q, dt
(35)
From Theorem 1 and the fact that J q¯t , φ
E
= q¯t , φ
,
d ¯ q¯ ) γ q¯ 2 . q¯ 2 = −BE (q, dt
(36)
q¯ eαt q¯ 0 ,
(37)
Then
where α = γ /2. The matrix B and its inverse have bounded norms so that [16]
q Ceαt q0 .
(38)
Remark 1. If the original problem is constant coefficient, then from Corollary 1 and the energy bound in Theorem 2, (38) becomes
q C q0 .
(39)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 6
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
3. A skew-symmetric DG approximation Eq. (21) is the form of the equation that we now approximate on each element. Let PN represent the space of polynomials of degree less than or equal to N in each space direction on an element and PNpw be the space of piecewise polynomials over the domain. For simplicity of exposition, we assume the same order polynomial in each direction, but this is not necessary in practice. We approximate the volume weighted solution by a polynomial J Q¯ ∈ PN and write it in the Lagrange form
JQ¯ =
N
(JQ¯ )i, j,k i (ξ ) j (η) k (ζ ),
(40)
i, j,k=0
where the i , j , k are the one-dimensional Lagrange interpolating polynomials and J ∈ PN is the polynomial interpolation of the Jacobian of the transformation. Similarly, we approximate all other quantities by their Lagrange form interpolants. For instance, the flux is approximated by N
f˜si ≈ F˜si =
f˜si Q¯ i, j,k , ξi, j,k i (ξ ) j (η) k (ζ ).
(41)
i, j,k=0
is also restricted to the piecewise polynomial space. We take φ ∈ PN , for which the basis on each element is The test function φ pw the product i j k . Next we approximate all integrals by Gauss–Lobatto quadrature, which we write as
N
N
udξ ≡
ui, j,k wi w j wk ,
(42)
i, j,k=0
where the w’s are the Gauss–Lobatto quadrature weights and ui, j,k = u ξi , η j , ζk with the nodes being the Gauss–Lobatto nodes. When approximating the inner products with Gauss–Lobatto quadrature we write ( ·, ·)N . As a discontinuous Galerkin method, we replace boundary fluxes by a numerical flux to couple the elements. We consider here either the exact upwind flux (λ = 1) or a centered flux (λ = 0),
F˜si,∗
L
R
Q ,Q
A˜ s R 1 i L R . L − Q = +λ F˜s Q + F˜si Q Q 2 2 i
(43)
Finally, we have options for how to compute the derivatives of the coefficient matrices expressed in (21). We will approximate the derivatives in the reference space by representing the coefficient matrices as we do other quantities as polynomial interpolants
∂ A˜ is ∂ IN A˜ is ≈ , ∂ξ i ∂ξ i
(44)
where IN : L2 (E ) → PN is the polynomial interpolation operator. An alternative is to compute the derivatives analytically in physical space, if that is practical. The elemental skew-symmetric spectral approximation is therefore derived from (21) as
JQ¯t , φ
N
3 T 1 i 1 ∂ φ φ dS − F˜si , 2 ∂ξ i N −1 ∂ E i ,N i=1 i=1 3 ∂ IN A˜ is ∂ IN A˜ is φ 1 . ¯ ¯ φ + φ , − , Q¯ = IN (S)Q, Q i i N 2 ∂ξ ∂ξ i=1 +
3
F˜si,∗
N
(45)
N
To re-write (45) in a form suitable for computation, see [12].
3.1. Stability of the skew-symmetric discontinuous Galerkin approximation The skew-symmetric approximation is stable, and has energy bounds that approximate those of the analytical solution. The procedure mirrors that of finding energy bounds on the analytic solution. Therefore we re-write the approximation (45) on element em as
JQ¯t , φ
N
+
1 m ¯ φ = 0, B Q, 2 N
(46)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
where
Bm N
=2 ¯ φ Q,
3 ∂ E i ,N
i=1 3
+
i=1
T F˜si,∗
1
φ
−1
dS
i
3
−
∂ φ ∂ξ i N
7
F˜si ,
i=1
N i ∂ IN A˜ is ∂ I A˜ s φ φ , − , Q¯ Q¯ i ∂ξ ∂ξ i
¯ φ − 2 IN (S)Q,
N
N
.
(47)
N
We then start with the discrete equivalent to Theorem 1: Theorem 3. If the boundary values are specified so that the approximation at the boundaries is dissipative, i.e.
3
∂ E i ,N
i=1
F˜si,∗
T
1 i T 1 − F˜ V dSi 2 s −1
0
(48)
∈ PN along physical boundaries, then for V pw
,V = BN V
m m −γ V 2 Bm N V ,V
where
(49)
J,N
m
m
3 N ˜i Mm A I ∂ s , Mm = γ = max − 2IN (S) m,x J m ∂ξ i 2
(50)
i=1
and Nel 2 m 2 V = V J,N J,N
(51)
m=1
is the discrete norm whose elemental contributions are N m 2 m 2 V = wi w j wk . Ji,mj,k V i, j,k J,N
(52)
i, j,k=0
Remark 2. Boundary approximations that are dissipative with the Riemann solver (43) are found by setting trivial external state conditions. Strong stability conditions when setting nontrivial external states can be found in [16]. Proof. Replacing the arguments,
Bm N
m = 2 m, V V
3 i=1
+
3
∂ E i ,N
m, V
i=1
i
−1
∂ I A˜ s m V ∂ξ i N
1
T m F˜si,∗ V
dS
i
−
3
i=1
− N
∂I
N
F˜si ,
A˜ isV
∂ V m ∂ξ i m
N
m m, V − 2IN (S)V
m ,V
∂ξ i
N
.
(53)
N
m = IN A˜ i V m . After summation by parts, the volume terms of the We then use summation by parts [23,24] and note that F˜si V s second sum cancel the second term of the third sum leaving
m
m, V BN V
m
=2
3
∂ E i ,N
i=1
i,∗ T
F˜s
When we sum over all elements,
,V = BN V
Nel
Bm N
m
m
V ,V
= 2{ I + B } +
m=1
where
I =
interior faces
1 i T m 1 − F˜ V dSi 2 s −1
Nel
m
V ,
m=1
∂ E i ,N
F˜si,∗
T
−1 V 2
+
m, V
∂ξ i
T 1 i i ˜ Fs V dS
∂ξ i
− 2IN (S)
N
− 2I
m V
.
(54)
N
m
3 ∂ IN A˜ is
m
3 ∂ IN A˜ is i=1
i=1
−1
(S)
m
,
V
(55)
N
(56)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 8
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
represents the contribution of the interior interfaces in terms of the jumps in quantities represented by [ · ] and
B =
boundary faces
∂ E i ,N
F˜si,∗
T
−
1 i T 1 F˜ V dSi 2 s −1
(57)
contains the contributions for the flux at physical boundaries. From [16],
1 λ (F ∗ )T [V ] = [F T V ] + [V ]2 .
(58)
2
2
1 2
λ 2 [V ] ≥ 0
Therefore,
(F ∗ )T [V ] − [F T V ] =
(59)
2
and I ≥ 0 for either the upwind or central Riemann solvers. So
,V 2 B + BN V
Nel m m m ,M V V
(60)
N
m=1
where
m
M =
m
3 ∂ IN A˜ is
∂ξ i
i=1
N
− 2I
(S)
(61)
is the matrix element in element em . We then bound the inner product in (60) from below by
m , M mV m V
where
N
2
m , −γ m V
(62)
J,N
Mm γ m = max m . J
ξ ∈E
(63)
2
Next, we can bound globally as
,V 2 B − γˆ V 2 , BN V
(64)
γˆ = max γ m .
(65)
J,N
where m
With dissipative boundary approximations, B ≥ 0, we are left with
,V −γˆ V 2 . BN V
(66)
J,N
Remark 3. As in the analytical case, note that if
BN
Mm
= 0 then γˆ = 0 and
,V 0. V
(67)
In parallel with the well-posedness result, with Theorem 3, we can show the stability result Theorem 4. With dissipative boundary conditions, the skew-symmetric discontinuous Galerkin spectral element approximation is stable, i.e.,
Q Ceαˆ t Q0 ,
where C and αˆ are constants independent of N and t. = Q¯ so that on element em Proof. Starting from (46), we set φ
m
J Q¯tm , Q¯ m
N
+
1 m m m = 0. B Q¯ , Q¯ 2 N
(68)
Then by Theorem 3 Nel m 2 2 d Q¯ = d Q¯ 2 = −BN Q, ¯ Q¯ γˆ Q¯ . J,N J,N J,N dt dt
(69)
m=1
Therefore,
2 Q¯ eγˆ t Q¯ 0 2 . J,N J,N
(70)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
9
The discrete and continuous norms are equivalent [25], which leads to the bound
Q¯ C1 eαˆ t Q¯ 0 ,
(71)
where αˆ = γˆ /2 and C1 is a constant independent of N and t. Finally, we use the fact that the symmetrization matrix, B, and its inverse are bounded [16] again to see that
Q C2 eαˆ t Q0 ,
(72)
where C2 is also a constant independent of N and t.
Remark 4. The growth rate, αˆ = γˆ /2, for the approximation is the same as the growth rate in the energy estimate for the exact solution (33) to within interpolation errors. If analytical derivatives are used to compute the terms that go into M, then the growth rate of energy of the approximate solution is identical to that of the analytical. Otherwise, we see how close the two are depends on the convergence of an interpolant and its derivatives in the maximum norm, which puts strong conditions on the smoothness of the coefficient and symmetrization matrices. From Remark 3, it follows that Corollary 2. If Mm = 0, then
d Q 0 J,N dt
(73)
and the approximation is stable and the energy is bounded by the initial energy,
Q C Q0 .
(74)
4. Approximation of constant coefficient systems on curved elements As seen in (7), having curved elements in a mesh leads to a non-constant coefficient problem even if the original system of equations has constant coefficients when J ai ξ = const. It is interesting, then, to see what affects the mapping alone has on the approximate solution, and conversely, what properties the transformations from the reference element must have to retain stability. For constant coefficient problems the energy of the skew-symmetric approximation should be bounded by the initial energy, as it is for the analytic solution (Remark 1). Constant coefficient problems also preserve constant states – “free stream preservation” in fluid dynamics parlance. We show that the skew symmetric approximation preserves constant states under the same conditions as for stability. The influence of the element transformation on the differential equations is seen in (6) and (7). If the coefficient matrices of the system and the state = c are constant, then q
∇
3 1 · f = ∇ξ · f˜ = J i=1
3 ∂ Jai · A cxˆ . ∂ξ i j=1 j j
(75)
Analytically, the metric terms satisfy three metric identities [26, pg. 158], in component form written as 3 i=1
∂ Jain = 0 n = 1, 2, 3, ∂ξ i
(76)
= c, q t = −∇ f˜ = 0 and the where Jain = J ai · xˆn is the component of J ai in the nth Cartesian coordinate direction. Then for q constant state is preserved by the PDE on the reference domain as well as in physical space. It was shown in [14] that constant states are preserved in the standard discontinuous Galerkin spectral element approximation iff the interpolant of the metric terms is divergence free, that is,
3 ∂ IN Jain i=1
∂ξ i
=0
n = 1, 2, 3.
(77)
We will show that the skew-symmetric approximation is also constant state preserving if this condition holds. 4.1. Stability of the skew-symmetric approximation for constant coefficient systems We start by showing that stability is ensured, and the energy norm does not grow for the skew symmetric approximation, when the discrete metric identities (77) are satisfied. Theorem 5. If the boundary conditions are dissipative, the skew-symmetric approximation (45) of the constant coefficient system is stable and
Q C Q0
(78)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 10
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
if the discrete metric identities,
3 ∂ IN Jain
=0
∂ξ i
i=1
n = 1, 2, 3
(79)
are satisfied. Proof. By Corollary 2, we have the result (78) if we can show that
3
Mm =
i=1
m ∂ IN A˜ is − 2IN S = 0. ∂ξ i
(80)
When the matrices Ai are constant, the symmetrization matrix B is also constant, so the matrix
S = B−1
3 i=1
∂ B ˜i A = 0. ∂ξ i s
(81)
Also when B is constant, 3
i=1
Now, 3
i=1
3 ∂ IN A˜ is ∂ IN A˜ i −1 = B B. ∂ξ i ∂ξ i i=1
(82)
3 3 ∂ IN A˜ i ∂ = IN Jain An . ∂ξ i ∂ξ i n=1 i=1
(83)
Using the fact that the An are constants, 3
i=1
3 3 ∂ IN Jain ∂ IN A˜ is −1 =B An ∂ξ i ∂ξ i n=1 i=1
(84)
Therefore, if the sums,
3 ∂ IN Jain i=1
∂ξ i
=0
n = 1, 2, 3,
(85)
i.e., the discrete metric identities hold, then Mm = 0.
4.1.1. Example As an example of the stabilization that the use of the skew-symmetric approximation provides, we approximate a constant solution to the symmetric wave equation,
⎡ ⎤ p
⎡
0
⎣u⎦ + ⎣ c v t 0
⎤⎡ ⎤ p
⎡
⎤⎡ ⎤
c
0
0
0
c
0
0⎦⎣u⎦ + ⎣0
0
0⎦⎣u⎦ = 0,
0
0
0
0
v
x
c
p
v
(86)
y
where c = 1 is the constant wave speed, on the mesh shown in Fig. 1. The boundaries are approximated as polynomials of degree five. We choose the initial and external boundary states so that the solution is the Gaussian plane wave
⎡ ⎤
⎡
1
⎤
⎢ ⎥ p ⎢ kx ⎥ (kx x2 +ky y2 −ct )2 ⎣u⎦ = ⎢ c ⎥e− w2 ⎢ ⎥ ⎣ ⎦ v ky
(87)
c
√ with kx = ky = 1/ 2 and w = 1.6. t ) as a function of time step number for the central Riemann solver, λ = 0, and N = 5. Fig. 2 shows the maximum residual (Q The original, non-skew-symmetric, DGSEM approximation is weakly unstable in that the maximum residual grows slowly with time, whereas the skew-symmetric approximation formulation stabilizes the approximation. Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
[m3Gsc;September 3, 2015;10:37] 11
Fig. 1. Grid with curved elements for testing.
Log10(Max. Residual)
0
1
2 Skew-Symmetric DGSEM 3
4 0
20000 Time Step
40000
Fig. 2. Residual growth in the DGSEM shows weak instability for the central Riemann solver whereas the skew-symmetric approximation is stable as a plane wave propagates through and then leaves the domain.
4.2. Constant state preservation Next we show that a constant state is preserved by the skew-symmetric approximation of a constant coefficient system if the metric identities hold. We start with the following result that is needed to show that the flux contributions at element interfaces cancel when the solution is constant. Lemma 1. The numerical flux is consistent if the mesh is conforming, i.e.
Q = F˜ i Q . F˜si,∗ Q, s
(88)
Proof. Recall that
ai · F˜si Q¯ = A˜ is Q¯ = B−1 A˜ i BQ¯ = B−1 J
3
¯ A j xˆ j BQ.
(89)
j=1
The numerical flux is
A˜ s L R 1 i L R = R . ,Q + F˜ i Q L − Q +λ F˜si,∗ Q F˜ Q Q s 2 s 2 i
(90)
Then for either the upwind or central numerical flux
Q = 1 F˜ i Q + F˜ i Q . F˜si,∗ Q, s 2 s
(91)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 12
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
Expanding out,
i,∗
F˜s
L i
J a
Q = B−1 Q,
R
+ J ai 2
·
3
A j xˆ j BQ.
(92)
j=1
But J ai = ∇ξ i is continuous across element boundaries if the mesh is conforming, because
∇ξ
i
1 = J
∂ X ∂ X × ∂ξ j ∂ξ k
i, j, k cyclic
(93)
contains only derivatives tangent to the surface. For conforming boundaries, these tangential derivatives match exactly and the result follows. Constant state preservation is then guaranteed by Theorem 6. For constant coefficient systems, if the discrete metric identities (79) hold, then the skew-symmetric approximation (45) preserves a constant state Q¯ = c, that is, the time derivative at each nodal point (i, j, k) in any element m, Qtm = 0. i, j,k
Proof. We write the approximation (45) as
JQ¯t , φ
N
where
1 1 1 Term2 + Term3 − Term4 + Term5 = Term1 2 2 2
−
¯ φ Term1 = IN (S)Q, Term2 =
3
∂ φ
F˜si ,
N
∂ξ i N 3 ∂ IN A˜ is Term3 = φ , Q¯ ∂ξ i i=1 N 3 ∂ IN A˜ is φ Term4 = , Q¯ ∂ξ i i=1 N 3 i,∗ T 1 i Term5 = φ dS F˜s i i=1
(94)
∂ E ,N
i=1
(95)
−1
and study each term individually. • Term1 When the coefficients are constant, the matrix B is constant so that
S = B−1
3 i=1
Therefore,
∂ B ˜i A = 0. ∂ξ i s
¯ φ Term1 = IN (S)Q, • Term2 Next, 3
∂ φ ∂ξ i
F˜si ,
i=1
Now,
= N
N
(96)
= 0.
3 ∂ E i ,N
i=1
(97)
3 i ∂ F˜ i T 1 i s ˜ ,φ . Fs φ dS − ∂ξ i −1 N i=1
F˜si = IN A˜ is B−1 c = IN B−1 A˜ i BB−1 c = IN B−1 A˜ i c, so 3 ∂ F˜si i=1
∂ξ i
=
3 i=1
3 ∂ c IN B−1 J ai · A j xˆ j ∂ξ i j=1
(98)
(99)
(100)
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
13
Then for constant coefficients, 3 i=1
3 3 3 i ∂ IN Jai ∂ −1 N −1 B I J a · A j xˆ j c = B · A j xˆ j c, i i ∂ξ ∂ξ j=1 i=1 j=1
(101)
and the term in the parentheses vanishes by the assumption that the metric identities are satisfied discretely. Therefore, only the boundary integral term in (98) remains and
Term2 = I1, where 3
I1 =
∂ E i ,N
i=1
(102)
i T 1 i dS . F˜s φ
(103)
−1
• Term3 From (84) and the assumption that the discrete metric identities hold, 3 ∂ IN A˜ is
∂ξ i
i=1
= 0.
(104)
Therefore,
Term3 = 0.
(105)
• Term4 If we sum by parts,
i T 1 i ¯ ∂ IN A˜ is φ dS − A˜ i φ , ∂Q ¯ ˜s φ , Q = . F s ∂ξ i ∂ξ i N −1 ∂ E i ,N
(106)
N
But by hypothesis, Q¯ = Bc so the volume term vanishes, leaving
Term4 = I1.
(107)
At this point, we have
JQ¯t , φ
N
+ Term5 − I1 = 0,
(108)
so the time derivative depends only on the contribution of the boundary terms:
Term5 − I1 =
3 i=1
By Lemma 1,
∂ E i ,N
F˜si,∗
T
T 1 i φ dS .
− F˜si
−1
¯ Q¯ = F˜si Q¯ F˜si,∗ Q,
(109)
(110)
so Term5 − I1 = 0, and
JQ¯t , φ
N
= 0.
(111)
is arbitrary, we can take it to be the Lagrange basis function φ = (ξ ) (η) (ζ )eˆn , where eˆn is the unit vector Finally, since φ i k k whose nth component is one and the rest are zero. The discrete orthogonality then implies that
Ji, j,k Q¯t i.e.,
Q¯t
i, j,k
i, j,k
wi w j wk = 0,
=0
(112)
(113)
so the constant state is preserved and the result follows since B is constant. 4.2.1. Example As an example of what happens when we do or do not satisfy the metric identities, we approximate a constant solution to the symmetric wave equation, (86), on the mesh shown in Fig. 1. A constant is a solution to the wave equation and that constant should remain constant for all time. For this example, we compute the metric terms by
J a1 = Yη xˆ − Xη yˆ ˆ J a2 = −Yξ xˆ + Xξ y.
(114)
= X xˆ + Y y, ˆ is the linear transfinite interpolation [12] of the element edges to the reference element. The mapping, X Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 14
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
2 N=3 2 4
Log10(Max. Error)
Log10(Max Error)
4 6
8
10
N=4
6
8
10 N=5-8
12
12 3
4
5
6
7
8
0
N
20000 Time Step
40000
Fig. 3. Maximum error for the constant state solved on the mesh of Fig. 1 (left) and time history of the maximum error for the initially constant state (right). For N ≥ 5 the discrete metric identities are satisfied.
∈ PN , where N is the order at which It was shown in [14] that the discrete metric identities are satisfied if and only if X the solution is approximated. It was noted there that with isoparametric approximation of the boundaries by polynomials the condition is satisfied, but if the boundaries are approximated by polynomials of higher order than the solution then it is not. Therefore, if the boundaries are represented at a fixed order then the constant state is preserved if and only if the solution is approximated at that order or higher. For the example, we approximate the curved boundaries of the domain in Fig. 1 by polynomials of degree five and set the external boundary states and initial conditions to be a constant. The metric identities are satisfied therefore only for N ≥ 5 and we would not expect the constant solution to satisfy the discrete equations for solution approximation orders less than five. Fig. 3 shows that the constant state is preserved when the discrete metric identities are satisfied, but otherwise not. Fig. 3 shows the maximum error as a function of approximation order for the constant state solution. We see large errors for N = 3 and N = 4, i.e., when the metric identities are not satisfied, but as soon as they are at N = 5, the errors drop almost eight orders of magnitude to the level of rounding. The time histories of the errors computed with the skew-symmetric approximation are also shown in Fig. 3. When the discrete metric identities are not satisfied (N = 3, 4) the constant state is not a solution of the discrete equations and the solution adjusts. For N ≥ 5 the metric identities are satisfied and the time derivatives remain near rounding error for all time. 5. Conservation Finally, we show that the addition of the non-conservative form of the equation to the DGSEM does not change the conservation properties of the approximation. Theorem 7. If the symmetric equation (13) is strongly conservative, i.e., S = 0, the approximation (45) is discretely conservative in that Nel
Jm Q¯tm dξ =
m=1E,N
∂ E i ,N boundary faces
1
F˜si,∗
−1
dSi .
(115)
Remark 5. For constant coefficient PDEs, S = 0 since it depends on derivatives of the covariant coefficient matrices and the derivatives of the symmetrization matrix. = eˆn , then each component of the system satisfies Proof. If we set φ
JQ¯t , eˆn
3
N
3 T 1 i ˆn 1 i ∂e ˜ + eˆn dS − Fs , 2 ∂ξ i N −1 ∂ E i ,N i=1 i=1 3 ∂ IN A˜ is eˆn 1 ∂ IN A˜ is ¯ ¯ eˆn + − , Q¯ = IN (S)Q, eˆn , Q i i N 2 ∂ξ ∂ξ N i=1
F˜si,∗
(116)
N
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
First, when S = 0,
¯ eˆn IN (S)Q,
N
15
= 0.
(117)
Next, it is immediate that 3
∂ eˆn ∂ξ i
F˜si ,
i=1
= 0.
(118)
N
The most important fact is that the inclusion of the “nonconservative” terms have no effect since eˆn is constant and A˜ is is symmetric so that
N i ∂ IN A˜ is ∂ I A˜ s eˆn ∂ IN A˜ is ∂ IN A˜ is − , Q¯ = eˆn , − = 0. eˆn , eˆn , Q¯ Q¯ Q¯ ∂ξ i ∂ξ i ∂ξ i ∂ξ i N
N
N
(119)
N
Therefore, conservation depends only on the surface integral of the fluxes, as required:
JQ¯t , eˆn
N
+
3
∂ E i ,N
i=1
F˜si,∗
T 1 i eˆn dS = 0. −1
(120)
This says that on each element, the solutions satisfy a discrete conservation property element-wise and that N
JQ¯
i jk
wi w j wk
(121)
i jk
is conserved. Finally, global conservation is ensured by the continuity of the numerical fluxes, which cancel leaving only the integral of the fluxes at physical boundaries. 5.1. Example As an example showing the conservation properties of the skew-symmetric approximation, we integrate the wave equation (86) with the initial condition
p(x, y, 0) = exp − ln 2
x2 +y2 2.32
(122)
u(x, y, 0) = v(x, y, 0) = 0 which produces an expanding p wave
p(x, y, t ) = −
∞ 0
e−ω /4b ωJ0 (rω) cos (ωt )dω, 2b 2
(123)
where J0 is the Bessel function of the first kind of order zero, and b = ln 2/w2 with w = 2.3. We compute the solution on the circular mesh shown in Fig. 4 on the left. A snapshot of p is shown in Fig. 4 on the right.
Fig. 4. Mesh and contours of the expanding p wave to test conservation of the skew-symmetric approximation.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
ARTICLE IN PRESS
JID: AMC 16
[m3Gsc;September 3, 2015;10:37]
D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
1.0
Normalized Total
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
Time Fig. 5. Normalized total p for the expanding wave as a function of time showing conservation until the wave is lost through the outer boundary
Fig. 5 shows the normalized total amount of p over the entire domain. That is, we plot Nel N
Ji,mj p¯ m i, j wi w j
m=1 i, j=0
normalized to its initial value as a function of time. We see that the total remains constant until it is lost through the outer mesh boundary, as expected. Fig. 5 supports the result of Theorem 7 that the skew-symmetric approximation is conservative. 6. Conclusions In this work, we analyze the stability of a skew symmetric DGSEM for constant coefficient advection problems on curved domains. The continuous problem is discretized with the novel DGSEM introduced in [16]. We show the importance of the so-called metric identities for the skew-symmetric DGSEM. Satisfying the metric identities is important to guarantee the socalled ‘free-stream preservation’ in a standard DGSEM. We prove that the same condition guarantees free-stream preservation for the new scheme and furthermore that the skew-symmetric formulation being the average of the conservative and the nonconservative form of the equation is discretely conservative when the metric identities hold and a conforming mesh is used. Thus, in conclusion, the new DGSEM formulation satisfies the same properties as the original: free-stream preservation and element-wise conservation but is provably stable on curved elements. Additionally, positive properties of a high order DGSEM, such as the good parallel efficiency, are still valid since the additional skew-symmetric terms depend on local element data only. References [1] G. Cohen, X. Ferrieres, S. Pernet, A spatial high-order hexahedral discontinuous Galerkin method to solve Maxwell’s equations in time domain, J. Comput. Phys. 217 (2) (2006) 340–363. [2] D.A. Kopriva, S.L. Woodruff, M.Y. Hussaini, Discontinuous spectral element approximation of Maxwell’s Equations, in: B. Cockburn, G. Karniadakis, C.-W. Shu (Eds.), Proceedings of the International Symposium on Discontinuous Galerkin Methods, Springer-Verlag, New York, 2000, pp. 355–361. [3] D.A. Kopriva, S.L. Woodruff, M.Y. Hussaini, Computation of electromagnetic scattering with a non-conforming discontinuous spectral element method, Int. J. Numer. Methods Eng. 53 (1) (2002) 105–122. [4] S. Deng, Numerical simulation of optical coupling and light propagation in coupled optical resonators with size disorder, Appl. Numer. Math. 57 (5–7) (2007) 475–485, doi:10.1016/j.apnum.2006.07.001. [5] S. Deng, W. Cai, V. Astratov, Numerical study of light propagation via whispering gallery modes in microcylinder coupled resonator optical waveguides, Opt. Express 12 (26) (2004) 6468–6480. [6] P. Rasetarinera, D. Kopriva, M. Hussaini, Discontinuous spectral element solution of acoustic radiation from thin airfoils, AIAA J. 39 (11) (2001) 2070–2075. [7] D. Stanescu, F. Farassat, M. Hussaini, Aircraft engine noise scattering – parallel discontinuous Galerkin spectral element method, Paper 2002-0800, AIAA, 2002. [8] D. Stanescu, J. Xu, F. Farassat, M. Hussaini, Computation of engine noise propagation and scattering off an aircraft, Aeroacoustics 1 (4) (2002) 403–420. [9] N. Castel, G. Cohen, M. Durufle, Application of discontinuous Galerkin spectral method on hexahedral elements for aeroacoustic, J. Comput. Acoust. 17 (2) (2009) 175–196. [10] G. Gassner, D.A. Kopriva, A comparison of the dispersion and dissipation errors of Gauss and Gauss–Lobatto discontinuous Galerkin spectral element methods, SIAM J. Sci. Comput. 33 (5) (2010) 2560–2579. [11] C. Altmann, A.D. Beck, F. Hindenlang, M. Staudenmaier, G.J. Gassner, C.-D. Munz, An efficient high performance parallelization of a discontinuous Galerkin spectral element method, in: R. Keller, D. Kramer, J.-P. Weiss (Eds.), Facing the Multicore-Challenge III, Lecture Notes in Computer Science, vol. 7686, Springer, Berlin, Heidelberg, 2013, pp. 37–47.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047
JID: AMC
ARTICLE IN PRESS D.A. Kopriva, G.J. Gassner / Applied Mathematics and Computation 000 (2015) 1–17
[m3Gsc;September 3, 2015;10:37] 17
[12] D.A. Kopriva, Implementing Spectral Methods for Partial Differential Equations, Scientific Computation, Springer, 2009. [13] F. Hindenlang, G.J. Gassner, C. Altmann, A. Beck, M. Staudenmaier, C.-D. Munz, Explicit discontinuous Galerkin methods for unsteady problems, Comput. Fluids 61 (2012) 86–93, doi:10.1016/j.compfluid.2012.03.006. [14] D.A. Kopriva, Metric identities and the discontinuous spectral element method on curvilinear meshes, J. Sci. Comput. 26 (3) (2006) 301–327. [15] J. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer, 2008. [16] D.A. Kopriva, G.J. Gassner, An energy stable discontinuous Galerkin spectral element discretization for variable coefficient advection problems, SIAM J. Sci. Comput. 36 (4) (2014) A2076–A2099. [17] J. Nordström, Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation, J. Sci. Comput. 29 (3) (2006) 375–404, doi:10.1007/s10915-005-9013-4. [18] J. Kozdon, E. Dunham, J. Nordström, Simulation of dynamic earthquake ruptures in complex geometries using high-order finite difference methods, J. Sci. Comput. 55 (1) (2013) 92–124, doi:10.1007/s10915-012-9624-5. [19] M. Deville, P. Fischer, E. Mund, High Order Methods for Incompressible Fluid Flow, Cambridge University Press, 2002. [20] W. Gordon, C. Hall, Construction of curvilinear co-ordinate systems and their applications to mesh generation, Int. J. Numer. Methods Eng. 7 (1973) 461–477. [21] K. Korczak, A. Patera, An isoparametric spectral element method for solution of the Navier–Stokes equations in complex-geometry, J. Comput. Phys. 62 (2) (1986) 361–382. [22] H.-O. Kreiss, J. Lorenz, Initial-Boundary Value Problems and the Navier–Stokes Equations, Pure and Applied Mathematics, 136, Academic Press, San Diego, CA, USA, 1989. [23] D.A. Kopriva, G. Gassner, On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods, J. Sci. Comput. 44 (2) (2010) 136–155. [24] G. Gassner, A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods, SIAM J. Sci. Comput. 35 (3) (2013) A1233–A1253. [25] C. Canuto, A. Quarteroni, Approximation results for orthogonal polynomials in Sobolev spaces, Math. Comput. 38 (157) (1982) 67–86, doi:10.2307/2007465. [26] J. Thompson, Z. Warsi, C. Mastin, Numerical Grid Generation, Elsevier Science Publishing, Amsterdam, 1985.
Please cite this article as: D.A. Kopriva, G.J. Gassner, Geometry effects in nodal discontinuous Galerkin methods on curved elements that are provably stable, Applied Mathematics and Computation (2015), http://dx.doi.org/10.1016/j.amc.2015.08.047