A stable discontinuous Galerkin method for the perfectly matched layer for elastodynamics in first order form

Abstract

We present a stable discontinuous Galerkin (DG) method with a perfectly matched layer (PML) for three and two space dimensional linear elastodynamics, in velocity-stress formulation, subject to well-posed linear boundary conditions. First, we consider the elastodynamics equation, in a cuboidal domain, and derive an unsplit PML truncating the domain using complex coordinate stretching. The hyperbolic structure of the underlying system enables the construction of continuous energy estimates, in the time domain for the elastic wave equation, and in the Laplace space for a sequence of PML model problems, with variations in one, two and three space dimensions, respectively. They correspond to PMLs normal to boundary faces, along edges and in corners. Second, we develop a DG numerical method for the linear elastodynamics equation using physically motivated numerical flux and penalty parameters, which are compatible with all well-posed, internal and external, boundary conditions. When the PML damping vanishes in all directions, by construction, our choice of penalty parameters yield an upwind scheme and a discrete energy estimate analogous to the continuous energy estimate. Third, to ensure numerical stability of the discretization when PML damping is present, it is necessary to extend the numerical DG fluxes, and the numerical inter-element and boundary procedures, to the PML auxiliary differential equations. This is crucial for deriving discrete energy estimates analogous to the continuous energy estimates. Numerical solutions are evolved in time using the high order arbitrary derivative (ADER) time stepping scheme of the same order of accuracy with the spatial discretization. By combining the DG spatial approximation with the high order ADER time stepping scheme and the accuracy of the PML we obtain an arbitrarily high-order accurate wave propagation solver in the time domain. Numerical experiments are presented in two and three space dimensions corroborating the theoretical results.

Appendices

Other receivers

In this section we present the seismograms for the remaining receivers, namely: Receiver 1, 2, 3, 7 and 8. In Fig. 12 we display the results for the HHS problem. Figure 13 shows the result for LOH1 problem with $N = 3$ degree polynomial approximation, and Fig. 14 shows the result for LOH1 problem with $N = 5$ degree polynomial approximation. Note the effective absorption properties of the PML, and the spectral accuracy of the DG method.

Fig. 12

The homogeneous half space problem

Fig. 13

The LOH1 benchmark problem with degree $N = 3$ polynomial approximation

Fig. 14

The LOH1 benchmark problem with degree $N = 5$ polynomial approximation

Source terms

We will derive a representation of the modified source terms in the physical space. Let $t\ge 0$ denote time, and f(t) and g(t) are real functions of exponential order. We denote their corresponding Laplace transforms ${\widetilde{f}}(s) = {\mathscr {L}}\left( f(t)\right)$, ${\widetilde{g}}(s) = {\mathscr {L}}\left( g(t)\right)$. We also define the convolution operator $*$,

$$\begin{aligned} g(t)*f(t) = \int _{0}^{t} g(t-\tau )f(\tau ) d\tau . \end{aligned}$$

Note that ${\mathscr {L}}\left( g(t)*f(t)\right) = {\widetilde{f}}(s){\widetilde{g}}(s)$. To succeed we will use the Dirac delta distribution $\delta (t)$

$$\begin{aligned} \int _{0}^{+\infty } \delta (t) dt = 1, \end{aligned}$$

with

$$\begin{aligned} \delta (t)*f(t) = f(t), \quad \delta ^{\prime }(t)*f(t) = f^{\prime }(t) + \delta (t)f(0). \end{aligned}$$

We also recall

$$\begin{aligned} {\mathscr {L}}\left( \delta ^{n}(t)\right) = s^n, \quad n = 0, 1, 2, \cdots , \end{aligned}$$

where $\delta ^{n}(t)$ is the n-th derivative of $\delta (t)$.

We have

$$\begin{aligned} {\mathscr {L}}\left( \delta ^{\prime }(t)*f(t) \right) = s{\widetilde{f}}(s). \end{aligned}$$

Consider the modified source terms in the Laplace domain

$$\begin{aligned} \widetilde{{\mathbf {F}}} = {\mathbf {P}}\left( S_x \widetilde{{\mathbf {F}}}_{Q} - \sum _{\xi = x,y, z} \frac{d_{\xi }S_x }{s + \alpha _{\xi } + d_{\xi }} \widetilde{{\mathbf {F}}}_{w_\xi }\right) . \end{aligned}$$

(117)

If $d_{\xi } = d_x = d \ge 0$ and $\alpha _{\xi } = \alpha _{x} = \alpha \ge 0$, then

$$\begin{aligned} S_x \widetilde{{\mathbf {F}}}_{Q}= \left( 1 + \frac{d}{s + \alpha } \right) \widetilde{{\mathbf {F}}}_{Q}, \quad \frac{S_x }{s + \alpha _{\xi } + d_{\xi }} \widetilde{{\mathbf {F}}}_{w_\xi } = \frac{1 }{s + \alpha } \widetilde{{\mathbf {F}}}_{w_\xi }, \end{aligned}$$

