Space-Time Discontinuous Galerkin Methods for Weak Solutions of Hyperbolic Linear Symmetric Friedrichs Systems

We study weak solutions and its approximation of hyperbolic linear symmetric Friedrichs systems describing acoustic, elastic, or electro-magnetic waves. For the corresponding first-order systems we construct discontinuous Galerkin discretizations in space and time with full upwind, and we show primal and dual consistency. Stability and convergence estimates are provided with respect to a mesh-dependent DG norm which includes the L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{L}}_2$$\end{document} norm at final time. Numerical experiments confirm that the a priori results are of optimal order also for solutions with low regularity, and we show that the error in the DG norm can be closely approximated with a residual-type error indicator.


Introduction
Linear wave equations are hyperbolic, and the formulation as first-order symmetric Friedrichs system provides a well established setting for analyzing and approximating solutions.A specific feature of hyperbolic systems is the transport of discontinuities along characteristics.Our goal is to provide a numerical scheme which is efficient for smooth solutions as well as for weak solutions with discontinuities.
For smooth solutions of linear symmetric Friedrichs systems O(h s−1/2 ) convergence can be established for discontinuous Galerkin approximations in space with respect to suitable mesh-dependent DG norm [9,Chap. 57], [5,Chap. 7].For acoustics, the convergence analysis of a space-time approximation in a DG semi-norm provides estimates for all discrete time steps [2,Prop. 6.5].
Finite volume convergence O(h 1/2 ) for hyperbolic linear symmetric Friedrichs systems is established in [18] combined with first-order time-stepping.Discontinuous Galerkin methods in time are analyzed in [12] for tent-type space-time meshes.This is adapted to space-time discontinuous Galerkin methods on general space-time meshes with upwind flux for acoustics in [2], where the convergence is established for sufficiently smooth solutions based on estimates in a suitable DG semi-norm.In particular, the analysis includes the adaptive approximation of corner singularities.
Here, we consider a DG method in space and time for linear symmetric Friedrichs systems, and we show inf-sup stability and convergence in the DG norm.Therefore we transfer our results for space-time Petrov-Galerkin methods in [6,7] with continuous approximations in time and for the DPG method in [10,11], where convergence in a stronger graph norm is considered.Our analysis includes bounds for the consistency error in the case that piecewise discontinuous material parameters are not aligned with the mesh.Convergence in the limit for piecewise discontinuous solutions of Riemann problems is established only in L 2 .
The space-time method is realized in the parallel finite element system M++ [4].In our numerical examples we confirm the a priori estimates for weak as well as for smooth solutions, and we demonstrate the efficiency of the p-adaptive scheme.
Space-time computations have a long history in practical engineering applications and in parallel time integration [26,13].The space-time approach allows for large-scale parallel computing and in case of point sources the reduction to the time cone within the space-time cylinder.Moreover, it allows for dual-primal goal-oriented error control and applications to inverse and optimal control problems where the adjoint problem is backward in time and relies on the forward solution in the full space-time cylinder.Space-time discretizations for the wave equation are constructed within a second-order approach in [25,19], with isogeometric methods in [27], a very weak approach is presented in [15], a quasi-Trefftz method is considered in [17], and a new approach to spacetime boundary integral equations for the wave equation is developed in [24].
In comparison with these methods the first-order DG approach is numerically expensive.On the other hand, convergence can be established with minimal regularity assumptions, the method easily extends to more general material laws and to more general hyperbolic conservation laws.
The paper is organized as follows.In Sect. 2 we introduce the notation and the formulation of wave equations as first-order systems, in Sect. 3 we introduce the DG discretization in time and in space.In Sect. 4 we consider well-posedness and stability, in Sect. 5 we prove existence of weak solutions and convergence estimates, in Sect.5.3 we introduce an a posteriori error indicator, and in Sect.6 we present numerical results.In Sect.7 we conclude with a discussion of possible extensions and open problems.

