Space-time discontinuous Petrov-Galerkin

methods for linear wave equations in heter-

ogeneous media

Johannes Ernesti, Christian Wieners

CRC Preprint 2018/13, August 2018

KARLSRUHE INSTITUTE OF TECHNOLOGY

KIT – The Research University in the Helmholtz Association www.kit.edu

Participating universities

Funded by

ISSN 2365-662X
 SPACE-TIME DISCONTINUOUS PETROV-GALERKIN METHODS
FOR LINEAR WAVE EQUATIONS IN HETEROGENEOUS MEDIA

Johannes Ernesti and Christian Wieners 1

Abstract. We establish an abstract space-time DPG framework for the approximation of linear waves
in heterogeneous media. The estimates are based on a suitable variational setting in the energy space.
The analysis combines the approaches for acoustic waves in Gopalakrishnan / Sepulveda (A space-
time DPG method for acoustic waves, arXiv 2017) and in Ernesti / Wieners (RICCAM proceedings,
submitted 2017) and is based on the abstract definition of traces on the skeleton of the time-space sub-
structuring. The method is evaluated by large-scale parallel computations motivated from applications
in seismic imaging, where the computational domain can be restricted substantially to a subset of the
full space-time cylinder.

1991 Mathematics Subject Classification. 65N30.

August 14, 2018.

1. Introduction

Space-time finite elements aim for a unified analysis of discretization and solution methods in space and time.
In particular they allow for an efficient combined error control and for scaling of the solution scheme to the next
generation of massively parallel computers.
The discontinuous Petrov-Galerkin method (DPG) is a well-suited finite element class for space-time applications
which provides robust a-priori estimates, reliable error control, and the efficient hybridization to a symmetric
positive definite Schur complement system. This is attractive for hyperbolic systems and allows to transfer
features of discretizations for elliptic problems to wave-type equations. The long-term goal is, as it is discussed
in [3] for the transport equation, to establish optimality of the solution process and of adaptive schemes. For a
general discussion on the DPG technology we refer to [7].
First results of space-time DPG methods are established in [8] for the Schrödinger equations and in [13, 15] for
acoustic waves. Here, we show that the analysis transfers to general wave equations in heterogeneous media and
provides robust estimated in the energy norm. Therefore, we recall in Lem. 3 and Lem. 4 the abstract DPG
analysis based on the technique introduced in [15] which avoids explicit traces. Then, following the arguments
in [3] we show that a test space exists which guarantees discrete inf-sup stability for general wave equations,
and we extend the analysis for the simplified DPG method with nonconforming traces as in [13] to this more
general setting. Finally, we apply a Strang-type argument to estimate the consistency error of the DPG method
due to inexact quadrature in heterogeneous media.

Keywords and phrases: space-time methods, discontinuous Petrov-Galerkin finite elements, linear hyperbolic systems, heteroge-
nous media
1 Institut für Angewandte und Numerische Mathematik, KIT, Karlsruhe, Germany. Email: [email protected]
2

The analysis is complemented by numerical results for wave propagation in heterogeneous media. Here we
discuss an application scenario motivated from seismic measurements, were the wave signal is initiated by a
point source and the results are only measured at selected points. In this application class the finite propagation
speed of wave solutions results into an a priori information about the region of interest within the space-time
cylinder and which allows to truncate the computational domain substantially.

2. Linear hyperbolic systems

A semigroup framework. We consider the evolution equation

M ∂t u + Au = f in (0, T ) × Ω0 ⊂ R × Rd

subject to homogeneous initial conditions u(0) = 0 in Ω0 , where Ω0 is a Lipschitz domain, f (t) ∈ L2 (Ω0 ; Rm ) is
a source function, and with
a) a symmetric positive definite operator M in L2 (Ω0 ; Rm ) represented by M ∈ L∞ (Ω0 ; Rm×m
sym );
b) a hyperbolic operator A with domain D(A) ⊂ L2 (Ω0 ; Rm ) such that

(Av, z)0,Ω0 = −(v, Az)0,Ω0 , v, z ∈ D(A) (1)

and such that M + A is surjective, i.e.,

(M + A)(D(A)) = L2 (Ω0 ; Rm ) . (2)

Then, M −1 A generates a semigroup in L2 (Ω0 ; Rm ) and for f ∈ C0 (0, T ); D(A) the solution is given by

Z t
exp (t − s)M −1 A M −1 f (s) ds .


u(t) = (3)
0

The solution belongs to the Banach space

V = v ∈ C1 [0, T ]; L2 (Ω0 ; Rm ) ∩ C0 [0, T ]; D(A) : v(0) = 0 ,



see, e.g., [19, Thm. 12.22]. In particular, we obtain for the range of the space-time operator L = M ∂t + A

C1c (0, T ) × Ω0 ; Rm ⊂ C0 (0, T ); D(A) ⊂ L(V) ⊂ L2 (0, T ) × Ω0 ; Rm )

(4)


so that L(V) is dense in L2 (0, T ) × Ω0 ; Rm , cf. [19, Thm. 12.16].
Linear wave equations. Our basic example is the acoustic wave equation for velocity and pressure with
(v, p) ∈ D(A) = H(div, Ω0 ) × H10 (Ω0 ), M (v, p) = (ρv, κ−1 p), and A(v, p) = (∇p, ∇ · v).
Note that in this example the definition of the domain D(A) includes homogeneous Dirichlet boundary conditions
for the pressure p; more general boundary conditions can be included into the domain of the operator A, provided
the conditions (1) and (2) are satisfied (see also [16, Sect. 2.2] for various examples).
This framework also applies to linear elastic waves described by (v, σ) ∈ D(A) = H10 (Ω0 ; Rd ) × H(div, Ω0 ; Rsymd×d
),
−1
M (v, σ) = (ρv, C σ), and A(v, σ) = (div σ, ε(v)), and to electro-magnetic waves described by Maxwell’s
equations with (E, H) ∈ H0 (curl, Ω0 ) × H1 (curl, Ω0 ), M (E, H) = (εE, µH), and A(E, H) = (−∇ × H, ∇ × E).
Note that in all cases 12 M u(t), u(t) Ω0 is the free energy, i.e., for elastic waves the kinetic and potential energy,

and for the Maxwell case the electro-magnetic energy.

 3

A main property of the linear wave equation is the finite speed of propagation cmax > 0, which allows – in case
of local support of the source function f – to restrict the computation of u ∈ V solving Lu = f to the cone

C+ (supp f ) = (t, x) ∈ (0, T ) × Ω0 : |x − x0 | ≤ cmax (t − t0 ) for all (t0 , x0 ) ∈ supp f , (5)

i.e., supp u ∈ C+ (supp f ), cf. [14, Chap. 7.2.4]. The maximal wave speed can be determined by the equivalent
formulation as symmetric Friedrichs system, i.e., by the representation of the linear operator in the form
m×m
P
Av = Aj ∂j v with symmetric matrices Aj ∈ Rsym . Then, the maximal speed of propagation in heterogeneous
media is given by

d
w > An w X
cmax = kck∞,Ω0 , c(x) = max max , An = n j Aj . (6)
|n|=1 |w|=1 w > M (x)w
j=1

p
E.g., in the acoustic case we have c(x) = κ(x)/ρ(x).

x
x (t
− t 0)
= c ma
|x − x 0| P
|xP − x0 |
P0
x0
tP − t0

C+ {P0 }
|x − x
0| = cm
ax (t − t0 )
t
t0 T


Figure 1. The grey area depicts the cone of dependence C+ {P0 } ⊂ (0, T ) × Ω0 in 1D for
a single point source P0 . Due to the limit wave speed, S information originating