(118)

and

$$\begin{aligned}&{\mathscr {L}}^{-1}\left( S_x \widetilde{{\mathbf {F}}}_{Q}\right) = \left( \delta (t) + {d} e^{-\alpha t} \right) * {{\mathbf {F}}}_{Q}\left( x,y,z, t\right) , \nonumber \\&{\mathscr {L}}^{-1}\left( \frac{1 }{s + \alpha } \widetilde{{\mathbf {F}}}_{w_\xi }\right) = e^{-\alpha t} * {{\mathbf {F}}}_{w_\xi } \left( x,y,z, t\right) , \quad \text {Re}\{s\} > - \alpha . \end{aligned}$$

(119)

Proof of Theorem 2

Proof

As before, multiply (53) with $\widetilde{{\mathbf {Q}}}^{\dagger }$ from the left and integrate over the whole spatial domain. Integration-by-parts gives

$$\begin{aligned} \begin{aligned}&\int _{\varOmega }\left( \left( sS_x\right) \widetilde{{\mathbf {Q}}}^{\dagger } {\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}} \right) dxdy dz \\&\quad = \sum _{\xi = x, y}\frac{1}{2}\int _{\varOmega }\left( \left[ \widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi }\frac{\partial {\widetilde{{\mathbf {Q}}}}}{\partial \xi }\right) -\left( {\mathbf {A}}_{\xi }\frac{\partial {\widetilde{{\mathbf {Q}}}}}{\partial \xi }\right) ^\dagger {\mathbf {Q}}\right] \right) dxdydz \\&\qquad + \sum _{\xi = x, y}\int _{{\widetilde{\varGamma }}}\left( \left[ \widetilde{{\mathbf {Q}}}^{\dagger }{A}_\xi \widetilde{{\mathbf {Q}}}\right] \Big |_{-1}^{1}\right) \frac{dxdydz}{d\xi } + \int _{\varOmega }\left( \widetilde{{\mathbf {Q}}}^{\dagger } {\mathbf {P}}^{-1}\widetilde{{\mathbf {F}}} \right) dxdydz. \end{aligned} \end{aligned}$$

(120)

Adding (120) to its complex conjugate, the spatial derivative terms vanish, yielding

$$\begin{aligned} \begin{aligned} \Vert \sqrt{{\mathrm{Re}}(s S_x)}\widetilde{{\mathbf {Q}}} (\cdot ,\cdot ,\cdot , s)\Vert _{P}^2&= \frac{1}{2}\left[ \left( \widetilde{{\mathbf {Q}}},\widetilde{{\mathbf {F}}} \right) _{{P}} + \left( \widetilde{{\mathbf {F}}},\widetilde{{\mathbf {Q}}} \right) _{{P}}\right] \\&\quad + \sum _{\xi = x, y}\int _{{\widetilde{\varGamma }}}\left( \left[ \widetilde{{\mathbf {Q}}}^{\dagger }{A}_x \widetilde{{\mathbf {Q}}}\right] \Big |_{-1}^{1}\right) \frac{dxdydz}{d\xi } . \end{aligned} \end{aligned}$$

(121)

Using Cauchy-Schwarz inequality and the fact (41) in the right hand side of (121) completes the proof. $\square$

Proof of Theorem 3

Proof

As above, multiply (55) with $\widetilde{{\mathbf {Q}}}^{\dagger }$ from the left and integrate over the whole spatial domain. Integration-by-parts gives

$$\begin{aligned}&\int _{\varOmega }\left( \left( sS_x\right) \widetilde{{\mathbf {Q}}}^{\dagger } {\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}} \right) dxdy dz \nonumber \\&\quad = \sum _{\xi = x, y,z}\frac{1}{2}\int _{\varOmega }\left( \left[ \widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi } \frac{\partial {\widetilde{{\mathbf {Q}}}}}{\partial \xi }\right) -\left( {\mathbf {A}}_{\xi }\frac{\partial {\widetilde{{\mathbf {Q}}}}}{\partial \xi }\right) ^\dagger {\mathbf {Q}}\right] \right) dxdydz \nonumber \\&\qquad + \sum _{\xi = x, y,z}\int _{{\widetilde{\varGamma }}} \left( \left[ \widetilde{{\mathbf {Q}}}^{\dagger }{A}_\xi \widetilde{{\mathbf {Q}}}\right] \Big |_{-1}^{1}\right) \frac{dxdydz}{d\xi } + \int _{\varOmega }\left( \widetilde{{\mathbf {Q}}}^{\dagger } {\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} \right) dxdydz. \end{aligned}$$

(122)

Adding the complex conjugate of the product, the spatial derivative terms vanish, yielding

