Space-Time Discontinuous Galerkin Discretizations for Linear First-Order Hyperbolic Evolution Systems

1 Introduction

Space-time methods for time-dependent PDEs discretize the full problem in the space-time cylinder, and then the corresponding large algebraic system is also solved for the full problem. This is in contrast to the method of lines or Rothe’s method, which first use a discretization either in space or in time and then apply standard techniques for the other variable. Our methods are based on treating space and time simultaneously in a variational manner. Depending on the choice of the ansatz and the test spaces, the methods become either explicit or implicit. Explicit methods are computationally efficient but suffer from severe limitations for the time step size, where the length of the time edge of the space-time elements is restricted by the smallest local resolution scale in space. To circumvent these restrictions, we focus on implicit methods.

A fully implicit space-time approach allows for flexible adaptive discretizations which combine adaptivity in space with local time stepping. A further motivation for developing space-time methods is the design of modern computer facilities with an enormous number of processor cores, where the parallel realization of conventional methods becomes inefficient. Since these machines allow a fully implicit space-time approach, new parallel solution techniques are required to solve the huge linear systems, particularly for time-dependent applications in three spatial dimensions.

In recent years, discontinuous Galerkin (DG) methods in space have become very popular, see, e.g., [18] for time-dependent first-order systems, where this discretization is coupled with explicit time integration. An application of this method to acoustic and elastic waves is considered in [7] combined with an adaptive space-time hp-strategy. Here, we extend these spatial DG discretization by a Petrov–Galerkin method in time with continuous ansatz space and discontinuous test space (see, e.g., [3]). A space-time method for elastic waves with a second-order formulation in space and implicit discontinuous Galerkin time discretization is considered in [22].

An alternative discontinuous Petrov–Galerkin (DPG) approach is proposed by L. Demkowicz (see [6] for an overview and [9] for space-time applications) for general linear first-order systems, where weak approximations are constructed by introducing skeleton variables. The application of this technique to the time-harmonic case is analyzed in [32]. For acoustic and elastic waves, the hybridization in space (applied to the second-order formulation) is presented in [26], and a hybrid space-time discontinuous Galerkin method is proposed in [29]. Both methods are implicit in every time slab, and only Dirichlet traces are used for the hybrid coupling. Space-time (Trefftz) discontinuous Galerkin methods for wave problems are analyzed in [8, 23].

Error estimation for linear wave equations (considered as second-order equations) is studied in [3, 27, 17] and for more general hyperbolic systems in [21]. Simple residual error indications are not sufficient for wave problems since, in the hyperbolic case, the error is transported and thus not correlated to large local residuals. Reliable error control requires the adjoint problem, as it is introduced for goal-oriented techniques in [3], to be solved. This technique requires a variational approach, since this allows for an error representation with respect to a given linear error functional.

In principle, all parallel solution methods in space apply also to implicit time integration schemes. Parallel strategies in time are studied extensively on the basis of the ‘parareal’ idea [24, 2, 16]. A general overview over the most popular algorithms and software packages is given in [13]. Methods such as MGRID [12] and PFASST [10] were developed under the aspect that they can be easily incorporated into existing time sequential code. In addition, solution concepts specially adapted to the full space-time problem were proposed. E.g., the wavefront method extends a spatial domain decomposition into time slices, see [14] for an application to the one-dimensional wave equation. In [15] a space-time multigrid method for parabolic problems is analyzed. A multigrid method for higher order discontinuous Galerkin discretizations of advection problems is proposed in [28].