from P0 only
affects this space-time region, resulting in C+ (supp f ) = P0 ∈supp f C+ {P0 } for a right-hand
side f .

The adjoint equation. Let A∗ be the adjoint operator of A with domain D(A∗ ), and let L∗ be the adjoint
operator for the wave equation backward in time defined in

V ∗ = z ∈ C1 [0, T ]; L2 (Ω0 ; Rm ) ∩ C0 [0, T ]; D(A∗ ) : z(T ) = 0 .



Then, L∗ (V ∗ ) ⊂ L2 (Ω; Rm ) is dense, we have L∗ v = −Lv for v ∈ Cc1 ((0, T ) × Ω0 ; Rm ), and

(Lv, z)0,(0,T )×Ω0 = (v, L∗ z)0,(0,T )×Ω0 , v ∈ V , z ∈ V∗ . (7)

The corresponding backward cone for a domain of interest ω ⊂ [0, T ] × Ω0 is given by


C− (ω) = (t, x) ∈ (0, T ) × Ω0 : |x − x0 | ≤ cmax (t0 − t) for all (t0 , x0 ) ∈ ω . (8)

In the following, we consider applications where the source f has local support, and where the solution is
evaluated only in the domain of interest ω, so that the solution process can be restricted to C+ (supp f ) ∩ C− (ω).
4

Subsets ofSthe space-time cylinder. We consider Ω ⊂ (0, T ) × Ω0 ⊂ R × Rd combining time slices in the
form Ω = [tn−1 , tn ] × Ωn with open subsets Ωn ⊂ Ω0 for n = 1, . . . , N and a decomposition of the time
interval 0 = t0 < t1 < · · · < tN = T . In every time slice we select D(A; Ωn ) ⊂ D(A)
such that the conditions

(1) and (2) are satisfied in L2 (Ωn ; Rm
). This defines Vn = C1 [tn−1 , tn ]; L2 (Ωn ; Rm ) ∩ C0 [tn−1 , tn ]; D(A; Ωn ) .
Using vn ∈ C0 [tn−1 , tn ]; L2 (Ωn ; Rm ) for vn ∈ Vn we define

V(Ω) = v ∈ L2 (Ω; Rm ) : vn ∈ Vn , v(0) = 0 , vn (tn ) = vn+1 (tn ) on Ωn ∩ Ωn+1 and vn+1 (tn ) = 0 on Ωn+1\ Ωn ,


V ∗ (Ω) = z ∈ L2 (Ω; Rm ) : zn ∈ Vn , z(T ) = 0 , zn (tn ) = zn+1 (tn ) on Ωn ∩ Ωn+1 and zn (tn ) = 0 on Ωn \ Ωn+1


with vn = v|(tn−1 ,tn ) . By construction, we have

(Lv, z)0,Ω = (v, L∗ z)0,Ω , v ∈ V(Ω) , z ∈ V ∗ (Ω) , (9)

∗ ∗

and the ranges L V(Ω) , L V (Ω) are dense in the energy space W = L2 (Ω; Rm ). The corresponding energy
p
norm in space and time is given by kwkW = (M w, w)0,Ω for w ∈ W .
Lemma 1. We have for u ∈ V(Ω) and f = Lu

T
kukW ≤ √ kM −1 f kW .
2

This is a Poincaré type estimate since it relies on the initial condition u(0) = 0.

Proof. The estimate relies on the representation (3) of the solution in every

slice [tn−1 , tn ] × Ωn . In every slice
define Wn = L2 (Ωn ; Rm ) and we use the energy inner product vn , wn Wn = (M vn , wn )0,Ωn for vn , wn ∈ Wn .
Since the operator M −1 A is skew-adjoint in Wn , i.e., (M −1 Avn , vn )Wn = (Avn , vn )0,Ωn = 0 for vn ∈ D(A; Ωn ),
the spectrum is contained in iR which yields k exp(sM −1 A)vn kWn ≤ kvn kWn for all s ∈ R. Now, inserting
Z t
u(t) = exp (t − tn−1 )M −1 A u(tn−1 ) + exp (t − s)M −1 A M −1 f (s) ds ,



t ∈ (tn−1 , tn ] ,
tn−1

yields recursively for t ∈ (tn−1 , tn ]

Z t
ku(t)kWn ≤ ku(tn−1 )kWn + kM −1 f (s)kWn ds
tn−1
!1/2
Z t
−1
f (s)k2Wk
p
≤ ku(tn−1 )kWn + t − tn−1 kM ds
tn−1

n−1
!1/2 !1/2
X Z tk Z t
−1 −1
f (s)k2Wk f (s)k2Wn
p p
≤ tk − tk−1 kM ds + t − tn−1 kM ds
k=1 tk−1 tn−1

n−1
!1/2
√ X Z tk
−1
Z t
−1
≤ t kM f (s)k2Wk ds + kM f (s)k2Wn ds .
k=1 tk−1 tn−1

Together, this yields

N Z N Z
!
tn Z T tk
X X
−1 1 2
kuk2W = ku(t)k2Wn dt ≤ t kM f (s)k2Wk ds dt = T kM −1 f k2W .
n=1 tn−1 0 tk−1 2 
k=1
 5

A variational setting. We extend the operator L in V(Ω) ⊂ W to a suitable Hilbert space defining

H(L, Ω) = {v ∈ W : Lv ∈ W }
v ∈ W : w ∈ W exists such that (w, z)0,Ω = (v, L∗ z)0,Ω for all z ∈ C1c (Ω; Rm ) .

=
p p
We use the weighted graph norm kvkH(L,Ω) = kvk2W + kM −1 Lvk2W = (M v, v)0,Ω + (M −1 Lv, Lv)0,Ω . The
closure of C1c (Ω, Rm ) with respect to k · kH(L,Ω) is denoted by H0 (L, Ω).
In the hyperbolic case we have H(L, Ω) = H(L∗ , Ω). Nevertheless, since L is associated to the forward problem
and L∗ to the backward problem, we will need different subspaces in the following arguments.
Since C1 (Ω, Rm ) is dense in H(L, Ω), unique extensions L ∈ L(H(L, Ω), W ) and L∗ ∈ L(H(L∗ , Ω), W ) exist.
Let V be the closure of V(Ω) in H(L, Ω), and let V ∗ be the closure of V ∗ (Ω) in H(L∗ , Ω). Then, Lem. 1 yields

kvkW ≤ CL kM −1 LvkW , v∈V (10)

√
with CL = T / 2, and, correspondingly, kzkW ≤ CL kM −1 L∗ zkW for z ∈ V ∗ for the adjoint problem backwards
in time. Moreover, (10) implies that L(V ) ⊂ W is closed, and since L(V) ⊂ W is dense, we obtain L(V ) = W ,
so that L ∈ L(V, W ) is a bijection. Furthermore, (9) extends to

(Lv, z)0,Ω = (v, L∗ z)0,Ω , v ∈V , z ∈V∗. (11)

The solution in the restricted space-time domain Ω ⊂ (0, T ) × Ω0 is now compared with the solution in the
full space-time cylinder (0, T ) × Ω0 . Therefore, let V0 and V0∗
be the closures of V = V((0, T ) × Ω0 ) and
V ∗ = V ∗ ((0, T ) × Ω0 ). Furthermore, set W0 = L2 (0, T ) × Ω0 ; Rm .
Lemma 2. Let f0 ∈ W0 be a source function and let u0 ∈ V0 be the unique solution of Lu0 = f0 in (0, T ) × Ω0 .
Assume that Ω ⊃ C+ (supp f0 ) ∩ C− (ω) for a domain of interest ω ⊂ [0, T ] × Ω0 . Then, we have for the unique
solution u ∈ V of Lu = f with restricted source function f = f0 |Ω ∈ W the identical values in ω, i.e.,