$$\begin{aligned} \begin{aligned} \Vert \sqrt{{\mathrm{Re}}(s S_x)}\widetilde{{\mathbf {Q}}} (\cdot ,\cdot ,\cdot , s)\Vert _{P}^2&= \frac{1}{2}\left[ \left( \widetilde{{\mathbf {Q}}},\widetilde{{\mathbf {F}}} \right) _{{P}} + \left( \widetilde{{\mathbf {F}}},\widetilde{{\mathbf {Q}}} \right) _{{P}}\right] \\&\quad + \sum _{\xi = x, y, z}\int _{{\widetilde{\varGamma }}}\left( \left[ \widetilde{{\mathbf {Q}}}^{\dagger }{A}_\xi \widetilde{{\mathbf {Q}}}\right] \Big |_{-1}^{1}\right) \frac{dxdydz}{d\xi } , \quad \text {Re}\{s\} > 0. \end{aligned}\nonumber \\ \end{aligned}$$

(123)

Using Cauchy-Schwarz inequality and the fact (41) in the right hand side of (123) completes the proof. $\square$

Proof of Theorem 8

We introduce the boundary terms

$$\begin{aligned} \mathscr {BT}_s\left( \widehat{{\widetilde{v}}}_{\eta }^{-}, \widehat{{\widetilde{T}}}_{\eta }^{-}\right)&= \frac{\varDelta {x}}{2} \frac{\varDelta {z}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{r = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{r = -1}\right) _{i k} h_{i} h_{k} \\&\quad - \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{q = -1}\right) _{i k} h_{i} h_{k} \le 0,\\ \mathscr {BT}_s\left( \widehat{{\widetilde{v}}}_{\eta }^{+}, \widehat{{\widetilde{T}}}_{\eta }^{+}\right)&= \frac{\varDelta {x}}{2} \frac{\varDelta {z}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{r = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{r = -1}\right) _{i k} h_{i} h_{k} \\&\quad + \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{q = 1}\right) _{i k} h_{i} h_{k} \le 0, \end{aligned}$$

the interface term

$$\begin{aligned} \mathscr {IT}_s\left( \widehat{{\widetilde{v}}}^{\pm }, \widehat{{\widetilde{T}}}^{\pm } \right) = - \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( {\widehat{{\widetilde{T}}}}_\eta \llbracket {{\widehat{{\widetilde{v}}}}_\eta \rrbracket }\right) _{i k} h_{i} h_{k} \equiv 0, \end{aligned}$$

and the fluctuation term

$$\begin{aligned} {\mathscr {F}}_{luc}\left( {\widetilde{G}}, Z\right)&= - \sum _{\xi = x, y}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }}\\&\quad \sum _{\eta = x, y, z}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\left( \left( \frac{1}{Z_{\eta }}|{\widetilde{G}}_\eta |^2\right) \Big |_{ -1} + \left( \frac{1}{Z_{\eta }}|{\widetilde{G}}_\eta |^2\right) \Big |_{ 1} \right) _{i, k} \le 0. \end{aligned}$$

Proof

From the left, multiply (100) with $\widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}$, we have

$$\begin{aligned} \begin{aligned}&\left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}}\\&\quad = \sum _{\xi = x, y}\frac{1}{2}\widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi }\left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) - \left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) ^T{\mathbf {A}}_{\xi } \right) \widetilde{{\mathbf {Q}}} + \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} \\&\qquad +\widetilde{{\mathbf {Q}}}^\dagger {\mathbf {H}}{\mathbf {H}}_{\xi }^{-1} \left( \frac{1}{2}{\mathbf {A}}_{x}\left( {\mathbf {B}}_{\xi }\left( 1,1\right) - {\mathbf {B}}_{\xi }\left( -1,-1\right) \right) \widetilde{{\mathbf {Q}}} -\left( {{\mathbf {e}}_{\xi }(-1)} \widetilde{\mathbf {FL}}_{\xi } + {{\mathbf {e}}_{\xi }(1)} \widetilde{\mathbf {FR}}_{\xi } \right) \right) \end{aligned}\nonumber \\ \end{aligned}$$

(124)

Using the identity (95) gives

$$\begin{aligned} \left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}}&= \sum _{\xi = x, y}\frac{1}{2}\widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi }\left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) - \left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) ^T{\mathbf {A}}_{\xi } \right) \widetilde{{\mathbf {Q}}} + \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} \nonumber \\&\quad - \sum _{\xi = x, y}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }} \sum _{i=1}^{P+1}\sum _{k= 1}^{P+1}\left( \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 - \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{ 1} \right. \nonumber \\&\quad \left. + \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 + \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{-1}\right) _{ik} h_ih_k \end{aligned}$$

(125)

Next we add (125) to its complex conjugate, and the spatial derivative terms vanish, we have

$$\begin{aligned} {\mathrm{Re}}\left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}}= & {} \frac{1}{2}\left( \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} + \widetilde{{\mathbf {F}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}} \right) \nonumber \\&- \sum _{\xi = x, y}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }} \sum _{i=1}^{P+1}\sum _{k= 1}^{P+1}\left( \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 - \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{ 1} \right. \nonumber \\&\left. + \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 + \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{-1}\right) _{ik} h_ih_k \end{aligned}$$