In this paper we present a fully implicit and parallel adaptive space-time discontinuous Galerkin discretization for linear first-order hyperbolic problems. The paper is structured as follows. In Section 2 we introduce a setting for linear hyperbolic operators by reference to applications in the field of linear transport and acoustic and electro-magnetic waves, and we establish the well-posedness of the space-time variational problem based on a technique developed in [31]. In Section 3, following the setting established in [20], we consider a semi-discrete discontinuous Galerkin discretization in spatial direction with upwind flux. On this basis we define an implicit Petrov–Galerkin space-time discretization in Section 4, and we prove well-posedness of the discrete method and convergence on tensor product space-time meshes. In Section 5 we propose a goal-oriented space-time error indicator based on the explicit computation of the dual solution. In Section 6 a multilevel preconditioner with semi-coarsening first in time and then in space is defined. Within the parallel finite element software system M++ [30] the adaptive method and the multilevel solution method are realized in Section 7. Moreover, the efficiency of the full scheme is demonstrated for two models, the linear transport equation and Maxwell’s equations in 2D.

2 A Space-Time Setting for Linear Hyperbolic Operators

Let Ω⊆ℝD be a bounded Lipschitz domain, and let H⊆L2⁢(Ω)J be a Hilbert space with weighted inner product (𝐯,𝐰)H=(M⁢𝐯,𝐰)0,Ω, where M∈L∞⁢(Ω)J×J is uniformly positive and symmetric. We consider a linear operator A:D⁢(A)→H with domain D⁢(A)⊂H. For given initial function 𝐮0∈D⁢(A), final time T>0 and right-hand side 𝐟∈L2⁢(0,T;H), we study the evolution equation

We specialize A to the case of linear balance laws determined by a flux function𝐅⁢(𝐯)=[B1⁢𝐯,…,BD⁢𝐯] with symmetric matrices Bj∈L∞⁢(Ω)J×J such that

Since the matrices Bd are symmetric, any linear combination n1⁢B1+…+nD⁢BD for 𝐧=(n1,…,nD)⊤∈ℝD is diagonalizable with real eigenvalues, so that (2.1) is a linear hyperbolic system [11, Section 7.3]. We discuss the following examples.

3 A Semi-Discrete Discontinuous Galerkin Discretization in Space

In this section we consider the semi-discrete evolution equation

in a finite dimensional subspace Hh⊂H associated to the mesh size h of the underlying mesh defined below. The discrete operator Ah∈ℒ⁢(Hh,Hh) will be constructed from a discontinuous Galerkin discretization. The discrete mass operator Mh∈ℒ⁢(Hh,Hh) is the Galerkin approximation of M defined by

Note that the discrete mass operator Mh is represented by a block diagonal positive definite matrix.

We assume that Ω is a bounded polyhedral Lipschitz domain decomposed into a finite number of open elements K⊂Ω such that Ω¯=⋃K∈𝒦K¯, where 𝒦 is the set of elements in space. Let ℱK be the set of faces of K, and for inner faces f∈ℱK let Kf be the neighboring cell such that f=∂⁡K∩∂⁡Kf, and let 𝐧K be the outer unit normal vector on ∂⁡K. The outer unit normal vector field on ∂⁡Ω is denoted by 𝐧.

Integration by parts of A⁢𝐯=div⁡𝐅⁢(𝐯) gives for smooth ansatz functions 𝐯 and smooth test functions ϕK,

This formulation is now the basis for the discretization. We select polynomial degrees pK, and we define the local spaces Hh,K=ℙpK⁢(K)J and the global discontinuous Galerkin space

For 𝐯h∈Hh we define 𝐯h,K=𝐯h|K∈Hh,K for the restriction to K.

We then define the discrete linear operator Ah∈ℒ⁢(Hh,Hh) for 𝐯h∈Hh and ϕh,K∈Hh,K by

where 𝐧K⋅𝐅Knum⁢(𝐯h) is the upwind flux obtained from local solutions of Riemann problems, see [20, Section 2]. Again using integration by parts, we obtain

On inner faces f=∂⁡K∩∂⁡Kf it is a consistency requirement that the difference 𝐧K⋅(𝐅Knum⁢(𝐯h)-𝐅⁢(𝐯h,K)) only depends on [𝐯h]K,f=𝐯h,Kf-𝐯h,K, and that 𝐧K⋅(𝐅Knum⁢(𝐯)-𝐅⁢(𝐯))=0 on all faces f∈ℱK for 𝐯∈D⁢(A). In particular, this yields