u|ω = u0 |ω .

Proof. Let z0 ∈ V0∗ be the dual solution with respect to the goal functional (u − u0 )|ω defined by

(L∗ z0 , w0 )0,(0,T )×Ω = (u − u0 , w0 )0,ω , w0 ∈ W0 .

Then we observe for the dual solution supp z0 ⊂ C− (ω).

Let χ ∈ C1c (0, T ) × Ω0 be a function with supp χ ⊂ Ω and χ ≡ 1 in C = C+ (supp f0 ) ∩ C− (ω). This extends
u ∈ V to χu ∈ V0 with u(t, x) = 0 for (t, x) ∈ (0, T ) × Ω0 \ Ω. Using supp u ∈ C+ (supp f0 ) and (11) we obtain

(L∗ z0 , χu)0,Ω = (L∗ z0 , χu)0,(0,T )×Ω0 = (z0 , Lχu)0,(0,T )×Ω0 = (z0 , Lu)0,C = (z0 , Lu)0,Ω

and (L∗ z0 , χu0 )0,(0,T )×Ω0 = (z0 , Lu0 )0,(0,T )×Ω0 . This yields the assertion by

ku − u0 k20,ω


= u − u0 , χ(u − u0 ) 0,ω
= (L∗ z0 , χu)0,Ω − (L∗ z0 , χu0 )0,(0,T )×Ω0
= (z0 , Lu)0,Ω − (z0 , Lu0 )0,(0,T )×Ω0
= (z0 , f )0,Ω − (z0 , f0 )0,(0,T )×Ω0 = (z0 , f )0,Ω − (z0 , f0 )0,C = 0 .

6

Integration by parts defines the operators DΩ ∈ L H(L, Ω), H(L∗ , Ω)0 and DΩ
0
∈ L H(L∗ , Ω), H(L, Ω)0 with

hDΩ v, zi = (Lv, z)0,Ω − (v, L∗ z)0,Ω = hDΩ

0
z, vi , v ∈ H(L, Ω) , z ∈ H(L∗ , Ω) .

The kernel of DΩ is denoted by N (DΩ ) = y ∈ H(L, Ω) : hDΩ y, zi = 0 for all z ∈ H(L∗ , Ω) ⊃ H0 (L, Ω). In


fact, the boundary conditions H0 (L, Ω) and in V are characterized by duality.

Lemma 3. We have H0 (L, Ω) = N (DΩ ) and V = v ∈ H(L, Ω) : hDΩ v, zi = 0 for all z ∈ V ∗ = ⊥ DΩ 0
(V ∗ ).


The proof is based on [11, Lem. 2.4] and [5] and uses properties of polar sets, cf. [20, Sect. 4.5].