(126)

Using Cauchy-Schwarz inequality for the source term in (126) and collecting contributions from both sides of the elements completes the proof. $\square$

Proof of Theorem 9

We introduce the boundary terms

$$\begin{aligned} \mathscr {BT}_s\left( \widehat{{\widetilde{v}}}_{\eta }^{-}, \widehat{{\widetilde{T}}}_{\eta }^{-}\right)&= \frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{s = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{s = -1}\right) _{i k} h_{i} h_{k} \\&\quad + \frac{\varDelta {x}}{2} \frac{\varDelta {z}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{r = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{r = -1}\right) _{i k} h_{i} h_{k}\\&\quad - \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{-*} \widehat{{\widetilde{v}}}_\eta ^{-} \right) \Big |_{q = -1}\right) _{i k} h_{i} h_{k} \le 0,\\ \mathscr {BT}_s\left( \widehat{{\widetilde{v}}}_{\eta }^{+}, \widehat{{\widetilde{T}}}_{\eta }^{+}\right)&= \frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{s = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{s = -1}\right) _{i k} h_{i} h_{k} \\&\quad + \frac{\varDelta {x}}{2} \frac{\varDelta {z}}{2} \sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{r = 1} - \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{r = -1}\right) _{i k} h_{i} h_{k} \\&\quad + \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( \left( \widehat{{\widetilde{T}}}_\eta ^{+*} \widehat{{\widetilde{v}}}_\eta ^{+} \right) \Big |_{q = 1}\right) _{i k} h_{i} h_{k} \le 0, \end{aligned}$$

the interface term

$$\begin{aligned} \mathscr {IT}_s\left( \widehat{{\widetilde{v}}}^{\pm }, \widehat{{\widetilde{T}}}^{\pm } \right) = - \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\sum _{\eta = x,y,z}\left( {\widehat{{\widetilde{T}}}}_\eta \llbracket {{\widehat{{\widetilde{v}}}}_\eta \rrbracket }\right) _{i k} h_{i} h_{k} \equiv 0, \end{aligned}$$

and the fluctuation term

$$\begin{aligned} {\mathscr {F}}_{luc}\left( {\widetilde{G}}, Z\right)&= - \sum _{\xi = x, y, z}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }}\\&\quad \sum _{\eta = x, y, z}\sum _{i = 1}^{P+1}\sum _{k = 1}^{P+1}\left( \left( \frac{1}{Z_{\eta }}|{\widetilde{G}}_\eta |^2\right) \Big |_{ -1} + \left( \frac{1}{Z_{\eta }}|{\widetilde{G}}_\eta |^2\right) \Big |_{ 1} \right) _{i, k} \le 0. \end{aligned}$$

Proof

From the left, multiply (102) with $\widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}$, we have

$$\begin{aligned} \begin{aligned}&\left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}} \\&\quad = \sum _{\xi = x, y, z}\frac{1}{2}\widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi }\left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) - \left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) ^T{\mathbf {A}}_{\xi } \right) \widetilde{{\mathbf {Q}}} + \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} \\&\qquad + \widetilde{{\mathbf {Q}}}^\dagger {\mathbf {H}}{\mathbf {H}}_{\xi }^{-1} \left( \frac{1}{2}{\mathbf {A}}_{x}\left( {\mathbf {B}}_{\xi }\left( 1,1\right) - {\mathbf {B}}_{\xi }\left( -1,-1\right) \right) \widetilde{{\mathbf {Q}}} -\left( {{\mathbf {e}}_{\xi }(-1)} \widetilde{\mathbf {FL}}_{\xi } + {{\mathbf {e}}_{\xi }(1)} \widetilde{\mathbf {FR}}_{\xi } \right) \right) \end{aligned}\nonumber \\ \end{aligned}$$

(127)

Using the identity (95) gives

$$\begin{aligned} \left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}}= & {} \sum _{\xi = x, y, z}\frac{1}{2}\widetilde{{\mathbf {Q}}}^\dagger \left( {\mathbf {A}}_{\xi }\left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) - \left( {\mathbf {H}}{\mathbf {D}}_{\xi }\right) ^T{\mathbf {A}}_{\xi } \right) \widetilde{{\mathbf {Q}}} + \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} \nonumber \\&- \sum _{\xi = x, y}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }} \sum _{i=1}^{P+1}\sum _{k= 1}^{P+1}\left( \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 - \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{ 1} \right. \nonumber \\&\quad \left. + \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 + \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{-1}\right) _{ik} h_ih_k \end{aligned}$$

(128)

Next we add (128) to its complex conjugate, and the spatial derivative terms vanish, we have