Symmetric Friedrichs systems
We consider weak solutions of linear hyperbolic first-order systems in the form of symmetric Friedrichs systems.Let Ω ⊂ R d be a bounded domain in space with Lipschitz boundary ∂Ω, I = (0, T ) a time interval, and we denote the space-time cylinder by Q = (0, T )×Ω.Boundary conditions will be imposed on Γ k ⊂ ∂Ω for k = 1, . . ., m depending on the model, where m is the dimension of the first-order system.
For S ⊂ Q the L 2 norm and inner product are denoted by • S and (•, •) S .
Let L = M ∂ t + A be a linear differential operator in space and time, where (M v)(t, x) = M (x)v(t, x) defines the operator M with a uniformly positive definite matrix-valued function M ∈ L ∞ (Ω; R m×m sym ), and where Av = d j=1 A j ∂ j v is a differential operator in space with matrices A j ∈ R m×m sym .Since M is uniformly positive definite, constants C M ≥ c M > 0 exists such that c M y y ≤ y M (x)y ≤ C M y y , y ∈ R m and a.a.x ∈ Ω .
We observe so that L * = −L is the adjoint differential operator.This is now complemented by initial and boundary conditions.For the unit normal vector n ∈ L ∞ (∂Ω; R d ) we define the matrix Correspondingly, we get for the operator L in space and time In order to define weak solutions, we include initial values for t = 0 and boundary conditions on Γ k for k = 1, . . ., m in the right-hand side.Therefore, we use a test space The property (1) characterizes adjoint boundaries Γ * k ⊂ ∂Ω for k = 1, . . ., m, so that the test space is defined by with homogeneous final values at t = T and homogenous values at the adjoint boundaries.
Our aim is to find a weak solution with This is now specified for acoustic, elastic and electro-magnetic waves.

Acoustic waves
The second-order wave equation is considered as first-order system with p = ∂ t φ and q = −κ∇φ, i.e., for volume data b, boundary data g N , p D , initial data q 0 , p 0 , positive parameters , κ, and the disjoint decomposition of the boundary ∂Ω = Γ D ∪ Γ N into Dirichlet and Neumann part.The corresponding Friedrichs system with m = 1 + d components is given by so that for smooth functions ϕ, ψ with ϕ = 0 on (0, T ) × Γ D and n In two space dimensions, this corresponds to the boundary parts Elastic waves Linear elastic waves are described by the first-order system for velocity v and stress σ with mass density , the symmetric gradient ε = ε(v) of v, and, in isotropic media, with Cε = 2µε + λ trace(ε)I 3 depending on the Lamé parameters µ, λ > 0. This corresponds to the Friedrichs system with For d = 3 we have m = 9 and Electro-magnetic waves The first-order system for the electric field E and the magnetic field intensity H with permittivity ε, permeability µ, and boundary decomposition ∂Ω = Γ E ∪ Γ M corresponds to a Friedrichs system with For d = 3 we have m = 6 and Remark 1 We only consider the case that the symmetric matrices A j , j = 1, . . ., d, are constant in Ω.In general, A j may depend on x ∈ Ω, e.g., for the linear transport equation Lu = ∂ t u + a • ∇u with m = 1 and transport vector a(x) ∈ R d .Then, Γ 1 is the inflow boundary, and for the adjoint equation we obtain For the DG analysis of this case we refer to [5,Chap. 2] in the steady case and to [6] for a Petrov-Galerkin space-time method.
The suitable choice of the subsets Γ k ⊂ ∂Ω for k = 1, . . ., m for the boundary conditions in general Friedrichs systems is discussed in [5,Chap. 7.2].Here we consider the special case for wave systems.The property (1) characterizes the adjoint boundaries Γ * k ⊂ ∂Ω for k = 1, . . ., m, and we observe m k=1 ) and w = (w 1 , . . ., w m ) ∈ V * and thus, defining with homogeneous initial value at t = 0 and homogeneous boundary values on Γ k , we obtain Boundary conditions are required in order to obtain uniqueness and wellposedness of the solution.Therefore, we require for the subsets Γ k ⊂ ∂Ω, for k = 1, . . ., m, that the operators L and L * are injective on V and V * , respectively, i.e., where the relatively open adjoint boundaries Γ * k ⊂ ∂Ω for k = 1, . . ., m are determined by property (1).Now we show that both conditions in (8) are necessary.The first condition for Γ k is required for uniqueness for strong solutions: if v ∈ V \ {0} exists with Lv = 0, then this is a non-trivial homogeneous strong solution, i.e., v solves (3) with u 0 = 0, f = 0, and g = 0. On the other hand, if the second condition is violated, weak solutions do not exist for all volume data: if w ∈ V * \ {0} and f ∈ L 2 (Q; R m ) exists with L * w = 0 and (f , w) Q = 0, no weak solution of (2) with homogeneous initial and boundary data u 0 = 0 and g = 0 exists.Remark 3 For smooth domains and data, the solution is also smooth, e.g., for acoustics φ(t) ∈ H s (Ω) for all t ∈ [0, T ] with s ≥ 2. This allows for improved approximation orders O(h s ) for φ.On the other hand, the necessary regularity requirements are quite restrictive [21], and the second-order formulation does not allow for the convergence analysis of piece-wise discontinuous solutions.
Remark 4 Waves in real media are dissipative and dispersive; e.g., modeling electro-magnetic waves in matter needs to include conductivity and impedance.The DG analysis can be extended to this case; see, e.g., [5,Chap. 7] for the steady case and [8] for visco-elastic waves with impedance boundary conditions.In the elastic model for Rayleigh damping or for the Kelvin-Voigt model, the linear operator takes the form L = M ∂ t + D + A with (Dv)(t, x) = D(x)v(t, x) and D ∈ L ∞ (Ω; R d×d sym ) symmetric positive semi-definite; then, All our subsequent results extend to this case, but for simplicity we only consider the case D = 0.