and

The upwind flux guarantees that the discrete operator is non-negative, i.e., (Ah⁢𝐯h,𝐯h)0,Ω≥0 for 𝐯h∈Hh. For the examples in Section 2, the numerical upwind flux in homogeneous media is given as follows (see [20] for the explicit solution of Riemann problems in heterogeneous media).

4 A Petrov–Galerkin Space-Time Discretization

Let Q¯=⋃R∈ℛR¯ be a decomposition of the space-time cylinder into space-time cells R=K×I with K⊂Ω and I=(t-,t+)⊂(0,T); ℛ denotes the set of space-time cells. For every R we choose local ansatz and test spaces Vh,R,Wh,R⊂L2⁢(R)J with Wh,R⊂∂t⁡Vh,R, and we define the global ansatz and test space

By construction, functions in Wh are discontinuous in space and time, and functions in Vh are continuous in time, i.e., 𝐯h⁢(𝐱,⋅) is continuous on [0,T] for a.a. 𝐱∈Ω.

In addition we aim for dim⁡(Vh)=dim⁡(Wh), which restricts the choice of Vh,R. In the most simple case this can be achieved for a tensor product space-time discretization with a fixed mesh 𝒦 in space and a time series 0=t0<t1<⋯<tN=T, i.e., ℛ=⋃K∈𝒦⋃n=1NK×(tn-1,tn). Then, we can select a discrete space Hh with Hh,K=ℙp⁢(K)J independently of t, and in every time slice we define Wh,R=Hh,K constant in time on R=K×(tn-1,tn). This yields in this case piecewise linear approximations in time,

In the general case, we select locally in space and time polynomial degrees pR and qR in R, and we set for the local test space Wh,R=(ℙpR⁢(K)⊗ℙqR-1)J. Then we define for R∈ℛ,

This yields 𝐯h,R⁢(𝐱,⋅)∈ℙqRJ for 𝐯h,R∈Vh,R and (𝐱,⋅)∈R.

The discontinuous Galerkin operator in space is extended to the space-time setting defining Ah⁢𝐯h∈Wh by

for 𝐯h∈Vh and 𝐰h∈Wh. We define the discrete space-time operator Lh∈ℒ⁢(Vh,Wh) and the corresponding discrete bilinear form bh(⋅,⋅)=(Lh⋅,⋅)0,Q by (Lh⁢𝐯h,𝐰h)0,Q=(Mh⁢∂t⁡𝐯h+Ah⁢𝐯h,𝐰h)0,Q.

In order to show that a solution to our Petrov–Galerkin scheme exists, we check the inf-sup stability of the discrete bilinear form bh⁢(⋅,⋅) with respect to the discrete norm

By construction, bh⁢(⋅,⋅) is bounded in Vh×Wh, i.e.,

For the verification of the inf-sup stability, we introduce the L2-projection Πh:W→Wh which is defined by (Πh⁢𝐯,𝐰h)0,Q=(𝐯,𝐰h)0,Q for 𝐰h∈Wh. Then, by construction, Πh⁢Ah=Ah and Πh⁢Lh=Lh. Moreover, we define the non-negative weight function in time dT⁢(t)=T-t, and we observe

Referring to Theorem 2.2, Lemma 4.1 shows the existence of a unique discrete Petrov–Galerkin solution (provided that the assumptions in Lemma 4.1 are satisfied).

The convergence will be analyzed with respect to the discrete norm ∥⋅∥Vh. For 𝐯∈V the consistency of the numerical flux in (4.1) yields (Ah⁢𝐯,𝐰h)0,Q=(div⁡𝐅⁢(𝐯),𝐰h)0,Q so that Ah⁢𝐯=Πh⁢div⁡𝐅⁢(𝐯). This shows that Ah and thus also ∥⋅∥Vh can be evaluated in V+Vh and that bh⁢(⋅,⋅) is continuous with respect to this extension.