Proof. For a given functional ` ∈ C1c (Ω; Rm )⊥ = η ∈ H(L, Ω)0 : hη, vi = 0 for all v ∈ C1c (Ω; Rm ) we define


u` ∈ H(L, Q) solving
M −1 Lu` , Lv 0,Ω + M u` , v 0,Ω = h`, vi , v ∈ H(L, Q) .



(12)
Then, M u` , v 0,Ω = − M −1 Lu` , Lv 0,Ω for all v ∈ C1c (Ω; Rm ), which yields z` = M −1 Lu` ∈ H(L∗ , Q) and

L∗ z` = −M u` . This is inserted in (12), and we obtain for the adjoint operator DΩ0

0
z` , vi = hDΩ v, z` i = (Lv, z` )0,Ω − (v, L∗ z` )0,Ω = Lv, M −1 Lu`



hDΩ 0,Ω
+ v, M u` 0,Ω
= h`, vi , v ∈ H(L, Q) ,

0
z` = `. This proves C1c (Ω; Rm )⊥ ⊂ DΩ
0
H(L∗ , Ω) , and by duality we conclude the first assertion


i.e., DΩ

N (DΩ ) = ⊥ DΩ
0
H(L∗ , Ω) ⊂ ⊥ C1c (Ω; Rm )⊥ = H0 (L, Q) .

Since ⊥ DΩ
0
V
∗ is closed in H(L, Q) and V ⊂ ⊥ DΩ 0
V ∗ by (9), we have V ⊂ ⊥ DΩ
0
(V ∗ ). On the other hand, for



⊥ 0 ∗
v ∈ DΩ V ⊂ H(L, Q) set f = Lv and let u ∈ V be the unique solution of Lu = f . Together, we have by
construction u − v ∈ ⊥ DΩ0
(V ∗ ) and L(u − v) = 0, so that

0
z, u − v = L∗ z, u − v Q − z, L(u − v) Q = L∗ z, u − v Q , z ∈ V∗ .




0 = DΩ

Since L∗ (V ∗ ) ⊂ W is dense, we obtain u = v ∈ V . 

3. Space-time substructuring
S
Let Ωh = K ⊂ Ω be a decomposition into openQ convex subsets K with skeleton ∂Ωh satisfying Ω = Ωh ∪ ∂Ωh .
∗ ∗ 0 ∗ 0
We consider the broken spacep Z = H(L , Ω h ) = K H(L , K) and its dual space Z = H(L , Ωh ) , again using
−1 ∗
the weighted norm kzkZ = (M z, z)0,Ω + (M L z, L z)0,Ωh ∗

Integration by parts defines the operators DK ∈ L H(L, K), H(L∗ , K) and DΩh ∈ L H(L, Ωh ), Z 0 by

X X
hDΩh v, zi = (Lv, z)0,Ωh − (v, L∗ z)0,Ωh = (Lv, z)0,K − (v, L∗ z)0,K = hDK v|K , z|K i , v ∈ H(L, Ωh ) , z ∈ Z .
K K
Q
Lem. 3 yields H0 (L, K) = N (DK ) and thus H0 (L, Ωh
) = H0 (L, K) = N (DΩh ). This allows to identify traces
on the skeleton ∂Ωh with functionals DΩh H(L, Ωh ) ⊂ Z 0 . By construction, we have

= (M −1 Lv, M z)0,Ωh − (M v, M −1 L∗ z)0,Ωh  ≤ kvkH(L,Ωh ) kzkZ ,

 
|hDΩh v, zi| v ∈ H(L, Ωh ) , z ∈ Z ,

i.e., kDΩh vkZ 0 ≤ kvkH(L,Ωh ) .

 7

For given f ∈ W let u ∈ V be the unique solution of Lu = f . Then, we observe

(Lu, z)0,Ωh = (u, L∗ z)0,Ωh + hDΩh u, zi = (f, z)0,Ωh , z ∈Z.

This allows for a weak characterization in the energy space W and the abstract trace space V̂ = DΩh (V ) ⊂ Z 0 .
Therefore, we define the bilinear form in (W × Z 0 ) × Z by

b(w, v̂; z) = (v, L∗ z)0,Ωh + hv̂, zi , (w, v̂) ∈ W × Z 0 , z ∈ Z .

We now show that a unique weak solution (u, û) ∈ W × V̂ exists with of b(u, û; z) = (f, z)0,Ωh for all z ∈ Z.
In order to obtain improved optimal estimates, we use in W × Z 0 the norm k(w, v̂)kW ×Z 0 = max{kwkW , kv̂kZ 0 }.
Then, the bilinear form b(·; ·) is continuous with |b(w, v̂; z)| ≤ k(w, v̂)kW ×Z 0 kzkZ , and we obtain the following
estimate for the inf-sup stability for the ideal DPG method, cf. [7].
Lemma 4. The bilinear form b(·; ·) is injective on W × V̂ and inf-sup stable with β = (2CL2 + 2)−1/2 , i.e.,

b(w, v̂; z)
sup ≥ β k(w, v̂)kW ×Z 0 , (w, v̂) ∈ W × Ŷ . (13)
z∈Z kzkZ

Proof. Let (w, v̂) ∈ W × V̂ be in the kernel of b, i.e, b(w, v̂; z) = 0 for z ∈ Z. Then we have

0 = (w, L∗ z)0,Ωh + hv̂, zi = w, L∗ z 0,Ω , z ∈ C1c (Ωh , Rm ) ,

h

i.e., w ∈ H(L, Ωh ) and Lw = 0 in Ωh . This yields for all v ∈ V with DΩh v = v̂

hDΩh (w − v), zi = Lw, z 0,Ω − w, L∗ z 0,Ω − hv̂, zi = 0 − b(w, v̂; z) = 0 , z ∈ H(L∗ , Ωh ) ,

h h

i.e., w − v ∈ N (DΩh ) = H0 (L, Ωh ). Thus we have w ∈ v + H0 (L, Ωh ) ⊂ V , and together with Lw = 0 and using
(10) we obtain w = 0. Thus, b is injective on W × V̂ .
In the second step we show inf-sup stability. For given z ∈ Z ⊂ W we select uz ∈ V with Luz = M z, which
yields

b(uz , DΩh uz ; z) = uz , L∗ z 0,Ω + Luz , z 0,Ω − uz , L∗ z 0,Ω = Luz , z 0,Ω = kzk2W .

h h h h

Inserting (10) we obtain

 
k(uz , DΩh uz )kW ×Z 0 = max kuz kW , kDΩh uz kZ 0 ≤ max kuz kW , kuz kH(L,Ωh )
q q
≤ kuz kH(L,Ωh ) ≤ CL2 + 1kM −1 Luz kW = CL2 + 1kzkW .

This establishes for all z ∈ Z by selecting (uz , DΩh uz ) and (M −1 L∗ z, 0) in W × Z 0

b(M −1 L∗ z, 0; z)
 
b(v, v̂; z) b(uz , DΩh uz ; z)
sup ≥ max , (14)
(v,v̂)∈W ×V̂ k(v, v̂)kW ×Z
0 k(uz , DΩh uz )kW ×Z 0 k(M −1 L∗ z, 0)kW ×Z 0
( )
kzk2W kM −1 L∗ zk2W
≥ max p 2 ,
CL + 1kzkW kM −1 L∗ zkW
max kzkW , kM −1 L∗ zkW

kzkZ
≥ p
2
≥√ p 2 .
CL + 1 2 CL + 1

Since b(·; ·) is injective in W × V̂ , we obtain (13) by duality [1, Lem. 4.4.2]. 

8

4. Petrov-Galerkin estimates

In Lem. 4 we established a variational space-time setting for weak solutions in the product space Y = W × V̂
with a broken test space Z. Introducing the trial-to-test operator TZ , see Tab. 1, we observe

b(y; z)
kTZ ykZ = sup ≥ βkykY , y∈Y ,
z∈Z kzkZ

i.e., inf-sup stability of the bilinear form b(·; ·) in Y × Z implies ellipticity kTZ yk2Z = hSy, yi ≥ β 2 kyk2Y of the
corresponding symmetric Schur complement problem in Y .

Table 1. Operators for the Petrov-Galerkin analysis.

B ∈ L(Y, Z 0 ) hBy, zi = b(y; z) y∈Y , z∈Z

Riesz operator AZ ∈ L(Z, Z 0 ) hAZ z, ψi = (z, ψ)Z z, ψ ∈ Z
trial-to-test operator TZ = A−1Z B ∈ L(Y, Z) (T y, z)Z = b(y; z) y∈Y , z∈Z
Schur complement S = B 0 A−1 0
Z B ∈ L(Y, Y )
natural embedding EYh ∈ L(Yh , Y ) EYh yh = yh yh ∈ Yh
EZh ∈ L(Zh , Z) EZh zh = zh zh ∈ Zh
AZh = EZ0 h AZ EZh ∈ L(Zh , Zh0 ) hAZh zh , ψh i = (zh , ψh )Z zh , ψh ∈ Zh
Galerkin projection PZh = A−1 0
Zh EZh AZ ∈ L(Z, Zh ) (PZh z, ψh )Z = (z, ψh )Z z ∈ Z , ψ h ∈ Zh
Bh = EZ0 h BEYh ∈ L(Yh , Zh0 ) hBh yh , zh i = b(yh ; zh ) yh ∈ Yh , zh ∈ Zh
TZh = A−1Zh Bh = PZh TZ EYh ∈ L(Yh , Zh ) (TZh yh , zh )Z = b(yh ; zh ) yh ∈ Yh , zh ∈ Zh
Sh = Bh0 A−1 0
Zh Bh ∈ L(Yh , Yh )

Now we select discrete spaces Wh ⊂ W and Vh ⊂ V , and we set Yh = Wh × DΩh (Vh ) ⊂ Y . An appropriate
discrete test space always exists (see, e.g., [3, Thm. 4.8] for the transport equation and [24, Sect. 6]).

Lemma 5. For given βh ∈ (0, β) a discrete test space Zh ⊂ Z exists such that

b(yh ; zh )
kTZh yh kZ = sup ≥ βh kyh kY , yh ∈ Yh . (15)
zh ∈Zh kzh kZ

Again this implies ellipticity hSh yh , yh i ≥ βh2 kyh k2Y of the corresponding discrete Schur complement problem.

Proof. Let Zhopt = TZ EYh (Yh ) be the optimal test space, and let Zh,k , k ∈ N be a dense family of discrete spaces
so that z = lim PZh,k z for all z ∈ Z. Since Zhopt is discrete, for all ε ∈ (0, 1) some k = k(ε) ∈ N exists such
k−→∞
that kPZh,k zh − zh kZ ≤ εkzh kZ for all zh ∈ Zhopt . This yields

kPZh,k zh k2Z = kzh k2Z − kPZh,k zh − zh k2Z ≥ kzh k2Z − ε2 kzh k2Z = (1 − ε2 ) kzh k2Z , zh ∈ Zhopt .
p
Now, for given βh ∈ (0, β), select ε = 1 − βh2 /β 2 > 0 and set Zh = Zh,k(ε) . This yields for all yh ∈ Yh

b(yh ; zh ) p p
sup = kTZh yh kZ = kPZh TZ EYh yh kZ ≥ 1 − ε2 kTZ EYh yh kZ ≥ 1 − ε2 βkyh kY = βh kyh kY .
zh ∈Zh kzh kZ 
 Table 2. Operators and bilinear forms for the local Petrov-Galerkin analysis with restricted 9
spaces YK = L2 (K; Rm ) × H(K, L) and ZK = H(K, L∗ ) for K ⊂ Ωh , and local finite element
spaces WK,h = Wh |K , VK,h = Vh |K , YK,h = WK,h × VK,h , and ZK,h = Zh |K .

0
BK ∈ L(YK , ZK ) hBK yK , zK i = bK (yK ; zK ) yK ∈ YK , zK ∈ ZK
0
BK,h ∈ L(YK,h , ZK,h ) hBK,h yK,h , zK,h i = bK (yK,h ; zK,h ) yK,h ∈ YK,h , zK,h ∈ ZK,h
0
AZK,h ∈ L(ZK,h , ZK,h ) hAZK,h zK,h , ψK,h i = (zK,h , ψK,h )ZK zK,h , ψK,h ∈ ZK,h
0
AYK,h ∈ L(YK,h , YK,h ) hAYK,h yK,h , φK,h i = (yK,h , φK,h )YK yK,h , φK,h ∈ YK,h
EYK,h ∈ L(YK,h , YK ) EYK,h yK,h = yK,h yK,h ∈ YK,h
0
SK,h = BK,h A−1 0
ZK,h BK,h ∈ L(YK,h , YK,h )

Since the optimal test space Zhopt is not accessible, the proof in Lem. 5 is not constructive, and the norm in
Yh = Wh × V̂h ⊂ W × Z 0 cannot be evaluated. Nevertheless, for broken test spaces Zh =
Q
ZK,h the well-
posedness of the discrete problem can be tested by a local criterion, and norm estimates can be evaluated in
WK,h × VK,h = (Wh × Vh )|K using the norm in VK,h ⊂ H(L, K), so that for the computable estimates it is not
0
required to evaluate the norm in the trace space V̂K,h = DK (VK,h ) ⊂ ZK,h .
Therefore we introduce local operators (see Tab. 2), the local bilinear form bK (·; ·) defined by

bK (wK , vK ; zK ) = (wK , L∗ zK )0,K + hDK vK , zK i , (wK , vK ) ∈ YK , zK ∈ ZK ,

and for all zK ∈ ZK the local affine spaces


ZK,h (zK ) = zK,h ∈ ZK,h : bK (yK,h , zK − zK,h ) = 0 for all yK,h ∈ YK,h .

Lemma 6. If ZK,h (zK ) is not empty for all zK ∈ ZK and K ⊂ Ωh , the operator Bh is injective in Yh and a
Fortin operator exists, i.e., a projection Πh ∈ L(Z, Zh ) with b(yh ; z − Πh z) = 0 for all yh ∈ Yh and z ∈ Z.

This provides discrete stability and the estimate βh ≥ β/kΠh kZ for the inf-sup constant [2, Prop. II.2.8].

Proof. For zK ∈ ZK and K ⊂ Ωh we define ΠK,h z ∈ ZK,h (zK ) as the element with minimal norm. Therefore,
we compute a critical point (zK,h , yK,h ) ∈ ZK,h × YK,h of the corresponding local Lagrange functional

1
FK,h (zK,h , yK,h ) = kzK,h k2ZK + bK (yK,h ; zK − zK,h )
2
1
= hAZK,h zK,h , zK,h i + hBK EYK,h yK,h , zK i − hBK,h yK,h , zK,h i ,
2
0
i.e., AZK,h zK,h = BK,h yK,h and BK,h zK,h = EY0 K,h BK
0
zK . This yields

0
SK,h yK,h = BK,h A−1 0 0
ZK,h BK,h yK,h = EYK,h BK zK . (16)

We have

hEY0 K,h BK
0
zK , ψK,h i = bK (ψK,h ; zK ) = hBK,h ψK,h , EZK,h zK i = 0 , ψK,h ∈ N (BK,h ) ,

i.e., EY0 K,h BK

0
zK ∈ N (BK,h )⊥ . Since SK,h is self-adjoint and N (SK,h ) = N (BK,h ), this shows that a Lagrange
multiplier yK,h ∈ YK,h solving the local Schur complement problem (16) exists, but in general, the solution is
not unique. Thus we select the solution with minimal norm in YK,h . This is determined by the pseuso-inverse
with respect to the topology in YK
 −1
+
SK,h = lim SK,h A−1
YK,h SK,h + δAYK,h SK,h A−1 0
YK,h ∈ L(YK,h , YK,h ) .
δ−→0
10

+
Then, yK,h = SK,h EY0 K,h BK
0
zK is a Lagrange multiplier and the minimizer is given by zK,h = ΠK,h zK with

ΠK,h = A−1 + 0 0
ZK,h BK,h SK,h EYK,h BK ∈ L(Z, ZK,h ) .
P  P
This defines Πh z = K,h EZK,h ΠK,h z|K , and thus b(yh , z − Πh z) = K bK (yh |K , z|K − ΠK,h z|K ) = 0 for
K
all yh ∈ Yh , i.e., Πh is a Fortin operator. Then, for all yh ∈ N (Bh )

0 = b(yh ; Πh z) = b(yh ; z) = hBEYh yh , zi , z ∈Z.

Since B is injective in Y , this implies yh = 0 and thus the assertion. 

Remark 7. If Bh is injective, a Fortin operator exists (see, e.g., [10] for a general Banach space case and [6]
for the application to DPG). An optimal Fortin operator Πopt h can be determined as follows: for given z ∈ Z
find zh ∈ Zh with minimal norm kzh kZ subject to the constraint b(yh ; z − zh ) = 0 for all yh ∈ Yh . Again, this
can be computed from a critical point (zh , yh ) ∈ Zh × Yh of the corresponding Lagrange functional
1
F (zh , yh ) = kzh k2Z + b(yh ; z − zh ) ,
2
i.e., AZh zh = Bh yh , Bh0 zh = EY0 h B 0 z, and thus zh = A−1 0 −1 0 0 0
Zh Bh yh and Sh yh = Bh AZh Bh yh = Bh zh = EYh B z.
This yields zh = Πopt
h z with

Πopt −1 −1 0 0
h = AZh Bh Sh EYh B .

Bounding Πopt
h involves an estimate for Sh−1 which requires a bound for the inf-sup constant. The construction
+
in Lem. 6 requires only local estimates for SK,h which can be computed by local discrete symmetric eigenvalue
problems
+
AYK,h SK,h AYK,h yK,h = λK,h AYK,h yK,h , (λK,h , yK,h ) ∈ [0, ∞) × YK,h . (17)


Let k(wK , vK )kYK = max kwK kWK , kvK kVK be the norm in YK = WK × VK .
Lemma 8. If the constants αK,h > 0 satisfy
+
kSK,h `K,h kYK ≤ αK,h kA−1
YK,h `K,h kYK ,
0
`K,h ∈ YK,h , K ⊂ Ωh , (18)

the Fortin operator constructed in Lem. 6 is bounded by kΠh kZ ≤ max αK,h .

√
Solving the eigenvalue problem (17) yields the estimate αK,h = 2 max λK,h .
+
Proof. For zK ∈ ZK we define `K,h = EY0 K,h BK
0
zK and yK,h = SK,h `K,h , so that

ΠK,h zK = A−1 + 0 0
ZK,h BK,h SK,h EYK,h BK zK

= A−1 +
ZK,h BK,h SK,h `K,h

= A−1
ZK,h BK,h yK,h .

The definition of the norm in YK yields |bK (yK,h ; zK,h )| = |hBK,h yK,h , zK,h i| ≤ kyK,h kYK kzK,h kZK , zK,h ∈ ZK.h ,
which gives

kΠK,h zK kZK = kA−1

ZK,h BK,h yK,h kZK ≤ kyK,h kYK .
 11

+
From (18) we we obtain kyK,h kYK = kSK,h `K,h kYK ≤ αK,h kA−1
YK,h `K,h kYK , and finally

kA−1
YK,h `K,h kYK = kA−1 0 0
YK,h EYK,h BK zK kYK ≤ kzK kZK .

Together this yields kΠK,h zkZK,h ≤ αK,h kzkZK . Then, kΠh zk2Z ≤ 2
kz|K k2ZK ≤ max αK,h
2
kzk2Z yields the
P
αK,h
assertion. 
Remark 9. Scaling arguments on uniformly shape regular meshes in [18] and [13, Sect. 5.2] show that a bound
for αK,h can be estimated on the reference cell. Results for the acoustic wave equation on tensor-product
space-time cells are presented in [13, Sect. 7.1]

5. The realization of the DPG method

In heterogeneous materials the finite element error also depends on the approximation error of the PDE, and the
realization of the DPG method uses an approximation Mh ∈ L∞ (Ω0 ; Rm×m sym ) of M ; then we set Lh = Mh ∂t + A
∗ 1 m
and Lh z = −Lh z for z ∈ Cc (Ωh ; R ).
The estimates for the DPG analysis use functionals in DΩh (V ) ⊂ Z 0 . In the implementation, we use represen-
f = L2 (∂Ωh ; Rm ). For this purpose we introduce trace mappings for sufficiently
tations of these functionals in W
smooth functions.
Since N (DK ) = H0 (L, K) ⊃ C1c (K; Rm ), we can define trace mappings trK , tr∗K ∈ L C1 (K; Rm ), L2 (∂K; Rm )

such that

(Lh vK , zK )0,K − (L∗h vK , zK )0,K = (trK vK , tr∗K zK )0,∂K , vK , zK ∈ C1 (K; Rm ) .

f and tr∗ ∈ L C1 (Ωh ; Rm ), W

This extends to trh ∈ L C1 (Ωh ; Rm ) ∩ V, W h
f such that

(Lh v, z)0,Ωh − (L∗h v, z)0,Ωh = (trh v, tr∗h z)0,∂Ωh , v ∈ C1 (Ωh ; Rm ) ∩ V , z ∈ C1 (Ωh ; Rm )

by the selection of an orientation on inner faces F = ∂K ∩ ∂KF and the restriction trh v|F = trK v|F for
conforming functions, and the jump term tr∗h z|F = trK z|F − trKF z|F for functions in the broken space (see
[13, Sect. 6.1] for the acoustic case).
The trace defines the approximation of the bilinear form b(·; ·)

bh (wh , ṽh ; zh ) = (wh , L∗h zh )0,Ωh + (ṽh , tr∗h zh )0,∂Ωh , wh ∈ W , ṽh ∈ W

f , zh ∈ Zh .

f , and Zh ⊂ C1 (Ωh ; Rm ) such that

Now we select Wh ⊂ W , Veh ⊂ W

(ṽh , tr∗h zh )0,∂Ωh

kṽh kZh0 = sup , vh ∈ Veh
zh ∈Zh kzh kZ

is a norm in Veh and such that the bilinear form bh (·; ·) is inf-sup stable, i.e., β0 > 0 exists such that

bh (wh , ṽh ; zh )
sup ≥ β0 k(wh , ṽh )kWh ×Zh0 , (wh , ṽh ) ∈ Wh × Veh (19)
zh ∈Zh kzh kZh

with respect to the norms for w ∈ W , z ∈ Z, and (wh , ṽh ) ∈ Wh × Veh

q q
kwkWh = (Mh w, w)0,Ω0 , kzkZh = kzk2Wh + kMh−1 L∗ zk2Wh , k(wh , ṽh )kW ×Zh0 = max{kwh kW , kṽh kZh0 } .
12

For the error analysis we construct a conforming extension Vh ⊂ V of the trace space Veh so that for all ṽh ∈ Veh
a reconstruction vh ∈ Vh exist satisfying

(trh vh , tr∗h zh )0,∂Ωh = (ṽh , tr∗h zh )0,∂Ωh , zh ∈ Zh . (20)

b h vh ∈ Zh0 the corresponding functional defined by

Then, we denote by tr

b h vh , zh i = (trh vh , tr∗h zh )0,∂Ωh ,

htr zh ∈ Zh , (21)

b h (Vh ) ⊂ Zh0 .
and we set V̂h = tr
The DPG approximation (uh , ũh ) ∈ Wh × Veh for right-hand side f ∈ W minimizes the residual

bh (wh , ṽh ; zh ) − (f, zh )0,Ω

(uh , ũh ) = arg min sup .
eh zh ∈Zh
(wh ,ṽh )∈Wh ×V kzh kZh

It is computed by

bh (uh , ũh ; zh ) = (f, zh )0,Ω , zh ∈ T̃h (Wh × Veh ) , (22)

where the discrete trial-to-test operator T̃h ∈ L(Wh × Veh , Zh ) is defined by

T̃h (wh , ṽh ), zh Z
= bh (wh , ṽh ; zh ) , (wh , ṽh ) ∈ Wh × Veh , zh ∈ Zh .

Now, the first Strang Lemma [9, Lem. 2.27] takes the following form.

Theorem 10. For f ∈ W let u ∈ V be the solution of Lu = f , and let (uh , ũh ) ∈ Wh × Veh be the DPG
approximation solving (22). We assume that a conforming reconstruction space Vh ⊂ V satisfying (20) and
ûh ∈ V̂h exist with

hûh , zh i = (ũh , tr∗h zh )0,∂Ωh , zh ∈ Zh . (23)

b h u ∈ Zh0 be defined by
Let tr

b h u, zh i = (Lh u, zh )0,Ωh − (u, L∗h zh )0,Ωh ,

htr zh ∈ Zh , (24)

and we set û = DΩh u ∈ Z 0 . Then, the error is bounded by

k(u − uh , tr
b h u − ûh )kWh ×Z 0
h
≤ inf k(u − wh , tr
b h u − tr
b h vh )kWh ×Z 0
h
(wh ,vh )∈Wh ×Vh
!
bh (wh , trh vh ; zh ) − b(u, û; zh )
+ β0−1 sup .
zh ∈Zh kzh kZh

Furthermore, if u ∈ V, the error is bounded by

1 + β0−1


k(u − uh , tr
b h u − ûh )kWh ×Z 0
h
≤ inf k(u − wh , tr
b h u − tr
b h vh )kWh ×Z 0
h
(wh ,vh )∈Wh ×Vh

+ β0−1 k id −M −1 Mh k∞ k∂t ukW .

 13

Proof. For all (wh , vh ) ∈ Wh × Vh we have

k(u − uh , tr
b h u − ûh )kWh ×Z 0
h
≤ k(u − wh , tr
b h u − tr
b h vh )kWh ×Z 0 + k(wh − uh , tr
h
b h vh − ûh )kWh ×Z 0 ,
h

and inserting (21), (23), and (19) gives

k(wh − uh , tr
b h vh − ûh )kWh ×Z 0
h
= k(wh − uh , trh vh − ũh )kWh ×Zh0
bh (wh − uh , trh vh − ũh ; zh )
≤ β0−1 sup .
zh ∈Zh kzh kZh

Using bh (uh , ũh ; zh ) = (f, zh )0,Ω = (Lu, zh )0,Ω = b(u, û; zh ) for zh ∈ T̃h (Wh × Veh ) and

bh (wh , ṽh ; zh ) bh (wh , ṽh ; zh )

sup = sup , (wh , ṽh ) ∈ Wh × Veh
zh ∈Zh kzh kZh zh ∈T̃h (Wh ×V
eh ) kzh kZh

yields first estimate. The second estimate follows from bh (wh , trh vh ; zh ) = (wh , L∗h zh )0,Ωh + htr
b h vh , zh i and

bh (wh − uh , trh vh − ũh ; zh ) = bh (wh , trh vh ; zh ) − (Lu, zh )0,Ω

= (wh , L∗h zh )0,Ωh + htr


b h vh , zh i + (Mh − M )∂t u, zh − (Lh u, zh )0,Ω
0,Ω

= (wh − u, L∗h zh )0,Ωh

+ htr
b h vh − tr
b h u, zh i + (Mh − M )∂t u, zh
0,Ω
.

Remark 11. In the case that the solution u ∈ V has the additional regularity u ∈ H1+s (Ω; Rm ) with s > k ≥ 0,
a trace function ũ ∈ L2 (∂Ωh ; Rm ) with

(ũ, tr∗h z)0,∂Ωh = (Lu, z)0,Ωh − (u, L∗ z)0,Ωh , z ∈ C1 (Ωh ; Rm )

exists. In this case, the proof of Thm. 10 only relies on the traces in L2 , and no conforming reconstruction space
Vh is required for the estimates. For Wh ⊃ Pk (Ωh ) we obtain

inf ku − wh kWh ≤ kMh k1/2

∞ inf ku − wh k0,Ω ≤ C1 kMh k1/2 s
∞ h kuks,Ω ,
wh ∈Wh wh ∈Wh

and for Veh ⊃ Pk (Γh ) with Γh = F ⊂ ∂Ωh , we obtain, provided that M, Mh ∈ L∞ (∂Ωh ; Rm×m
S
sym ), the estimate

inf ktr b h vh kZ 0 = inf kũ − ṽh kZ 0 ≤ C2 h−1/2 inf kũ − ṽh k0,∂Ωh ≤ C3 hs−1 kuks,Ω
b h u − tr
vh ∈Vh h h
ṽh ∈V
eh ṽh ∈V
eh

with constants C1 , C2 , C3 > 0 depending on the shape regularity of Ωh and on Mh . Together, this results in the
convergence estimate
b h u − ûh )kWh ×Z 0 ≤ (1 + β0−1 ) max{C1 , C3 }hs−1 kuks,Ω + β0−1 k id −M −1 Mh k∞ k∂t ukW .
k(u − uh , tr h

Remark 12. Choosing an extended space Zhext ⊃ Zh allows for a conforming reconstruction Vh ⊂ Zhext ∩ V
satisfying (20), see [13, Sect. 6] for the local construction and for examples in case of acoustic waves.
In the nonconforming case Ve 6⊂ trh (Vh ) the construction of the reconstruction space Vh depends on Zh , and
the asymptotic arguments in Lem. 5 do not apply. In particular, the reconstruction space may have strong
oscillations at the corners and edges of ∂K. The numerical experiments in [12] indicate that the optimal choice
of Zh has to be well balanced to ensure discrete inf-sup stability on the one hand and to limit the nonconformity
on the other hand. Nevertheless, for a given choice of Wh , Veh , and Zh the well-posedness of the discrete system
can be guaranteed by Lem. 6, and explicitly computing a reconstruction Vh a lower bound for the inf-sup
constant can be provided by Lem. 8.
14

6. A numerical example for an application in geophysics

Many applications rely on accurate numerical simulations of waves through complex material structures. For
instance, geophysical structures like the earth’s crust below the sea bed feature complex varying material
properties. A typical example is the problem of full waveform inversion (FWI), where the material distribution
is reconstructed from measurements of the wave field close to the surface [21]. Here, in a field survey a
wave is excited at some point S0 ∈ Ω0 , and the scattered wave field is measured by receiver devices located at
R0 , . . . , RN ∈ Ω0 . During the experiment, each receiver records a time series of approximate point measurements
u(t, Rn ), t ∈ [0, T ]. The collection of all these measurements is called a seismogram, see Fig. 5 for examples.
The recorded seismogram contains information about the material structure the wave has traveled through.
Full waveform inversion techniques attempt to reconstruct from this information by applying iterative schemes
of Newton-type (see, e.g., [17]). During the iteration, a large number of wave equations has to be solved nu-
merically for different right-hand sides and varying material parameters.
x

0 9
x
R0 RN
S0
Ω0
Ω0 (T, RN )
S C+ ({S}) ∩ C− (R)
(T, R0 )

3 t
y 0 T
Figure 2. The left figure shows the spatial domain Ω0 = (0, 9) × (0, 3), the source position
S0 = (2.602, 0.802), and the first and last receiver R0 = (2.002, 0.303) and RN = (4.002, 0.303).
In the space-time cylinder, the source is located
at S = (0.4, S0 ) and the signal is measured in
R = conv {(t, R0 ), . . . , (t, RN ) : t ∈ (0, T )} , which results in the domain of interest intersecting
of the forward and backward cone C+ ({S})∩C− (R). This is illustrated on the right for {y = 0}.

To demonstrate the flexibility and the accuracy of the space-time DPG method for heterogeneous media, we
consider a numerical example corresponding to the forward problem within FWI. We use the acoustic wave
equation for the Marmousi benchmark, a synthetic model problem for geophysical structures in two space
dimensions featuring a material distribution that is similar to what is located inside the earth crust (see,
e.g., [4]). Fig. 3 shows rescaled variants of the spatial varying mass density and bulk modulus that have similar
structure to the materials in this benchmark. We select the spatial domain Ω0 ⊂ R2 , the receiver positions
R0 , . . . , RN , and the source position S = (t0 , S0 ) ∈ (0, T ) × Ω0 as shown in Fig. 2.
For the inversion, we only require the wave field in the space-time region that contributes to the transport of
information from the source to the receiver array. Thus, using Lem. 2 for the maximal wave speed cmax = 1.37,
we find a superset of the relevant region of the full space-time domain (0, T ) × Ω0 , see Fig. 2. At the source
position S, the wave is excited by the right-hand side f (x, t) = 100 φt (t0 − t)φx (|S0 − x|), (t, x) ∈ (0, T ) × Ω0
with a Ricker wavelet for f0 = 100 and r0 = 0.05 in time

φt (τ ) = (1 − 2aτ 2 ) exp(−aτ 2 ), a = π 2 f02 τ 2 ,

 6
π r
and in space φx (r) = χ[0,r0 ) cos 2 r0 .
 15

The numerical simulations were performed using the DPG approximations with local finite element spaces

WK,h = Q3 (K)3 , Ṽh |(tn−1 ,tn )×∂K = Q4 (F )3 , Ṽh |{tn }×K0 = Q4 (F )2 , ZK,h = Q6 (K)3


on a quadrilateral mesh, yielding a scheme that converges with order 4 in L2 (0, T ) × Ω0 for smooth solutions,
see [12, Sec. 4.7]. The material parameters κ and ρ are cell-wise constant as shown in Fig. 3.
The method is realized in the parallel Finite Element system M++ [22], the linear systems are solved with a
preconditioned cg method using the reduction to the symmetric positive definite Schur complement system for
the skeleton approximation Veh , cf. [23]. The implementation of the space-time DPG method is been evaluated
systematically in [12] for configurations where the analytic solution is known. Moreover, basic applications of
space-time DPG to FWI are discussed in [12].

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 1.9 2.0 2.1 2.2 2.3 2.4 2.5

Figure 3. Approximations of the spatial material distributions for Marmousi, κ on the left
and ρ on the right. The first row corresponds to the discrete material parameters on a coarse
quadrilateral mesh with Nx · Ny = 144 · 48 = 6 912 cells, the second row shows the material on
a refined mesh with Nx · Ny = 288 · 96 = 27 648 cells.

The pressure component of the full numerical space-time solution is shown in Fig. 4, also see Fig. 6. This
is compared with the solution on a mesh that has been truncated to the space-time region of interest. The
seismograms corresponding to the materials are shown in Fig. 5.
The numerical experiments
with the truncated space-time cylinder correspond to our results in Lem. 2. Up to
a small relative L2 0, T ; R3 difference of 0.00318, the seismograms on the full mesh and the truncated mesh
coincide. The linear system of all space-time unknowns is solved approximately using an iterative solver, which
can explain the small difference. On the other hand, only explicit time stepping methods guarantee finite speed
of the discrete wave propagation. Since DPG is an implicit scheme, this also may explain a small difference of
the results for the full space-time cylinder and the truncated domain.
The comparison of the solutions on different refinement levels of the mesh shows that both seismograms have a
similar structure at the beginning but differ strongly at later times. Here, we have two possible reasons for this
difference: on the one hand, the approximation quality for the wave field increases on mesh refinement, and on
the other hand, on the finer mesh we have a better resolution of the material parameters Mh , i.e., a different
wave equation is solved, see Thm. 10.
16

Figure 4. The left picture shows the pressure component of the numerical solution on a regular
grid with Nt · Nx · Ny = 128 · 144 · 48 = 884 736 space-time cells resulting in 156 576 000 space-
time DoFs for all components in the skeleton space Ṽh . On the right, a truncated version of the
space-time cylinder is depicted that respects the position of the source and the receivers. The
truncation reduced the amount of cells by approx. 43% to 503 024 cells and 89 394 900 DoFs.

Conclusion and Outlook. In this example we show that the method also yields accurate results for application
driven simulations. We demonstrate that for wave-type equations the finite wave speed can be exploited to
reduce the size of the linear system considerably while yielding the same results. Further numerical experiments
are necessary to compare the roles of the discretization error for fixed material M and h → 0 to the model error
resulting from the approximate material Mh . For that reason, we intend to evaluate the performance of an
established discretization like discontinuous Galerkin with time stepping in comparison to the space-time DPG
method.

Acknowledgment. We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft

(DFG) through CRC 1173.

References
[1] D. Braess. Finite Elements. Theory, fast solvers, and applications in solid mechaics. 3th ed. Cambridge
University Press, 2007.
[2] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer, 1991.
[3] D. Broersen, W. Dahmen, and R. Stevenson. On the stability of DPG formulations of transport equations.
Mathematics of Computation, 87(311):1051–1082, 2018.
[4] A. Brougois, M. Bourget, P. Lailly, M. Poulet, P. Ricarte, and R. Versteeg. Marmousi, model and data. In
EAEG Workshop-Practical Aspects of Seismic Data Inversion, 1990.
[5] C. Carstensen, L. Demkowicz, and J. Gopalakrishnan. Breaking spaces and forms for the DPG method
and applications including Maxwell equations. Computers & Mathematics with Applications, 2016.
 17

R28 R28

R24 R24

R20 R20

R16 R16

R12 R12

R8 R8

R4 R4

R0 R0
0 1 2 3 4 5 6 7 8 t 0 1 2 3 4 5 6 7 8 t

full regular
truncated refined

R0 R0

0 1 2 3 4 5 6 7 8 t 0 1 2 3 4 5 6 7 8 t
R
Figure 5. All plots show measurements of the form t 7−→ φ(s − t, x − Rn ) · u(s, x) d(s, x) of
the pressure component at the corresponding receivers. Here φ is a measurement kernel with
small support in space and time. On the left-hand side, a comparison of the seismograms on
the full mesh (solid) and the truncated mesh (dotted) is depicted. At the bottom, the recording

at R0 is shown for illustration purposes. Since the difference is very small (relative L2 0, T ; R3
difference is 0.00318), both seismograms look very similar. On the right, we compare the
seismogram corresponding to the truncated mesh (solid) to the seismogram obtained on a
mesh that has been uniformly refined (dotted). While being similar to the measurements on
the coarser grid, we see a different fine structure due to the different material approximation Mh .
18

Figure 6. Pressure distribution for three different numerical simulations at selected time steps
t = 0.4, 1.49, 2.57, 3.66, 4.74, 5.83, 6.91, 8. The first column shows the results for the space-
time DPG method using the full space-time cylinder on level 4. The second and third column
correspond to space-time DPG in the truncated space-time cylinder on level 4 and level 5
respectively.

[6] C. Carstensen and F. Hellwig. Low-order discontinuous Petrov–Galerkin finite element methods for linear
elasticity. SIAM Journal on Numerical Analysis, 54(6):3388–3410, 2016.
[7] L. F. Demkowicz and J. Gopalakrishnan. An overview of the discontinuous Petrov–Galerkin method. In
Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations,
 19

pages 149–180. Springer, 2014.

[8] L. F. Demkowicz, J. Gopalakrishnan, S. Nagaraj, and P. Sepulveda. A spacetime DPG method for the
Schrödinger equation. SIAM J. Numer. Anal., 55(4):1740–1759, 2017.
[9] A. Ern and J.-L. Guermond. Theory and practice of finite elements, volume 159. Springer Science &
Business Media, 2013.
[10] A. Ern and J.-L. Guermond. A converse to Fortin’s Lemma in Banach spaces. Comptes Rendus Mathema-
tique, 354(11):1092–1095, 2016.
[11] A. Ern, J.-L. Guermond, and G. Caplain. An intrinsic criterion for the bijectivity of Hilbert operators
related to Friedrichs systems. Communications in partial differential equations, 32(2):317–341, 2007.
[12] J. Ernesti. Space-Time Methods for Acoustic Waves with Applications to Full Waveform Inversion. PhD
thesis, Karlsruher Institut für Technologie (KIT), 2018.
[13] J. Ernesti and C. Wieners. A space-time DPG method for acoustic waves. In U. Langer and O. Stein-
bach, editors, Space-Time Methods. Applications to Partial Differential Equations, volume 21 of Radon
Series on Computational and Applied Mathematics. Walter de Gruyter, 2018. submitted, preprint
http://www.math.kit.edu/user/~wieners/SpaceTimeDPG_revised.pdf.
[14] L. C. Evans. Partial differential equations, 2. ed. American Mathematical Society, Providence, RI, 2010.
[15] J. Gopalakrishnan and P. Sepulveda. A spacetime DPG method for acoustic waves. arXiv:1709.08268,
2017.
[16] M. Hochbruck, T. Pazur, A. Schulz, E. Thawinan, and C. Wieners. Efficient time integration for discon-
tinuous Galerkin approximations of linear wave equations. ZAMM, 95(3):237–259, 2015.
[17] L. Métivier, R. Brossier, S. Operto, and J. Virieux. Acoustic multi-parameter fwi for the reconstruction of
p-wave velocity, density and attenuation: preconditioned truncated newton approach. In SEG Technical
Program Expanded Abstracts 2015, pages 1198–1203. Society of Exploration Geophysicists, 2015.
[18] S. Nagaraj, S. Petrides, and L. F. Demkowicz. Construction of DPG Fortin operators for second order
problems. Computers & Mathematics with Applications, 2017.
[19] M. Renardy and R. C. Rogers. An introduction to partial differential equations, volume 13. Springer Science
& Business Media, 2006.
[20] W. Rudin. Functional Analysis. International Series in Pure and Applied Mathematics. McGraw-Hill, Inc.,
New York, 1991.
[21] W. Symes. The seismic reflection inverse problem. Inverse problems, 25(12):123008, 2009.
[22] C. Wieners. A geometric data structure for parallel finite elements and the application to multigrid methods
with block smoothing. Computing and visualization in science, 13(4):161–175, 2010.
[23] C. Wieners. The skeleton reduction for finite element substructuring methods. In Numerical Mathematics
and Advanced Applications ENUMATH 2015, pages 133–141. Springer, 2016.
[24] C. Wieners and B. Wohlmuth. Robust operator estimates and the application to substructuring methods
for first-order systems. ESAIM: M 2 AN, 48:161–175, 2014.