The full-upwind discontinuous Galerkin discretization
In this section we introduce an upwind DG discretization for the first-order system.
3.1 The DG finite element space in the space-time cylinder For the discretization, we use tensor product space-time cells combining the mesh in space with a decomposition in time.For 0 = t 0 < t 1 < • • • < t N = T , we define time intervals I n,h = (t n−1 , t n ), time-step sizes t n = t n − t n−1 , and We set t = max t n , and we assume quasi-uniformity, i.e., t n ∈ [C sr t, t] with C sr ∈ (0, 1] independent of N .Let K h be a mesh so that of the space-time cylinder Q.Let F ∈ F K be the faces of the element K, and we set F h = K F K , so that ∂Ω h = F ∈F h F is the skeleton in space; ∂Q h = N n=0 {t n } × ∂Ω h is the corresponding space-time skeleton.For inner faces F ∈ F h ∩ Ω and K ∈ K h , let K F be the neighboring cell such that F = ∂K∩∂K F .On boundary faces F ∈ F h ∩∂Ω we set K F = K.Let n K be the outer unit normal vector on ∂K.We assume that Ω = Ω h ∪ ∂Ω h and that the boundary decomposition is compatible with the mesh, i.e., We set h K = diam K, h F = diam F , and h = max h K .We assume quasiuniform meshes and shape-regularity, i.e., h F ≥ c sr h K for F ∈ F K with c sr > 0 independent of h K .In the following, we use the mesh-dependent norms In order to calibrate the accuracy in space and time, we assume, depending on a reference velocity c ref > 0, that the mesh size in time and space are well balanced satisfying Since we only consider fully implicit methods, we have no restriction with respect to stability of the time integration.
Remark 5 For simplicity we use only tensor-product space-time meshes.For the extension to more general meshes in the space-time cylinder we refer to [14], see also the analysis in [2].General meshes in R 1+d are considered in [22].
Then, the condition (10) can be relaxed to a local condition.
The DG discretization is defined for a finite dimensional subspace On the space-time skeleton ∂Q h , we define For the positive definite matrix function sym ) be a piecewise constant approximation, and for K ∈ K h let M h,K ∈ R m×m sym be the continuous extension of M h | K to K; in case of material jumps this can result to different values on the left and right side of a face, i.e., Let L h = M h ∂ t +A be the corresponding linear differential operator, where the approximated operator M h is given by ( For our applications, we use a tensor-product construction of the finite element space.
For every space-time cell R = I n,h × K we select polynomial degrees p R = p n,K ≥ 0 in time and q R = q n,K ≥ 0 in space.With this we define the discontinuous finite element spaces where P p denotes the set of polynomials up to order p.For the following, we fix p = max p R and q = max q R , so that On the space-time skeleton ∂Q h = N n=0 {t n } × Ω ∪ I h × ∂Ω h , the inverse inequality and the discrete trace inequality [5, Lem.1.44 and Lem.1.46] yield with C inv , C tr > 0 depending on the space-time mesh regularity (and thus also on c ref ), the polynomial degrees in V h , and the material parameters. Let For In every time interval I n,h we use the projection In the following, we derive the discretizations in the infinite dimensional piecewise continuous spaces S h and V h , since several properties only rely on the mesh.We use the finite dimensional DG spaces V h ⊂ V h and S h ⊂ S h if we require additional properties of the discrete space such as inverse and trace inequalities.

A discontinuous Galerkin method in time
For v h , w h ∈ V h we obtain after integration by parts in all intervals

Introducing the jump terms [w
We have dual consistency by construction, i.e.,

Again integrating by parts and defining [v
Together, we obtain so that This extends to discontinuous functions in V h as follows.

Lemma 1 We have
Proof The assertion follows from

A discontinuous Galerkin method in space
For v h , w h ∈ S h we observe, integrating by parts for all elements For conforming functions v, we have for the flux v on inner faces F ⊂ Ω, and for discontinuous functions we define the jump term [w h ] K,F = w h,K F − w h,K .On boundary faces F ⊂ ∂Ω this depends on the boundary conditions, and we set We use the discontinuous Galerkin with full upwind discretization in space which is of the form where the upwind flux A up n K ∈ R m×m is obtained by solving local Riemann problems.
For the DG method we require dual consistency for the bilinear form and the right hand side for the boundary values for and for the inconsistency complement we require that C 1 ≥ c 1 > 0 exists such that We assume that C 1 > 0 only depends on the material parameters, and that For acoustic, elastic and electro-magnetic waves the upwind flux is explicitly evaluated, e.g., in [16,Sect. 4.3].Here, we only consider the dual representation; integration by parts yields the primal representation.

Acoustic waves
The full upwind DG approximation for the acoustic wave equation ( 5) is given by the piecewise constant approximations for the material parameters κ, > 0.
On inner boundaries material discontinuities can result in Z K = Z K F , on boundary faces we define Z h = Z K on ∂Ω ∩ ∂K.On Dirichlet boundary faces The right-hand side is complemented by the stabilization, so that Integration by parts gives .

Elastic waves
The full upwind DG approximation for the elastic wave equation ( 6) is given by are the impedance of compressional waves and shear waves, respectively.On Dirichlet boundary faces Integrating by parts yields Electro-magnetic waves The full upwind DG approximation for the electromagnetic wave equation ( 7) is given by

and on impedance boundary faces
Again, integration by parts yields .

A discontinuous Galerkin method in time and space
Combining the two semi-discretizations, we obtain the full DG discretization with right-hand side in the space-time cylinder for For the space-time DG method we have by construction dual consistency for the bilinear form and the right hand side and and positivity for the inconsistency complement by combining ( 18) and ( 21).Together with ( 19) and ( 22) we obtain for v h , w h ∈ V h , and ( 23) yields with Proof By inserting v h (t) into (21) and integrating over time we find and thus with Lem. 1 we get for all v h ∈ V h

Well-posedness and stability
We show that the discrete problem has a unique solution and is stable with respect to different norms.

Well-posedness of the space-time DG discretization
The well-posedness of the discrete equation is now established as in [2, Prop.5.1].
Proof Since dim V h < ∞, it is sufficient to show that u h = 0 is the unique solution of the homogeneous problem Since (37) implies b h (u h , u h ) = 0, we obtain by (32) for the jump terms , L h u h = 0. Now the assertion follows from Lem. 2 and (33) by The previous lemma shows that the discrete graph norm defined by is well defined and a norm in V h .Since the discrete graph norm is only a semi-norm in V h , we have to use stronger norms for the convergence analysis.

Stability in space and time
Let 0 = c p,0 < c p,1 < • • • < c p,p < 1 be the Radau Ia collocation points, so that (with quadrature weights ω p,k > 0 for k = 0, . . ., p), and let λ p,k ∈ P p be the corresponding Lagrange polynomials Together this is combined to the corresponding interpolation I h : For the interpolation we will use in the following the estimate Lemma 4 If p n,K = p n for all K ∈ K h and n = 1, . . ., N , we have for and together with Lem. 1 we obtain For the upwind DG discretization in space we obtain by ( 21) so that together we obtain the assertion by Remark 7 Together with (38) and (39) we obtain L 2 stability with respect to the discrete graph norm by the discrete solution (36), and assume homogeneous boundary data g = 0.If p n,K = p n for all K ∈ K h and n = 1, . . ., N , the solution is bounded by , so that together with Lem. 4 and I h (d T u h )(0) = T u h (0) we get the assertion by Remark 8 The estimate in Lem. 2 directly implies that the Petrov-Galerkin method with test space V * h = d T V h is well-defined and L 2 stable: the Petrov-Galerkin solution and thus, in case of homogeneous boundary data g = 0 we obtain This is proposed and analyzed in [1] in the semi-discrete case for the advectiondiffusion problem.Our numerical tests indicate, that the Petrov-Galerkin modification does not improve the approximation quality, and in the next section we show, that stability and convergence in the DG norm can be established also for the Galerkin method with ansatz and test space V h and with adaptively chosen p n,K .

Inf-sup stability in the DG norm
Suitable mesh-dependent DG semi-norms and norms can be defined for all 2

and 7]. Analogously to the proof of Lem. 3 we observe that
We have i.e., v h h,DG ≤ v h h,DG + , and continuity of the bilinear form b The inf-sup stability for the advection equation [5,Lem. 2.35] can be transferred to our setting.

and we obtain by the discrete trace inequality (13b)
and together with the inverse inequality (13a) this yields We observe, using (42), Using (43), we obtain the assertion with 5 Convergence of the DG space-time approximation In the first step, we show that stability in L 2 implies convergence in the limit of the DG approximation.Then, by assuming some regularity of the solution, qualitative convergence results are obtained in the DG norm.

Convergence in the limit
Let Q h h∈H be a shape-regular family of space-time meshes with mesh sizes For h ∈ H, let u h ∈ V h be the solution of the discrete problem (36).The proof of existence of a unique discrete solution in Lem. 3 only relies on the properties (32) and (33) of the DG bilinear form and thus only implicitly on the boundary parts Γ k ⊂ ∂Ω.In order to obtain a unique weak solution of (2) in the limit, constraints for the selection of Γ k ⊂ ∂Ω, k = 1, . . ., m, are necessary, cf.(8).This is used in the following.
Theorem 2 Assume that p n,K = p n ≥ 1 and q n,K ≥ 1.In case of homogeneous boundary data g = 0 and convergent approximations of the material parameters , the discrete solutions u h h∈H are converging to a weak solution u ∈ L 2 (Q; R m ) of (2).Moreover, u is a strong solution satisfying (3), and the strong solution is unique.
Proof By the assumption p n,K = p n we can apply Lem. 4 with the construction of the interpolation I h and Cor. 1, so that (u h ) h∈H is uniformly bounded by By (32) and the definition of h (with g = 0), this also implies that is uniformly bounded for h ∈ H, so that together with the asymptotic consistency of the material parameters we obtain with a constant C f ,u0 > 0 depending on the data Then we obtain for all v ∈ V h u, using dual consistency (31) for the last step.This extends to H 1 0 (Q; R m ), and by the assumption i.e., for the limit u the weak derivative Lu = f in L 2 (Q; R m ) exists.This extends to initial and boundary data.Therefore, let exists, and we get again by (31) , and thus u is indeed a strong solution with homogeneous boundary conditions at (0, T ) × ∂Ω.
Next, we show that the weak limit is unique.Therefore, select another subsequence H 1 ⊂ H with 0 ∈ H 1 and with a weak limit ũ Then, we also obtain ũ(0) = u 0 and (A n ũ) k = 0 for k = 1, . . ., m.A sequence (e h ) h∈H with e h ∈ V h exists such that lim h∈H e h = u − ũ, and we get so that u = ũ.This shows that the weak limit is unique, so that the full sequence is converging, i.e., lim h∈H u h = u.
The same argument applies to all strong solutions, i.e., u is the unique strong solution of (3).

Remark 9
The result extends to inhomogeneous boundary data In particular, the regularity result that the limit of the DG approximations is a strong solution requires sufficient regularity of the boundary data.

Convergence in the DG norm
We adapt the convergence result for the DG norm (41) in [5, Thm.2.37] to our setting.
Theorem 3 Assume that the strong solution of (3) is sufficiently smooth satisfying u ∈ H s (Q; R m ) with s ≥ 1 and s ≤ min n,K {p n,K , q n,K } + 1.Then, the error for the discrete solution u h ∈ V h of (36) is bounded by with C > 0 depending on the mesh regularity, the polynomial degrees in V h , and the material parameters.
Proof Since we assume for the solution u ∈ H Thus we obtain for the discrete solution u h ∈ V h Galerkin orthogonality up to data error By the trace estimate ( 13) we obtain w h Q , so that the consistency term can by bounded by For all v h ∈ V h this yields the estimate, using Thm. 1 and continuity of the bilinear form b h (•, •) in the DG norms and constants C 4 , C 5 depending on the mesh regularity and the polynomial degrees in V h .Using s ≤ min{p, q} + 1, Then, the result follows from interpolation estimates using [5, Lem.1.59] and This recovers the convergence result [2, Prop.6.5] for the DG semi-norm (41).
Corollary 2 Assume that the strong solution of (3) is sufficiently smooth Then, the error for the discrete solution u h ∈ V h of (36) is bounded in every time step by with C > 0 depending on the mesh regularity, the polynomial degree, and the material parameters.
For the proof Thm. 3 is applied with T = t n ; then, the assertion directly follows from 1  2 Remark 10 If M ∈ L ∞ (Ω; R m×m sym ) is smooth, the consistency term can be estimated by discontinuous and if the jumps of the material parameters are not resolved by the mesh, the consistency error can be estimated in case of higher regularity of the solution: if ∂ t u ∈ L 2 (0, T ; L q (Ω; R m )) with q > 2, we obtain Remark 11 For the continuous solution the energy is conserved, i.e., From Cor. 2 we obtain energy conservation in the limit

Remark 12
The constants in Thm. 1 and 3 depend on the mesh and polynomial degrees p.For triangulations and a quasi-uniform distribution of p it is known that C inv ∼ p 2 , C tr ∼ p [23,Thm. 4.7].Estimates of quasi-interpolations are considered in [20,Thm. 3.1] where it is shown that the classical Clément interpolation estimate holds with h replaced by h/p.

Error control
For the error u − u h in the DG semi-norm we obtain from ( 19) and ( 21) and in the DG norm Up to the error u h −u at final time T in (47) and the parameter approximation error M − M h in (48), this can be evaluated explicitly by the residual error given by the local contributions Lemma 5 Let u ∈ L 2 (Q; R m ) be the weak solution of (2) and u h ∈ V h the discrete solution of (36).
Then, if u is a strong solution, the error in the DG norm is bounded by .

Numerical experiments
The convergence estimates in the DG norm are illustrated by numerical experiments for acoustics (5) for cases where the exact solution is know which is then used for Dirichlet boundary conditions.The results for uniform refinement are compared with a simple adaptive strategy by increasing the polynomial degree for η res,R ≥ θ 1 max R η res,R and decreasing the polynomial degree for η res,R ≤ θ 0 max R η res,R , see [6] for details.In addition, we consider an example motivated from the application to seismic imaging where the exact solution is not known, and the convergence is demonstrated with respect to the residual error indicator.
Experiment 1 We test the convergence of the solution in Q = (0, 1) × (0, 1) 2 and f = 0 with smooth initial value and piecewise constant material so that the impedance is constant across the interface.We start with Then, the solution is given by u Case a) If the material interface is resolved by the mesh (M = M h ), we observe for linear approximations in space and time on uniformly refined meshes the expected convergence rate in the DG norm (Fig. 1).For this configuration also the dual problem is smooth which results in better convergence rates for the L 2 error, in particular in the adaptive case.Although the material interface cannot be resolved by the mesh, the solution is sufficiently smooth so that the approximation error of the material data M h − M can be estimated by Rem. 10.We observe nearly optimal convergence in the DG norm, but now the L 2 convergence gets worse in comparison with the first case.
In both cases, the convergence of u(T ) − u h (T ) in L 2 is faster than the convergence in the DG norm, and the residual error indicator yields results close to the error in the DG norm; this confirms the estimate in Lem. 5. We observe that adaptivity provides better solutions with a substantial reduction of the required problem size dim V h to achieve a certain accuracy.Therefore a single adaptive step is sufficient, where the polynomial degree in space and time is increased for η res,R ≥ ϑ 1 max R ∈R h η res,R and decreased for η res,R ≤ ϑ 0 max R ∈R h η res,R , depending on ϑ 1 > ϑ 0 > 0. Note that this results in a different refinement pattern in every time interval, and a simple refinement in space is not sufficient for a strong reduction of the required degrees of unknowns.Here, we select ϑ 1 = 0.3 and ϑ 0 = 0.02, and in the figures for the adaptive results the mesh size is logarithmically interpolated depending on the degrees of freedom.
Experiment 2 At next, we test the convergence of a Riemann problem in Q = (0, 1/2) × (−1, 1) × (0, 1) with f = 0, where the solution is given by Then, Lu = 0, so that u is a strong solution, and since the condition in Rem. 9 applies, we obtain convergence in the limit by Thm. 2. On the other hand, the solution is piecewise discontinuous, so that the smoothness assumption in Thm. 3 is not satisfied.We also observe convergence, cf.Fig. 3, but with a reduced rate O(h 1/3 ).In particular, the rate is not improved for the L 2 error, and simple adaptivity is not sufficient to increase the efficiency.
Here, the solution is not smooth, and the results do not improve if the material parameters are aligned with the mesh.Moreover, further tests show that the convergence order of approximately O(h 0.4 ) in the DG norm cannot be improved by adaptivity, which indicates that without sufficient regularity and jumps along the characteristics the DG norm is not appropriate for a qualitative convergence analysis, as it is possible for point singularities, see [2].Then, the convergence analysis requires high regularity in weighted Sobolev spaces.The configuration, the distribution of the the piecewise constant parameters and κ, and the parallel solution framework in M++ are described in detail in [8].Since in this application only the evaluation in a small measurement region (4.75, 7.25)×(0, 0.4) ⊂ Ω is of interest, the space-time domain can be truncated, see [10,Lem. 2].Here the convergence is only tested by evaluating the residual error indicator on uniformly refined meshes and for one and two p-adaptive steps with θ 0 = 0.01 and θ 1 = 0.1.Since all data are aligned with the mesh but discontinuous, the regularity of the solution is limited.We observe approximately linear convergence with respect to the estimate of the DG norm, and again we observe improved convergence by space-time adaptivity, cf.Fig. 4.

Conclusion and Outlook
The convergence analysis in the DG norm only assumes regularity of the spacetime solution u in H 1 (Q; R m ); this implies regularity of the solution u(t n ) at all time steps in H 1/2 (Ω; R m ).This clearyly extends convergence results with respect to the graph norm, where the analysis requires higher regularity.Moreover, the simple residual error indicator yields estimates very close to the error in the DG norm.On the other hand, for discontinuous Riemann problems we can prove only convergence in the limit, and the numerical experiments demonstrate that we obtain convergence in L 2 but with a reduced rate, which can be improved by adaptivity in L 2 but not in the DG norm.
All our estimates rely on a Hilbert space setting.This may be not appropriate for hyperbolic systems, and numerical tests demonstrate better convergence rates in L 1 (Q; R m ), but a corresponding analysis remains an open problem.

Fig. 4 :
Fig. 4: Convergence test for a forward problem in seismic imaging in a truncated space-time domain.