$$\begin{aligned} {\mathrm{Re}}\left( sS_x\right) \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}}= & {} \frac{1}{2}\left( \widetilde{{\mathbf {Q}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {F}}} + \widetilde{{\mathbf {F}}}^\dagger {{\mathbf {H}}}{\mathbf {P}}^{-1} \widetilde{{\mathbf {Q}}} \right) \nonumber \\&- \sum _{\xi = x, y}\frac{\varDelta {x}}{2} \frac{\varDelta {y}}{2} \frac{\varDelta {z}}{2} \frac{2}{\varDelta {\xi }} \sum _{i=1}^{P+1}\sum _{k= 1}^{P+1}\left( \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 - \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{ 1} \right. \nonumber \\&\left. + \sum _{\eta = x,y,z} \left( \frac{1}{Z_\eta }|{\widetilde{G}}_\eta |^2 + \widehat{{\widetilde{T}}}_\eta ^* \widehat{{\widetilde{v}}}_\eta \right) \Big |_{-1}\right) _{ik} h_ih_k \end{aligned}$$

(129)

Using Cauchy-Schwarz inequality for the source term in (129) and collecting contributions from both sides of the elements completes the proof. $\square$

The Arbitrary DERivative (ADER) time integration

In this section we will summarize the ADER time-stepping scheme. For more elaborate discussions, we refer the reader to [8, 10, 38]. The DG semi-discrete approximation (84)–(85) can be written as a system first order ODEs

$$\begin{aligned} \begin{aligned} \frac{d \bar{{\mathbf {Q}}} }{dt} = \underbrace{D\bar{{\mathbf {Q}}}}_{\text {PDE}} + \underbrace{F\bar{{\mathbf {Q}}}}_{\text {Num. flux} \rightarrow 0}, \end{aligned} \end{aligned}$$

(130)

where $D\bar{{\mathbf {Q}}}$ emanates from the PDE, is the discrete spatial operator involving derivatives and lower order PML terms, and $F\bar{{\mathbf {Q}}}$ is the numerical flux fluctuation term, incorporating the boundary and interface conditions.

The numerical flux fluctuation is a very small term, $F\bar{{\mathbf {Q}}} \approx 0$, and will vanish $F\bar{{\mathbf {Q}}} \rightarrow 0$ in the limit of mesh refinement $\varDelta {t} \rightarrow 0$.

We now introduce the discrete time variables $t_k \le t \le t_{k+1}$, $\varDelta {t}_k = t_{k+1}- t_{k}$, and the pseudo time variable $\tau = t -t_k$ such that $0 \le \tau \le \varDelta {t}_k$, and ${d}/{d\tau } = {d}/{dt}$. Going from the current time level $\tau =0$ to the next time level $\tau = \varDelta {t}_k$, we integrate the ODE (130), exactly having

$$\begin{aligned} \bar{{\mathbf {Q}}}(\varDelta {t})&= \bar{{\mathbf {Q}}}(0) + \int _0^{\varDelta {t}_k}{D\bar{{\mathbf {Q}}}} d\tau + \int _0^{\varDelta {t}_k}{F\bar{{\mathbf {Q}}}} d\tau , \quad \nonumber \\&= \bar{{\mathbf {Q}}}(0) + D\int _0^{\varDelta {t}_k}{\bar{{\mathbf {Q}}}} d\tau + F\int _0^{\varDelta {t}_k}{\bar{{\mathbf {Q}}}} d\tau , \end{aligned}$$

(131)

where the second equality follows from linearity. If we can evaluate the integrals $\int _0^{\varDelta {t}_k}{\bar{{\mathbf {Q}}}} d\tau$ in (131) exactly, then the time integration in (131) is exact. However, exact time integration is possible only in the most trivial case where the right hand side of (130) vanish identically for all components. Now, we will make an important approximation. We assume that the time step $\varDelta {t}$ is sufficiently small, such that

$$\begin{aligned} \frac{d\bar{{\mathbf {Q}}}(\tau )}{d\tau } \approx D\bar{{\mathbf {Q}}}, \quad F\bar{{\mathbf {Q}}} \approx 0, \end{aligned}$$

(132)

are reasonable approximations. Next we construct the predictor, $\widetilde{\bar{{\mathbf {Q}}}}(\tau )$, by Taylor expansions of the solution around $\tau = 0$ and replace the time derivatives with spatial operator in (132), we have

$$\begin{aligned} \widetilde{\bar{{\mathbf {Q}}}}(\tau ) = \bar{{\mathbf {Q}}}(0) + \tau \frac{d\bar{{\mathbf {Q}}}(0)}{d\tau } + \frac{\tau ^2}{2} \frac{d^2\bar{{\mathbf {Q}}}(0)}{d\tau ^2} + ... \approx \sum _{m = 0}^{P}\frac{\tau ^m}{m !} D^m \bar{{\mathbf {Q}}}(0), \end{aligned}$$

(133)

where P is the polynomial degree used in the spatial approximation. We can now approximate the integrals in (131) using the predictor. The result of this integration is called the time average, $\bar{\bar{{\mathbf {Q}}}}(0)$,

$$\begin{aligned} \int _0^{\varDelta {t}_k}{\bar{{\mathbf {Q}}}} d\tau \approx { \bar{\bar{{\mathbf {Q}}}}(0) = \int _{0}^{\varDelta {t}_k} \widetilde{\bar{{\mathbf {Q}}}}(\tau ) d\tau = \sum _{m = 0}^{P}\frac{\varDelta {t}_k^{(m+1)}}{(m +1) !} D^m \bar{{\mathbf {Q}}}(0)}. \end{aligned}$$

(134)

By replacing the integrals in (131) with the time average ${\bar{\bar{{\mathbf {Q}}}}(0)}$, we derive a high order accurate, explicit, one-step, time integration scheme

$$\begin{aligned} \bar{{\mathbf {Q}}}(\varDelta {t})&= \bar{{\mathbf {Q}}}(0) + D\bar{\bar{{\mathbf {Q}}}}(0) + F\bar{\bar{{\mathbf {Q}}}}(0). \end{aligned}$$

(135)

Note that the numerical flux fluctuation, $F\bar{\bar{{\mathbf {Q}}}}(0)$, is evaluated only once for any order of approximation. This is opposed to Runge–Kutta methods or standard Taylor series methods where the numerical flux fluctuation is included in the spatial operator, $D^m \rightarrow (D + {F})^m$, to approximate higher time derivatives in the Taylor series terms. For the ADER scheme, the fact that the numerical flux fluctuation is evaluated only once for any order implies that most of the computations are performed within the element to compute the predictor in (133) and the time average in (134). This tremendously decreases computational cost of high performance computing applications, since we can design efficient communication avoiding parallel algorithms. Since the predictor $\widetilde{\bar{{\mathbf {Q}}}}(\tau )$ is defined in the entire time interval $0 \le \tau \le \varDelta {t}_k$, the ADER scheme, (131)–(135), is also easily amenable to local time-stepping methods. When the ADER time stepping scheme is combined with the DG spatial approximation the fully discrete scheme is often referred as the ADERDG method [10]. For a DG polynomial approximation of degree P, a stable ADERDG method is $(P+1)$th order accurate in both space and time.

Hat-variables

We will now reformulate the boundary condition (21) and interface condition (24) by introducing transformed (hat-) variables so that we can simultaneously construct (numerical) boundary/interface data for particle velocities and tractions. The hat-variables encode the solution of the IBVP on the boundary/interface. The hat-variables will be constructed such that they preserve the amplitude of the outgoing characteristics and exactly satisfy the physical boundary conditions [17]. To be more specific, the hat-variables are solutions of the Riemann problem constrained against physical boundary conditions (21), and the interface condition (24). We refer the reader to [17, 18] for more detailed discussions. Once the hat-variables are available, we construct physics based numerical flux fluctuations by penalizing data against the incoming characteristics (20) at the element faces.

1.1 Boundary data

For $Z_\eta > 0$, we define the characteristics

$$\begin{aligned} q_\eta = \frac{1}{2}\left( Z_\eta v_\eta + T_\eta \right) , \quad p_\eta = \frac{1}{2}\left( Z_\eta v_\eta - T_\eta \right) , \quad \eta = x, y, z. \end{aligned}$$

(136)

Here, $q_\eta$ are the left going characteristics and $p_\eta$ are the right going characteristics. We will construct boundary data which satisfy the physical boundary conditions (21) exactly and preserve the amplitude of the outgoing characteristics $q_\eta$ at $\xi = -1$, and $p_\eta$ at $\xi = 1$.

To begin, define the hat-variables preserving the amplitude of outgoing characteristics

$$\begin{aligned}&{q}_\eta \left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = {q}_\eta \left( {v}_\eta , {T}_\eta , Z_\eta \right) , \quad \text {at} \quad \xi = -1, \nonumber \\&\quad \text {and} \quad {p}_\eta \left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = {p}_\eta \left( {v}_\eta , {T}_\eta , Z_\eta \right) , \quad \text {at} \quad \xi = 1. \end{aligned}$$

(137)

Since hat-variables also satisfy the physical boundary condition (21), we must have

$$\begin{aligned}&\frac{Z_{\eta }}{2}\left( {1-\gamma _\eta }\right) {\widehat{v}}_\eta -\frac{1+\gamma _\eta }{2} {\widehat{T}}_\eta = 0, \quad \text {at} \quad \xi = -1, \nonumber \\&\quad \text {and} \quad \frac{Z_{\eta }}{2} \left( {1-\gamma _\eta }\right) {\widehat{v}} _\eta + \frac{1+\gamma _\eta }{2}{\widehat{T}}_\eta = 0, \quad \text {at} \quad \xi = 1. \end{aligned}$$

(138)

The algebraic problem for the hat-variables, defined by equations (137) and (138), has a unique solution, namely

$$\begin{aligned}&{\widehat{v}}_\eta = \frac{(1+\gamma _\eta )}{Z_{\eta }}q_\eta , \quad {\widehat{T}}_\eta = {(1-\gamma _\eta )}q_\eta , \quad \text {at} \quad \xi = -1, \nonumber \\&\quad \text {and} \quad {\widehat{v}}_\eta = \frac{(1+\gamma _\eta )}{Z_{\eta }}p_\eta , \quad {\widehat{T}}_\eta = {-(1-\gamma _\eta )}p_\eta , \quad \text {at} \quad \xi = 1. \end{aligned}$$

(139)

The expressions in (139) define a rule to update particle velocities and tractions on the physical boundaries $\xi = -1, 1$. That is

$$\begin{aligned} v_\eta = {\widehat{v}}_\eta , \quad {T}_\eta = {\widehat{T}}_\eta , \quad \text {at} \quad \xi = -1, \quad \text {and} \quad v_\eta = {\widehat{v}}_\eta , \quad {T}_\eta = {\widehat{T}}_\eta , \quad \text {at} \quad \xi = 1. \end{aligned}$$

(140)

By construction, the hat-variables ${\widehat{v}}_\eta , {\widehat{T}}_\eta$, satisfy the following algebraic identities:

$$\begin{aligned}&{q}_\eta \left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = {q}_\eta \left( {v}_\eta , {T}_\eta , Z_\eta \right) , \quad \text {at} \quad \xi = -1, \nonumber \\&\quad {p}_\eta \left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = {p}_\eta \left( {v}_\eta , {T}_\eta , Z_\eta \right) , \quad \text {at} \quad \xi = 1,\end{aligned}$$

(141a)

$$\begin{aligned}&{q}_\eta ^2\left( {v}_\eta , {T}_\eta , Z_\eta \right) -{p}_\eta ^2\left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = Z_\eta {\widehat{T}}_\eta {\widehat{v}}_\eta , \quad \text {at} \quad \xi = -1, \nonumber \\&\quad {p}_\eta ^2\left( {v}_\eta , {T}_\eta , Z_\eta \right) -{q}_\eta ^2\left( {\widehat{v}}_\eta , {\widehat{T}}_\eta , Z_\eta \right) = -Z_\eta {\widehat{T}}_\eta {\widehat{v}}_\eta , \quad \text {at} \quad \xi = 1,\end{aligned}$$

(141b)

$$\begin{aligned}&{\widehat{T}}_\eta {\widehat{v}}_\eta = \frac{1-\gamma _\eta ^2}{Z_{\eta }}{q}_\eta ^2\left( {v}_\eta , {T}_\eta , Z_\eta \right) \ge 0, \quad \text {at} \quad \xi = -1, \nonumber \\&\quad {\widehat{T}}_\eta {\widehat{v}}_\eta = -\frac{1-\gamma _\eta ^2}{Z_{\eta }}{p}_\eta ^2\left( {v}_\eta , {T}_\eta , Z_\eta \right) \le 0 , \quad \text {at} \quad \xi = 1. \end{aligned}$$

(141c)

The algebraic identities (141a)–(141c) will be crucial in proving numerical stability. Please see [17] for more details.

1.2 Interface data

To begin, define the outgoing characteristics at the interface

$$\begin{aligned} q_\eta ^{+} = \frac{1}{2}\left( Z_\eta ^{+} v^{+}_\eta + T_\eta ^{+}\right) , \quad p_\eta ^{-} = \frac{1}{2}\left( Z_\eta ^{-} v^{-}_\eta - T_\eta ^{-}\right) , \quad \eta = x, y, z, \end{aligned}$$

(142)

where $Z_\eta ^{\pm } > 0$ are the impedances. We define the hat-variables preserving the amplitude of the outgoing characteristics at the interface

$$\begin{aligned}&{\widehat{q}}_\eta ^{+} \left( {\widehat{v}}_\eta ^{+}, {\widehat{T}}_\eta ^{+}, Z_\eta ^{+}\right) = {q}_\eta ^{+} \left( {v}_\eta ^{+}, {T}_\eta ^{+}, Z_\eta ^{+} \right) , \nonumber \\&{\widehat{p}}_\eta ^{-} \left( {\widehat{v}}_\eta ^{-}, {\widehat{T}}_\eta ^{-}, Z_\eta ^{-}\right) = {p}_\eta ^{-}\left( {v}_\eta ^{-}, {T}_\eta ^{-}, Z_\eta ^{-} \right) . \end{aligned}$$

(143)

The hat-variables also satisfy the interface conditions (24) exactly. For each setup, given ${q}_\eta ^{+}, {p}_\eta ^{-}$, the procedure will solve (143), and (24) for the hat-variables, ${\widehat{v}}_\eta ^{\pm }, {\widehat{T}}_\eta ^{\pm }$. We consider an interface separating two elastic solids. The hat-variables must satisfy (143) and the interface conditions, (24), that is, force balance, no opening and no slip conditions. Combining the two equations in (143) and ensuring force balance, ${\widehat{T}}_\eta ^{-} = {\widehat{T}}_\eta ^{+} = {\widehat{T}}_\eta$, we have

$$\begin{aligned}&{\widehat{T}}_\eta + \alpha _\eta \llbracket {{\widehat{v}}_\eta \rrbracket } = \varPhi _\eta , \quad \varPhi _\eta = \alpha _\eta \left( \frac{2}{Z_\eta ^{+}}q_\eta ^{+} - \frac{2}{Z_\eta ^{-}}p_\eta ^{-}\right) , \nonumber \\&\alpha _\eta = \frac{Z_\eta ^{+}Z_\eta ^{-}}{Z_\eta ^{+} + Z_\eta ^{-}} >0, \quad \eta = x, y, z. \end{aligned}$$

(144)

Furthermore, enforcing no opening and no slip conditions, $\llbracket {{\widehat{v}}_\eta \rrbracket } = 0$, in (144) gives

$$\begin{aligned} {\widehat{T}}_\eta = \varPhi _\eta , \quad \llbracket {{\widehat{v}}_\eta \rrbracket } = 0, \quad \eta = x, y, z. \end{aligned}$$

(145)

We can now define explicitly the hat-variables corresponding to the particle velocities and tractions

$$\begin{aligned} {\widehat{T}}_\eta ^{-} = {\widehat{T}}_\eta ^{+} = {\widehat{T}}_\eta , \quad {\widehat{v}}_\eta ^{+} = \frac{2p_\eta ^{-}+{\widehat{T}}_\eta }{Z_\eta ^{-}} , \quad {\widehat{v}}_\eta ^{-} = \frac{2q_\eta ^{+}-{\widehat{T}}_\eta }{Z_\eta ^{+}} , \quad \quad \eta = x, y, z. \end{aligned}$$

(146)

Note that we have equivalently redefined the physical interface condition (24) as follows

$$\begin{aligned} v_{\eta }^{\pm } = {\widehat{v}}_{\eta }^{\pm }, \quad T_{\eta }^{\pm } = {\widehat{T}}_{\eta }, \quad \eta = x, y, z. \end{aligned}$$

(147)

By construction, the hat-variables ${\widehat{v}}_\eta ^{\pm }$, ${\widehat{T}}_\eta ^{\pm }$ satisfy the following algebraic identities:

$$\begin{aligned}&{q}_\eta \left( {\widehat{v}}_\eta ^{+}, {\widehat{T}}_\eta ^{+}, Z_\eta ^{+}\right) = {q}_\eta \left( {v}_\eta ^{+}, {T}_\eta ^{+}, Z_\eta ^{+} \right) , \quad {p}_\eta \left( {\widehat{v}}_\eta ^{-}, {\widehat{T}}_\eta ^{-}, Z_\eta ^{-}\right) = {p}_\eta \left( {v}_\eta ^{-}, {T}_\eta ^{-}, Z_\eta ^{-} \right) , \end{aligned}$$

(148a)

$$\begin{aligned}&\left( q_\eta ^2 \left( {v}_\eta ^{+}, {T}_\eta ^{+}, Z_\eta ^{+}\right) \right) -{p}_\eta ^2 \left( {\widehat{v}}_\eta ^{+}, {\widehat{T}}_\eta ^{+}, Z_\eta ^{+}\right) = Z_\eta ^{+}{\widehat{T}}_\eta {\widehat{v}}_\eta ^{+}, \nonumber \\&\quad p^2_\eta \left( {v}_\eta ^{-}, {T}_\eta ^{-}, Z_\eta ^{-} \right) -{q}^2_\eta \left( {\widehat{v}}_\eta ^{-}, {\widehat{T}}_\eta ^{-}, Z_\eta ^{-}\right) = -Z_\eta ^{-}{\widehat{T}}_\eta {\widehat{v}}_\eta ^{-}, \end{aligned}$$

(148b)

$$\begin{aligned}&\frac{1}{Z_\eta ^+}\left( q^2_\eta \left( {v}_\eta ^{+}, {T}_\eta ^{+}, Z_\eta ^{+} \right) - {p}^2_\eta \left( {\widehat{v}}_\eta ^{+}, {\widehat{T}}_\eta ^{+}, Z_\eta ^{+}\right) \right) + \frac{1}{Z_\eta ^-}\left( p^2_\eta \left( {v}_\eta ^{-}, {T}_\eta ^{-}, Z_\eta ^{-} \right) \right. \nonumber \\&\quad \left. -{q}^2_\eta \left( {\widehat{v}}_\eta ^{-}, {\widehat{T}}_\eta ^{-}, Z_\eta ^{-}\right) \right) = {\widehat{T}}_\eta \llbracket {{\widehat{v}}_\eta \rrbracket } \equiv 0, \end{aligned}$$

(148c)

The identities (148a)–(148b) can be easily verified, see [17], and will be useful in the prove of numerical stability.