# A mass conservative Eulerian-Lagrangian Runge-Kutta discontinuous Galerkin
method for wave equations with large time stepping.††thanks:  Research of the
first author is supported by the China Scholarship Council for 2 years’ study
at the University of Delaware, CSC NO. 201906340216. Research of the second
author is supported by NSF grant NSF-DMS-2111253, Air Force Office of
Scientific Research FA9550-18-1-0257 and University of Delaware.

Xue Hong School of Mathematical Sciences, University of Science and Technology
of China, Hefei, Anhui, 230026, P.R. China. (xuehong1@mail.ustc.edu.cn). Jing-
Mei Qiu Department of Mathematical Sciences, University of Delaware, Newark,
DE, 19716, USA. (jingqiu@udel.edu).

###### Abstract

We propose an Eulerian-Lagrangian (EL) Runge-Kutta (RK) discontinuous Galerkin
(DG) method for wave equations. The method is designed based on the ELDG
method for transport problems [J. Comput. Phy. 446: 110632, 2021.], which
tracks solution along approximations to characteristics in the DG framework,
allowing extra large time stepping sizes with stability. The wave equation can
be written as a first order hyperbolic system. Considering each characteristic
family, a straightforward application of ELDG will be to transform to the
characteristic variables, evolve them on associated characteristic related
space-time regions, and transform them back to the original variables.
However, the mass conservation could not be guaranteed in a general setting.
In this paper, we formulate a mass conservative semi-discrete ELDG method by
decomposing each variable into two parts, each of them associated with a
different characteristic family. As a result, four different quantities are
evolved in EL fashion and recombined to update the solution. The fully
discrete scheme is formulated by using method-of-lines RK methods, with
intermediate RK solutions updated on the background mesh. Numerical results on
1D and 2D wave equations are presented to demonstrate the performance of the
proposed ELDG method. These include the high order spatial and temporal
accuracy, stability with extra large time stepping size, and mass conservative
property.

Key words: Eulerian-Lagrangian; discontinuous Galerkin; characteristic method;
wave equations; hyperbolic system; mass conservative property.

##  1 Introduction

In this paper, we propose an Eulerian-Lagrangian (EL) discontinuous Galerkin
(DG) method for the wave equation in the form of

$u_{tt}-\nabla\cdot(\mathbf{Q}(\mathbf{x},t)\nabla u)=0,\
(\mathbf{x},t)\in\mathbb{R}^{d}\times[0,T],$ (1.1)

where $d$ is the spatial dimension,
$u:\mathbb{R}^{d}\times[0,T]\rightarrow\mathbb{R}$, and
$\mathbf{Q}(\mathbf{x},t)=diag(q_{1}(\mathbf{x},t),\cdots,q_{d}(\mathbf{x},t))$.
Such a model could come from many fields of engineering, e.g. petroleum
engineering, defense industry, telecommunication and geoscience ([12], [19]).
There are many efficient DG methods for solving the second order wave
equation. These include the interior penalty DG (IPDG) method[32, 14, 15],
local DG method[1, 36, 8]. One can also reformulate the problem as a first-
order hyperbolic system

$U_{t}+\sum_{j=1}^{d}(A_{j}(\mathbf{x},t)U)_{x_{j}}=0,\
(\mathbf{x},t)\in\mathbb{R}^{d}\times[0,T],$ (1.2)

for which DG [10] and space-time DG [13, 29] methods are available.

One effective numerical approach for solving the hyperbolic equation is the
characteristic method, which updates time-dependent solution by tracking
characteristics. For scalar hyperbolic equations, there have been many
pioneering works. In [7], Celia etc. developed an Eulerian Lagrangian
Localized Adjoint Methods (ELLAM) [7], which introduces an adjoint problem for
each test function in the continuous finite element framework and has been
applied to different problems [34, 33]. On the other hand, the ELDG [35, 6]
and SLDG [4] are being developed in the discontinuous Galerkin finite element
framework with a similar introduction of adjoint problems for test functions.
Another line of development, that is closely related to this work, is the
Arbitrary Lagrangian Eulerian (ALE) DG method [20, 17]. Both ELDG and ALE-DG
evolve the DG solution on a dynamic moving mesh: the dynamic mesh movement of
the ELDG approximates characteristics for the potential of using larger time
stepping sizes with stability, whereas the mesh movement of ALE-DG is designed
for better shock resolution.

In this paper, we propose a mass conservative EL RK DG method for hyperbolic
systems (1.2) with large time stepping sizes. We start from 1D cases, for
which we consider characteristic variables and the associated characteristic
space-time regions. For each characteristic family, a straightforward
application of ELDG will be to transform to characteristic variables, evolve
them on associated space-time regions, and transform them back to original
variables. However, the mass conservation could not be guaranteed in a general
setting. We decompose each variable into two parts, each of them associated
with different characteristic families; as a result, four different quantities
are evolved in EL fashion and recombined to update the solution. The fully
discrete scheme is formulated by using method-of-lines RK methods, with
intermediate RK solutions updated on the background mesh. For 2D hyperbolic
systems, characteristic variables are no longer constant along the
bicharacteristics, characteristic Galerkin [30] or evolution Galerkin [3, 22,
23] methods have been proposed by taking into account information propagated
in all bicharacteristic directions. However, the algorithm implementation is
very complex. In this paper, we use the splitting method for higher
dimensional problem, which maintain the simplicity of ELDG method for 1D
cases, as well as other great properties such as mass conservation, high order
spatial and temporal accuracy, and allows for extra large time steps with
stability.

This paper is organized as follows. In Section 2, we review the ELDG for one-
dimensional (1D) linear transport problems. In Section 3, we develop the ELDG
method for a first-order hyperbolic system by evolving each component
associated with each characteristic families, recombining them to update the
solution. We also develop it to 2D problems by dimensional splitting. Mass
conservation of ELDG schemes are proved. In Section 4, performance of the
proposed ELDG method is shown through extensive numerical tests. Finally,
concluding remarks are made in Section 5.

##  2 Review of ELDG formulation for 1D linear transport problems

To illustrate the key idea of the ELDG scheme, we start from a 1D linear
transport equation in the following form

$u_{t}+(a(x,t)u)_{x}=0,\quad x\in[x_{a},x_{b}].$ (2.1)

For simplicity, we assume periodic boundary conditions, and the velocity field
$a(x,t)$ is a continuous function of space and time. We perform a partition of
the computational domain
$x_{a}=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{N+\frac{1}{2}}=x_{b}$ as the
background mesh partition. Let $I_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$
denote an element of length $\Delta x_{j}=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$
and define $\Delta x=\max_{j}\Delta x_{j}.$ We define the finite dimensional
approximation space, $V_{h}^{k}=\\{v_{h}:v_{h}|_{I_{j}}\in P^{k}(I_{j})\\}$,
where $P^{k}(I_{j})$ denotes the set of polynomials of degree at most $k$. We
let $t^{n}$ be the $n$-th time level and $\Delta t=t^{n+1}-t^{n}$ to be the
time-stepping size.

The key idea in the ELDG formulation is, design adjoint problems for test
functions to take advantage of information propagation along characteristics.
The ELDG method proposed in [6] is formulated on a space-time region
$\Omega_{j}=\tilde{I}_{j}(t)\times[t^{n},t^{n+1}]$ with

$\tilde{I}_{j}(t)=[\tilde{x}_{j-\frac{1}{2}}(t),\tilde{x}_{j+\frac{1}{2}}(t)],\quad
t\in[t^{n},t^{n+1}]$

being the dynamic interval, see Figure 2.1. Here
$\tilde{x}_{j\pm\frac{1}{2}}(t)=x_{j\pm\frac{1}{2}}+(t-t^{n+1})\nu_{j\pm\frac{1}{2}}$
are straight lines emanating from cell boundaries $x_{j\pm\frac{1}{2}}$ with
slopes $\nu_{j\pm\frac{1}{2}}=a(x_{j\pm\frac{1}{2}},t^{n+1})$ and
${I}^{\star}_{j}\doteq\tilde{I}_{j}(t^{n})=[x^{*}_{j-\frac{1}{2}},x^{*}_{j+\frac{1}{2}}]$
is the upstream cell of $I_{j}$ at $t^{n}$. The dynamic interval of
$\tilde{I}_{j}(t)$ can always be linearly mapped to a reference cell $\xi$ in
$I_{j}$ by a mapping.

$\Omega_{j}$
$t^{n}$$t^{n+1}$$x_{j-\frac{1}{2}}$$x_{j+\frac{1}{2}}$$x_{j-\frac{1}{2}}^{\star}$
$x_{j+\frac{1}{2}}^{\star}$$\widetilde{I}_{j}(t^{n})=I_{j}^{\star}$$\widetilde{I}_{j}(t^{n+1})=I_{j}$$\widetilde{I}_{j}(t)$$\nu_{j-\frac{1}{2}}$
$\nu_{j+\frac{1}{2}}$ $\xi=x_{j-1/2}$$\xi=x_{j+1/2}$ Figure 2.1: Illustration
for dynamic element $\widetilde{I}_{j}(t)$ of ELDG.

A local adjoint problem of (2.1) for the test function is defined as:

$\begin{cases}\psi_{t}+\alpha(x,t)\psi_{x}=0,\ (x,t)\in\Omega_{j},\\\
\psi(t=t^{n+1})=\Psi(x).\quad\forall\Psi(x)\in P^{k}(I_{j}).\end{cases}$ (2.2)

Here $\alpha(x,t)$ is a bilinear function of $(x,t)$ with $\forall
t\in[t^{n},t^{n+1}]$,

$\alpha(x_{j-\frac{1}{2}},t^{n+1})=a(x_{j-\frac{1}{2}},t^{n+1})\doteq\nu_{j-\frac{1}{2}},\quad{\alpha}(x_{j+\frac{1}{2}},t^{n+1})=a(x_{j+\frac{1}{2}},t^{n+1})\doteq\nu_{j+\frac{1}{2}},$
(2.3)

and

$\alpha(x,t)=-\nu_{j-\frac{1}{2}}\frac{x-\tilde{x}_{j+\frac{1}{2}}(t)}{\Delta
x_{j}(t)}+\nu_{j+\frac{1}{2}}\frac{x-\tilde{x}_{j-\frac{1}{2}}(t)}{\Delta
x_{j}(t)}\in
P^{1}([\tilde{x}_{j-\frac{1}{2}}(t),\tilde{x}_{j+\frac{1}{2}}(t)]).$ (2.4)

Notice that the test function $\psi$ stays the same polynomial, if
$\tilde{I}_{j}(t)$ is mapped to a reference interval $I_{j}$. The ELDG [6]
scheme can be formulated by $\int_{\tilde{I}_{j}(t)}\eqref{ELDGadjoint1d}\cdot
u+\eqref{scalar1d}\cdot\psi$

$\frac{d}{dt}\int_{\tilde{I}_{j}(t)}(u\psi)dx=-\left(\hat{F}\psi\right)\left|{}_{\tilde{x}_{j+\frac{1}{2}}(t)}\right.+\left(\hat{F}\psi\right)\left|{}_{\tilde{x}_{j-\frac{1}{2}}(t)}\right.+\int_{\tilde{I}_{j}(t)}F\psi_{x}dx.$
(2.5)

where $F(u)\doteq(a-\alpha)u$ and $\hat{F}$ is the Lax-Friedrichs flux. A
method-of-lines RK discretization can be used for high order temporal accuracy
[6].

##  3 The ELDG algorithm for wave equation.

The design of the ELDG algorithm for wave equation shares a similar spirit as
the 1D scalar case. We start from a 1D wave equation. We first rewrite it as a
linear system and formulate the ELDG scheme by tracking information along
different characteristics families.

###  3.1 1D wave equation

We consider the wave equation

$u_{tt}=(a^{2}(x)u_{x})_{x}+f(x,t).$ (3.1)

For simplicity, we assume periodic boundary conditions, and the velocity field
$a(x)$ is a continuous and periodic function of space. Let $u^{1}=u_{t}$ and
$u^{2}=u_{x}$, then we can rewrite (3.1) as a linear system

$U_{t}+(AU)_{x}=F,$ (3.2)

where

$U=\begin{bmatrix}u^{1}\\\ u^{2}\end{bmatrix},\qquad
A(x)=\begin{bmatrix}0&-a^{2}(x)\\\ -1&0\end{bmatrix},\qquad
F=\begin{bmatrix}f(x,t)\\\ 0\end{bmatrix}.$

We use the following notations for the eigen-decomposition of $A(x)$:

  * •

eigenvalue: $\lambda^{(1)}(x),\quad\lambda^{(2)}(x),$ where
$\lambda^{(1)}(x)=a(x),\quad\lambda^{(2)}(x)=-a(x)$ for wave equation.

  * •

$A(x)=R(x)\Lambda R^{-1}(x)$, where
$\Lambda(x)=diag(\lambda^{(1)}(x),\lambda^{(2)}(x))$,

$R(x)=\begin{bmatrix}r_{11}(x)&r_{12}(x)\\\
r_{21}(x)&r_{22}(x)\end{bmatrix}=\begin{bmatrix}-a(x)&a(x)\\\
1&1\end{bmatrix}$ (3.3)

contains the right eigenvectors, and

$R^{-1}(x)=\begin{bmatrix}l_{11}(x)&l_{12}(x)\\\
l_{21}(x)&l_{22}(x)\end{bmatrix}=\begin{bmatrix}\frac{-1}{2a(x)}&\frac{1}{2}\\\
\frac{1}{2a(x)}&\frac{1}{2}\end{bmatrix}\doteq\begin{bmatrix}l^{(1)}(x)\\\
l^{(2)}(x)\end{bmatrix}$ (3.4)

contains the left eigenvectors.

In the following, we propose a mass conservative ELDG scheme for the system
(3.2) by the procedure below. The notion of its background meshes is the same
as the 1D scalar case.

(1) Two partitions of space-time regions $\Omega^{(1)}_{j}$ and
$\Omega^{(2)}_{j}$. According to the first and second characteristic families,
we partition the computational domain as two sets of space-time regions
${\Omega^{(1)}_{j}}$ and ${\Omega^{(2)}_{j}}$ respectively. Here
$\Omega^{(1)}_{j}=\tilde{I}^{(1)}_{j}(t)\times[t^{n},t^{n+1}]$ is related to
the first characteristic family.
$\tilde{I}^{(1)}_{j}(t)=[\tilde{x}^{(1)}_{j-\frac{1}{2}}(t),\tilde{x}^{(1)}_{j+\frac{1}{2}}(t)]$
is the dynamic interval emanating from cell boundaries $x_{j\pm\frac{1}{2}}$
with slopes ${\nu}^{(1)}_{j\pm\frac{1}{2}}$ approximating the first
characteristic velocity, see Figure 3.2 (left). In general, we choose
$\nu^{(1)}_{j\pm\frac{1}{2}}=\lambda^{(1)}(x_{j\pm\frac{1}{2}})$.
$I^{\star,(1)}_{j}\doteq\tilde{I}^{(1)}_{j}(t^{n})$ is the upstream cell of
$I_{j}$ from the first characteristic family at $t^{n}$. The dynamic interval
$\tilde{I}^{(1)}_{j}(t)$ can be linearly mapped to a reference cell $\xi\in
I_{j}$, (see Figure 2.1). Here, we let $\tilde{x}^{(1)}(\tau;(\xi,t^{n+1}))$
be the linear map from $\tilde{I}^{(1)}_{j}(t)$ to $I_{j}$. Similar definition
can be made to $\Omega^{(2)}_{j}$, $\tilde{I}^{(2)}_{j}(t)$ and
$I^{\star,(2)}_{j}$ for the second characteristic family. See Figure 3.2
(right) for illustration of $\Omega^{(2)}_{j}$.

$\Omega^{(1)}_{j}$
$t^{n}$$t^{n+1}$$x_{j-\frac{1}{2}}$$x_{j+\frac{1}{2}}$$x^{\star,(1)}_{j-\frac{1}{2}}$
$x^{\star,(1)}_{j+\frac{1}{2}}$$\tilde{I}^{(1)}_{j}(t^{n})=I^{\star,(1)}_{j}$$\tilde{I}^{(1)}_{j}(t^{n+1})=I_{j}$$\tilde{I}^{(1)}_{j}(t)$$\nu^{(1)}_{j-\frac{1}{2}}$
$\nu^{(1)}_{j+\frac{1}{2}}$ $\Omega^{(2)}_{j}$
$t^{n}$$t^{n+1}$$x_{j-\frac{1}{2}}$$x_{j+\frac{1}{2}}$$x^{\star,(2)}_{j-\frac{1}{2}}$
$x^{\star,(2)}_{j+\frac{1}{2}}$$\tilde{I}^{(2)}_{j}(t^{n})=I^{\star,(2)}_{j}$$\tilde{I}^{(2)}_{j}(t^{n+1})=I_{j}$$\tilde{I}^{(2)}_{j}(t)$$\nu^{(2)}_{j-\frac{1}{2}}$
$\nu^{(2)}_{j+\frac{1}{2}}$ Figure 3.2: Illustration for dynamic elements
$\widetilde{I}^{(1)}_{j}(t)$ (left) and $\widetilde{I}^{(2)}_{j}(t)$ (right)
of ELDG for the first and second characteristic families of the system.

(2) Adjoint Problems. We consider an adjoint problem for the first
characteristic family on $\Omega^{(1)}_{j}$:

$\begin{cases}(\psi^{(1)})_{t}+\alpha^{(1)}(\psi^{(1)})_{x}=0,\
t\in[t^{n},t^{n+1}],\\\ (\psi^{(1)})(t=t^{n+1})=\Psi^{(1)}(x),\end{cases}$
(3.5)

where

$\alpha^{(1)}(x,t)=-\nu^{(1)}_{j-\frac{1}{2}}\frac{x-\tilde{x}^{(1)}_{j+\frac{1}{2}}(t)}{\Delta
x_{j}^{(1)}(t)}+\nu^{(1)}_{j+\frac{1}{2}}\frac{x-\tilde{x}^{(1)}_{j-\frac{1}{2}}(t)}{\Delta
x_{j}^{(1)}(t)}\in P^{1}(\tilde{I}^{(1)}_{j}(t)).$ (3.6)

Similarly on $\Omega^{(2)}_{j}$:

$\begin{cases}(\psi^{(2)})_{t}+\alpha^{(2)}(\psi^{(2)})_{x}=0,\
t\in[t^{n},t^{n+1}],\\\ (\psi^{(2)})(t=t^{n+1})=\Psi^{(2)}(x),\end{cases}$
(3.7)

where

$\alpha^{(2)}(x,t)=-\nu^{(2)}_{j-\frac{1}{2}}\frac{x-\tilde{x}^{(2)}_{j+\frac{1}{2}}(t)}{\Delta
x_{j}^{(2)}(t)}+\nu^{(2)}_{j+\frac{1}{2}}\frac{x-\tilde{x}^{(2)}_{j-\frac{1}{2}}(t)}{\Delta
x_{j}^{(2)}(t)}\in P^{1}(\tilde{I}^{(2)}_{j}(t)).$ (3.8)

(3) Formulation of a mass-conservative semi-discrete ELDG scheme. Before we
present a mass-conservative semi-discrete ELDG scheme, we formulate the scheme
by a localized characteristic field. In particular, a piecewise constant
$a_{j}$ approximating $a(x)$ in (3.4) is defined on $I_{j}$, and the
corresponding

$R_{j}=\begin{bmatrix}r_{j}^{11}&r_{j}^{12}\\\
r_{j}^{21}&r_{j}^{22}\end{bmatrix}=\begin{bmatrix}-a_{j}&a_{j}\\\
1&1\end{bmatrix}$ (3.9)

and

$R_{j}^{-1}=\begin{bmatrix}\frac{-1}{2a_{j}}&\frac{1}{2}\\\
\frac{1}{2a_{j}}&\frac{1}{2}\end{bmatrix}\doteq\begin{bmatrix}l_{j}^{(1)}\\\
l_{j}^{(2)}\end{bmatrix}.$ (3.10)

Define $l_{j}^{(1)}$ is locally defined on $\Omega_{j}^{(1)}$ approximating
$l^{(1)}(x)$. For simplicity, we only present the first order ELDG scheme.
Take the vector product of $l^{(1)}_{j}$ from left with (3.2), we have a
scalar equation

$l^{(1)}_{j}(U_{t}+(A(x)U)_{x})=l^{(1)}_{j}F(x,t).$ (3.11)

Next, integrating over the space-time interval $\tilde{I}^{(1)}_{j}(t)$, then
we have

$\displaystyle\frac{d}{dt}\int_{\tilde{I}^{(1)}_{j}(t)}(l^{(1)}_{j}U)dx+l^{(1)}_{j}(A(x)U-\nu^{(1)}_{j+\frac{1}{2}}U)|_{\tilde{x}^{(1)}_{j+\frac{1}{2}}(t)}-l^{(1)}_{j}(A(x)U-\nu^{(1)}_{j-\frac{1}{2}}U)|_{\tilde{x}^{(1)}_{j-\frac{1}{2}}(t)}$
(3.12) $\displaystyle=\int_{\tilde{I}^{(1)}_{j}(t)}l^{(1)}_{j}F(x,t)dx.$

The first order ELDG discretization of eq. (3.12) is to find
$l^{(1)}_{j}U_{h}(x,t)\in P^{0}(\tilde{I}^{(1)}_{j}(t))$, so that

$\displaystyle\frac{d}{dt}\int_{\tilde{I}^{(1)}_{j}(t)}l^{(1)}_{j}U_{h}dx=-\left[l^{(1)}_{j}\widehat{(A(x)U_{h}-\nu^{(1)}_{j+\frac{1}{2}}U_{h})}|_{j+\frac{1}{2}}\right]+\left[l^{(1)}_{j}\widehat{(A(x)U_{h}-\nu^{(1)}_{j-\frac{1}{2}}U_{h})}|_{j-\frac{1}{2}}\right]$
(3.13) $\displaystyle+\int_{\tilde{I}^{(1)}_{j}(t)}l^{(1)}_{j}F(x,t)dx\doteq
L_{1}(U_{h}(t),t,\tilde{I}^{(1)}_{j}(t)).$

Here $\widehat{f^{11}}$ at a cell boundary can be taken as a monotone flux,
e.g. the Lax-Friedrichs flux

$\widehat{(AU-\nu^{(1)}_{j+\frac{1}{2}}U)}_{{j+\frac{1}{2}}}=\frac{1}{2}\left(A(x^{+}_{j+\frac{1}{2}})U^{+}_{j+\frac{1}{2}}-\nu^{(1)}_{j+\frac{1}{2}}U_{j+\frac{1}{2}}^{+}+A(x^{-}_{j+\frac{1}{2}})U_{j+\frac{1}{2}}^{-}-\nu^{(1)}_{j+\frac{1}{2}}U_{j+\frac{1}{2}}^{-}-\alpha_{1,2}(U_{j+\frac{1}{2}}^{+}-U_{j+\frac{1}{2}}^{-})\right)$

where
$\alpha_{1,2}=max\\{|\lambda^{(1)}(x_{j+\frac{1}{2}})-\nu^{(1)}_{j+\frac{1}{2}}|,|\lambda^{(2)}(x_{j+\frac{1}{2}})-\nu^{(2)}_{j+\frac{1}{2}}|\\}$.
Similarly, we can easily update $l^{(2)}_{j}U_{h}$ related to $\lambda^{(2)}$
in the following:

$\displaystyle\frac{d}{dt}\int_{\tilde{I}^{(2)}_{j}(t)}l^{(2)}_{j}U_{h}dx=-\left[l^{(2)}_{j}\widehat{(A(x)U_{h}-\nu^{(2)}_{j+\frac{1}{2}}U_{h})}|_{j+\frac{1}{2}}\right]+\left[l^{(2)}_{j}\widehat{(A(x)U_{h}-\nu^{(2)}_{j-\frac{1}{2}}U_{h})}|_{j-\frac{1}{2}}\right]$
(3.14) $\displaystyle+\int_{\tilde{I}^{(2)}_{j}(t)}l^{(2)}_{j}F(x,t)dx\doteq
L_{2}(U_{h}(t),t,\tilde{I}^{(2)}_{j}(t)),$

where $l_{j}^{(2)}$ is a constant vector, locally defined on
$\tilde{I}^{(2)}_{j}(t)$ approximating $l^{(2)}(x)$. A simple first order ELDG
scheme is composed by two evolution equations (3.13) and (3.14). That is, we
can update $u_{h}^{1}$ by (3.13),(3.14) and
$u_{h}^{1}=(r_{j}^{11}l^{(1)}_{j}+r_{j}^{12}l^{(2)}_{j})U_{h}$:

$\displaystyle\int_{I_{j}}u^{1,n+1}_{h}dx$
$\displaystyle=\int_{I_{j}}r_{j}^{11}l^{(1)}_{j}U_{h}^{n+1}dx+\int_{I_{j}}r_{j}^{12}l^{(2)}_{j}U_{h}^{n+1}dx$
(3.15)
$\displaystyle=r_{j}^{11}\int_{I_{j}^{*,(1)}}l^{(1)}_{j}U_{h}^{n}dx+r_{j}^{11}\int_{t^{n}}^{t^{n+1}}L_{1}(U^{(1)}_{h}(t),t,\tilde{I}^{(1)}_{j}(t))dt$
$\displaystyle+r_{j}^{12}\int_{I_{j}^{*,(2)}}l^{(2)}_{j}U_{h}^{n}dx+r_{j}^{12}\int_{t^{n}}^{t^{n+1}}L_{2}(U^{(2)}_{h}(t),t,\tilde{I}^{(2)}_{j}(t))dt,$

where $U_{h}^{n}$ and $U_{h}^{n+1}$ are defined on the background mesh
$I_{j}$, $U^{(1)}_{h}(t)$ and $U^{(2)}_{h}(t)$ are defined on the space-time
dynamic meshes $\tilde{I}^{(1)}_{j}(t)$ and $\tilde{I}^{(2)}_{j}(t)$
respectively. Similarly, we can update $u_{h}^{2}$.

We apply forward-Euler method for time discretization with above ELDG scheme
(3.15):

$\displaystyle\int_{I_{j}}u^{1,n+1}dx$ (3.16)
$\displaystyle=r_{j}^{11}\int_{I_{j}^{*,(1)}}l^{(1)}_{j}U^{n}dx-\Delta t\
r_{j}^{11}l^{(1)}_{j}\left[\widehat{(A(x)U^{n}-\nu^{(1)}_{j+\frac{1}{2}}U^{n})}|_{x^{\star,(1)}_{j+\frac{1}{2}}}-\widehat{(A(x)U-\nu^{(1)}_{j-\frac{1}{2}}U)}|_{x^{\star,(1)}_{j-\frac{1}{2}}}\right]$
$\displaystyle+r_{j}^{12}\int_{I_{j}^{*,(2)}}l^{(2)}_{j}U^{n}dx-\Delta t\
r_{j}^{12}l^{(2)}_{j}\left[\widehat{(A(x)U^{n}-\nu^{(2)}_{j+\frac{1}{2}}U^{n})}|_{x^{\star,(2)}_{j+\frac{1}{2}}}-\widehat{(A(x)U-\nu^{(2)}_{j-\frac{1}{2}}U)}|_{x^{\star,(2)}_{j-\frac{1}{2}}}\right].$

######  Remark 3.1.

The above ELDG scheme is not mass conservative for two reasons:

  * (1)

Flux terms can’t cancel each other as $r_{j}^{11}l^{(1)}_{j}$ and
$r_{j}^{12}l^{(2)}_{j}$ are discontinuous across cell boundary of
$\Omega_{j}^{(1)},\Omega_{j}^{(2)}$.

  * (2)

$\sum_{j}r_{j}^{11}\int_{I_{j}^{*,(1)}}l^{(1)}_{j}U^{n}dx+r_{j}^{12}\int_{I_{j}^{*,(2)}}l^{(2)}_{j}U^{n}dx\neq\sum_{j}\int_{I_{j}}U^{n}dx$
because of the inconsistency in characteristic transformations between
neighboring cells among $r_{j}^{11}l^{(1)}_{j}$ and $r_{j}^{12}l^{(2)}_{j}$.

To obtain the mass conservation, a critical point is that eigen-decomposition
has to be consistent among two characteristic families and independent of
partitions $\Omega_{j}^{(1)},\Omega_{j}^{(2)}$. We choose $a_{p}(x)$ as a
continuous approximation of $a(x)$ which leads that
$R_{p}(x),\Lambda_{p}(x),R_{p}^{-1}(x)$ are the approximations of
$R(x),\Lambda(x),R^{-1}(x)$ respectively, where

$R_{p}(x)=\begin{bmatrix}r^{p}_{11}(x)&r^{p}_{12}(x)\\\
r^{p}_{21}(x)&r^{p}_{22}(x)\end{bmatrix}=\begin{bmatrix}-a_{p}(x)&a_{p}(x)\\\
1&1\end{bmatrix}$ (3.17)

and

$R_{p}^{-1}(x)=\begin{bmatrix}l^{p}_{11}(x)&l^{p}_{12}(x)\\\
l^{p}_{21}(x)&l^{p}_{22}(x)\end{bmatrix}=\begin{bmatrix}\frac{-1}{2a_{p}(x)}&\frac{1}{2}\\\
\frac{1}{2a_{p}(x)}&\frac{1}{2}\end{bmatrix}\doteq\begin{bmatrix}l_{p}^{(1)}(x)\\\
l_{p}^{(2)}(x)\end{bmatrix}.$ (3.18)

The following equalities hold

$\displaystyle(r^{p}_{11}(x)l_{p}^{(1)}(x)+r^{p}_{12}(x)l_{p}^{(2)}(x))U(x)$
$\displaystyle=u^{1},$ (3.19)
$\displaystyle(r^{p}_{21}(x)l_{p}^{(1)}(x)+r^{p}_{22}(x)l_{p}^{(2)}(x))U(x)$
$\displaystyle=u^{2},$

from $R_{p}(x)R_{p}^{-1}(x)=R_{p}^{-1}(x)R_{p}(x)=I$. They are critical to
design a mass-conservative ELDG scheme. Take the vector product of
$r^{p}_{11}(x)l_{p}^{(1)}(x)$ from left with (3.2), we have a scalar equation

$r^{p}_{11}(x)l_{p}^{(1)}(x)(U_{t}+(A(x)U)_{x})=r^{p}_{11}(x)l_{p}^{(1)}(x)F(x,t).$
(3.20)

Multiply $\psi^{(1)}$ to the equation above,

$r^{p}_{11}(x)l_{p}^{(1)}(x)(U_{t}+(A(x)U)_{x})\psi^{(1)}=r^{p}_{11}(x)l_{p}^{(1)}(x)F(x,t)\psi^{(1)}.$
(3.21)

Meanwhile, from $r^{p}_{11}(x)l_{p}^{(1)}(x)U$ $\cdot$ (3.5), we have

$r^{p}_{11}(x)l_{p}^{(1)}(x)U(\psi^{(1)})_{t}+r^{p}_{11}(x)l_{p}^{(1)}(x)U\alpha^{(1)}(\psi^{(1)})_{x}=0.$
(3.22)

Next, sum equations (3.21) and (3.22), and integrate over the space-time
interval $\Omega^{(1)}_{j}$, then

$\displaystyle\int_{\Omega^{(1)}_{j}}\left(r^{p}_{11}(x)l_{p}^{(1)}(x)U_{t}\psi^{(1)}+r^{p}_{11}(x)l_{p}^{(1)}(x)U(\psi^{(1)})_{t}+r^{p}_{11}(x)l_{p}^{(1)}(x)(A(x)U)_{x}\psi^{(1)}\right)dxdt$
(3.23)
$\displaystyle+\int_{\Omega^{(1)}_{j}}r^{p}_{11}(x)l_{p}^{(1)}(x)U\alpha^{(1)}(\psi^{(1)})_{x}dxdt=\int_{\Omega^{(1)}_{j}}r^{p}_{11}(x)l_{p}^{(1)}(x)F(x,t)\psi^{(1)}dxdt.$

A further manipulation on the L.H.S. of (LABEL:eq:_system_f1) gives

$\displaystyle\int_{\Omega^{(1)}_{j}}\left((r^{p}_{11}(x)l_{p}^{(1)}(x)U\psi^{(1)})_{t}+r^{p}_{11}(x)l_{p}^{(1)}(x)(A(x)U)_{x}\psi^{(1)}+r^{p}_{11}(x)l_{p}^{(1)}(x)U\alpha^{(1)}(\psi^{(1)})_{x}\right)dxdt$
(3.24) $\displaystyle=$
$\displaystyle\int_{\Omega^{(1)}_{j}}\left((r^{p}_{11}(x)l_{p}^{(1)}(x)U\psi^{(1)})_{t}+(r^{p}_{11}(x)l_{p}^{(1)}(x)A(x)U\psi^{(1)})_{x}-(r^{p}_{11}(x)l_{p}^{(1)}(x))_{x}A(x)U\psi^{(1)}\right)dxdt$
$\displaystyle-$
$\displaystyle\int_{\Omega^{(1)}_{j}}\left(r^{p}_{11}(x)l_{p}^{(1)}(x)A(x)U(\psi^{(1)})_{x}-r^{p}_{11}(x)l_{p}^{(1)}(x)U\alpha^{(1)}(\psi^{(1)})_{x}\right)dxdt$
$\displaystyle=$
$\displaystyle\int^{t^{n+1}}_{t^{n}}\left(\frac{d}{dt}\int_{\tilde{I}^{(1)}_{j}(t)}(r^{p}_{11}(x)l_{p}^{(1)}U\psi^{(1)})dx+[r^{p}_{11}(x)l_{p}^{(1)}A(x)U\psi^{(1)}-\nu^{(1)}r^{p}_{11}(x)l_{p}^{(1)}U\psi^{(1)}]|_{j-\frac{1}{2}}^{j+\frac{1}{2}}\right)dt$
$\displaystyle-\int^{t^{n+1}}_{t^{n}}\int_{\tilde{I}^{(1)}_{j}(t)}\left((r^{p}_{11}(x)l_{p}^{(1)}(x))_{x}A(x)U\psi^{(1)}+r^{p}_{11}(x)l_{p}^{(1)}(x)(A(x)U-\alpha^{(1)}U)\psi^{(1)}_{x}\right)dxdt.$

Letting $f^{11}(U)=r^{p}_{11}l_{p}^{(1)}(AU-\alpha^{(1)}U)$, the time
differential form of (LABEL:eq:_system_f1) with (LABEL:eq:leftsystem_f1) gives

$\displaystyle\frac{d}{dt}\int_{\tilde{I}^{(1)}_{j}(t)}(r^{p}_{11}(x)l_{p}^{(1)}(x)U\psi^{(1)})dx+\left(f^{11}\psi^{(1)}\right)\left|{}_{\tilde{x}^{(1)}_{j+\frac{1}{2}}(t)}\right.-\left(f^{11}\psi^{(1)}\right)\left|{}_{\tilde{x}^{(1)}_{j-\frac{1}{2}}(t)}\right.-\int_{\tilde{I}^{(1)}_{j}(t)}f^{11}\psi^{(1)}_{x}dx$
(3.25)
$\displaystyle-\int_{\tilde{I}^{(1)}_{j}(t)}(r^{p}_{11}(x)l_{p}^{(1)}(x))_{x}A(x)U\psi^{(1)}dx=\int_{\tilde{I}^{(1)}_{j}(t)}r^{p}_{11}(x)l_{p}^{(1)}(x)F(x,t)\psi^{(1)}dx.$

Similarly, we have an equation related to $\lambda^{(2)}$

$\displaystyle\frac{d}{dt}\int_{\tilde{I}^{(1)}_{j}(t)}(r^{p}_{12}(x)l_{p}^{(2)}(x)U\psi^{(2)})dx+\left(f^{12}\psi^{(2)}\right)\left|{}_{\tilde{x}^{(2)}_{j+\frac{1}{2}}(t)}\right.-\left(f^{12}\psi^{(2)}\right)\left|{}_{\tilde{x}^{(2)}_{j-\frac{1}{2}}(t)}\right.-\int_{\tilde{I}^{(2)}_{j}(t)}f^{12}\psi^{(2)}_{x}dx$
(3.26)
$\displaystyle-\int_{\tilde{I}^{(2)}_{j}(t)}(r^{p}_{12}(x)l_{p}^{(2)}(x))_{x}A(x)U\psi^{(2)}dx=\int_{\tilde{I}^{(2)}_{j}(t)}r^{p}_{12}(x)l_{p}^{(2)}(x)F(x,t)\psi^{(2)}dx,$

where $f^{12}(U)=r^{p}_{12}l_{p}^{(2)}(AU-\alpha^{(2)}U)$. Then, we can update
$u^{1}$ by (3.19),(3.25), (3.26), taking $\Psi^{(1)}(x)=\Psi(x)$ in (3.5) and
$\Psi^{(2)}(x)=\Psi(x)$ in (3.7):

$\displaystyle\int_{I_{j}}u^{1,n+1}\Psi(x)dx\overset{\eqref{eq:
pickupu1u2}}{=}\int_{I_{j}}r^{p}_{11}l_{p}^{(1)}U^{n+1}\Psi(x)dx+\int_{I_{j}}r^{p}_{12}l_{p}^{(2)}U^{n+1}\Psi(x)dx$
(3.27)
$\displaystyle=\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}\psi^{(1)}dx-\int_{t^{n}}^{t^{n+1}}\left(f^{11}\psi^{(1)}\right)\left|{}_{\tilde{x}^{(1)}_{j+\frac{1}{2}}(t)}\right.+\left(f^{11}\psi^{(1)}\right)\left|{}_{\tilde{x}^{(1)}_{j-\frac{1}{2}}(t)}\right.dt$
$\displaystyle+\int_{t^{n}}^{t^{n+1}}\int_{\tilde{I}^{(1)}_{j}(t)}r^{p}_{11}(x)l_{p}^{(1)}(x)F(x,t)\psi^{(1)}+(r^{p}_{11}(x)l_{p}^{(1)}(x))_{x}A(x)U\psi^{(1)}+f^{11}\psi^{(1)}_{x}dxdt$
$\displaystyle+\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}\psi^{(2)}dx-\int_{t^{n}}^{t^{n+1}}\left(f^{12}\psi^{(2)}\right)\left|{}_{\tilde{x}^{(2)}_{j+\frac{1}{2}}(t)}\right.+\left(f^{12}\psi^{(2)}\right)\left|{}_{\tilde{x}^{(2)}_{j-\frac{1}{2}}(t)}\right.dt$
$\displaystyle+\int_{t^{n}}^{t^{n+1}}\int_{\tilde{I}^{(2)}_{j}(t)}r^{p}_{12}(x)l_{p}^{(2)}(x)F(x,t)\psi^{(2)}+(r^{p}_{12}(x)l_{p}^{(2)}(x))_{x}A(x)U\psi^{(2)}+f^{12}\psi^{(2)}_{x}dxdt.$

As the ELDG method for scalar, the dynamic interval of $I^{(1)}_{j}(t)$ and
$I^{(2)}_{j}(t)$ can always be linearly mapped to a reference cell $\xi\in
I_{j}$, the ELDG discretization of eq. (3.27) is to find $u_{h}^{1}(x,t)\in
P^{k}(I_{j}(t))$, so that

$\displaystyle\int_{I_{j}}u^{1}_{h}(x,t^{n+1})\Psi(x)dx$
$\displaystyle=\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}\psi^{(1)}dx+\int_{t^{n}}^{t^{n+1}}L_{11}(U_{h}(t),t,\tilde{I}^{(1)}_{j}(t))dt$
(3.28)
$\displaystyle+\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}\psi^{(2)}dx+\int_{t^{n}}^{t^{n+1}}L_{12}(U_{h}(t),t,\tilde{I}^{(2)}_{j}(t))dt,$

for $\psi^{(1)}(x,t)$ satisfying the adjoint problem (3.5) with
$\forall\Psi(x)=\psi(x,t^{n+1})\in P^{k}(I_{j})$. Here

$\displaystyle L_{11}(U_{h}(t),t,\tilde{I}^{(1)}_{j}(t))$
$\displaystyle=-\widehat{f^{11}_{j+\frac{1}{2}}}\psi^{(1),-}_{j+\frac{1}{2}}+\widehat{f^{11}_{j-\frac{1}{2}}}\psi^{(1),+}_{j-\frac{1}{2}}+\int_{\tilde{I}^{(1)}_{j}(t)}f^{11}\psi^{(1)}_{x}(x,t)dx$
(3.29)
$\displaystyle+\int_{\tilde{I}^{(1)}_{j}(t)}(r^{p}_{11}l_{p}^{(1)})_{x}AU_{h}\psi^{(1)}(x,t)+r^{p}_{11}l_{p}^{(1)}F\psi^{(1)}(x,t)dx,$
$\displaystyle L_{12}(U_{h}(t),t,\tilde{I}^{(2)}_{j}(t))$
$\displaystyle=-\widehat{f^{12}_{j+\frac{1}{2}}}\psi^{(2),-}_{j+\frac{1}{2}}+\widehat{f^{12}_{j-\frac{1}{2}}}\psi^{(2),+}_{j-\frac{1}{2}}+\int_{\tilde{I}^{(2)}_{j}(t)}f^{12}\psi^{(2)}_{x}(x,t)dx$
$\displaystyle+\int_{\tilde{I}^{(2)}_{j}(t)}(r^{p}_{12}l_{p}^{(2)})_{x}AU_{h}\psi^{(2)}(x,t)+r^{p}_{12}l_{p}^{(2)}F\psi^{(2)}(x,t)dx,$

where
$\psi^{(1),\pm}_{j\pm\frac{1}{2}}=\psi^{(1)}(x^{(1),\pm}_{j\pm\frac{1}{2}}(t),t)$
and $\widehat{f^{11}}$ at a cell boundary can be taken as a monotone flux,
e.g. the Lax-Friedrichs flux

$\widehat{f^{11}_{j+\frac{1}{2}}}(U)=r^{p}_{11}(x_{j+\frac{1}{2}})l_{p}^{(1)}(x_{j+\frac{1}{2}})\widehat{(A\cdot
U-\nu^{(1)}_{j+\frac{1}{2}}U)}_{{j+\frac{1}{2}}}.$

$\psi^{(2),\pm}_{j\pm\frac{1}{2}}=\psi^{(2)}(x^{(2),\pm}_{j\pm\frac{1}{2}}(t),t)$
and $\widehat{f^{12}}$ can be similarly defined at a cell boundary. We can
similarly obtain the ELDG scheme for $u^{2}_{h}$.

(4) Fully discrete ELDG scheme with method-of-lines RK schemes.

To update (3.28) from $U^{n}_{h}$ to $U^{n+1}_{h}$, we first apply the forward
Euler time discretization to get 1st order accuracy, then we generalize the
scheme to general RK methods. There are two main steps involved here. In order
to describe the implementation procedure of the fully discrete ELDG scheme, we
define the $L^{2}$ projection.

######  Definition 3.2.

($L^{2}$ projection) Let $u(x)\in L^{2}(\Omega)$, $M=\\{I_{j}\\}_{j=1}^{N}$
and $\tilde{M}=\\{\tilde{I}_{j}\\}_{j=1}^{N}$ be two different meshes of
$\Omega$. We have function spaces $V_{h}^{k}=\\{u(x):u|_{I_{j}}\in
P^{k}(I_{j}),\forall j\\}$ and
$\tilde{V}_{h}^{k}=\\{\tilde{u}(x):\tilde{u}|_{\tilde{I}_{j}}\in
P^{k}(\tilde{I}_{j}),\forall j\\}$ corresponding to meshes $M$ and
$\tilde{M}$. The $L^{2}$ projection of $u_{M}\in V_{h}^{k}$ onto space
$\tilde{V}_{h}^{k}$ can be defined as, find
$\tilde{u}_{\tilde{M}}(x)\in\tilde{V}_{h}^{k}$, s.t.

$\int_{\tilde{I}_{j}}\tilde{u}_{\tilde{M}}(x)\varphi(x)dx=\int_{\tilde{I}_{j}}u_{M}(x)\varphi(x)dx,\
\forall\varphi(x)\in\tilde{V}_{h}^{k}.$ (3.30)

We denote $\tilde{u}_{\tilde{M}}(x)=Proj[u_{M}(x);M,\tilde{M}]$. The
evaluation of the right hand side of (3.30) can be done in a subinterval-by-
subinterval fashion. The implementation details can be found in [16].

Then, we propose a fully discrete EL RK DG scheme with procedure as described:

  1. 1.

Obtain the initial condition on upstream meshes ${\tilde{I}^{(1)}_{j}(t^{n})}$
and ${\tilde{I}^{(2)}_{j}(t^{n})}$ of (3.28) by
$U^{(1)}_{h}(t^{n})=Proj[U^{n}_{h};I_{j},{\tilde{I}^{(1)}_{j}(t^{n})}]$ and
$U^{(2)}_{h}(t^{n})=Proj[U^{n}_{h};I_{j},{\tilde{I}^{(2)}_{j}(t^{n})}]$, which
are the $L^{2}$ projections of solutions from the background mesh to the
upstream mesh.

  2. 2.

Update (3.28) from $U^{n}_{h}$ to $U^{n+1}_{h}$, component-by-component.

    1. (a)

Get the mesh information of the dynamic element
$\tilde{I}^{(1)}_{j}(t^{(l)}),l=0,...,s$ on RK stages by
$\tilde{x}^{(1)}_{j\pm\frac{1}{2}}(t)=x_{j\pm\frac{1}{2}}+(t-t^{n+1})\nu^{(1)}_{j\pm\frac{1}{2}}$.
Here $s=1$ for forward-Euler method and $s=2$ for SSPRK2, see the blue domain
in Figure 3.2 and for general RK2 with intermediate stage in Figure 3.3.

    2. (b1)

For forward-Euler method, compute

$\displaystyle\int_{I_{j}}u^{1,n+1}_{h}\Psi(x)dx=\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}\psi^{(1),n}dx+\Delta
tL_{11}(U_{h}^{n},t^{n},\tilde{I}^{(1)}_{j}(t^{n}))$ (3.31) $\displaystyle+$
$\displaystyle\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}\psi^{(2),n}dx+\Delta
tL_{12}(U_{h}^{n},t^{n},\tilde{I}^{(2)}_{j}(t^{n}))$ $\displaystyle=$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}\psi^{(1),n}dx}+\Delta
t\left(-\widehat{f^{11}_{j+\frac{1}{2}}}(U_{h}^{(1)}(t^{n}))\psi^{(1),n,-}_{j+\frac{1}{2}}+\widehat{f^{11}_{j-\frac{1}{2}}}(U_{h}^{(1)}(t^{n}))\psi^{(1),n,+}_{j-\frac{1}{2}}\right.$
$\displaystyle\left.+\int_{I_{j}^{*,(1)}}f^{11}(U_{h}^{(1)}(t^{n}))\psi^{(1)}_{x}(x,t^{n})dx+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\int_{I_{j}^{*,(1)}}(r^{p}_{11}l_{p}^{(1)})_{x}AU^{n}_{h}\psi^{(1)}(x,t^{n})dx}\right.$
$\displaystyle\left.+\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}F(x,t^{n})\psi^{(1)}(x,t^{n})dx\right)$
$\displaystyle+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}\psi^{(2),n}dx}+\Delta
t\left(-\widehat{f^{12}_{j+\frac{1}{2}}}(U_{h}^{(2)}(t^{n}))\psi^{(2),n,-}_{j+\frac{1}{2}}+\widehat{f^{12}_{j-\frac{1}{2}}}(U_{h}^{(2)}(t^{n}))\psi^{(2),n,+}_{j-\frac{1}{2}}\right.$
$\displaystyle\left.+\int_{I_{j}^{*,(2)}}f^{12}(U_{h}^{(2)}(t^{n}))\psi^{(2)}_{x}(x,t^{n})dx+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\int_{I_{j}^{*,(2)}}(r^{p}_{12}l_{p}^{(2)})_{x}AU^{n}_{h}\psi^{(2)}(x,t^{n})dx}\right.$
$\displaystyle\left.+\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}F(x,t^{n})\psi^{(2)}(x,t^{n})dx\right)$

where $\Delta t^{n}=t^{n+1}-t^{n}$,
$U_{h}^{(1)}(t^{n})=Proj[U^{n}_{h},\\{I_{j}\\}_{j=1}^{N},\\{I_{j}^{*,(1)}\\}_{j=1}^{N}]$
and
$U_{h}^{(2)}(t^{n})=Proj[U^{n}_{h},\\{I_{j}\\}_{j=1}^{N},\\{I_{j}^{*,(2)}\\}_{j=1}^{N}]$.
We compute the four integration terms of
(3.31)$\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}\psi^{(1),n}dx$,
$\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}\psi^{(2),n}dx$,
$\int_{\tilde{I}^{(1)}_{j}(t^{n})}(r^{p}_{11}l_{p}^{(1)})_{x}AU^{n}_{h}\psi^{(1)}(x,t^{n})dx$
and
$\int_{\tilde{I}^{(2)}_{j}(t^{n})}(r^{p}_{12}l_{p}^{(2)})_{x}AU^{n}_{h}\psi^{(2)}(x,t^{n})dx$
highlighted in blue with $U_{h}^{n}$ on background meshes in order to makes
the scheme satisfy mass conservation. Thus they have to be evaluated
subinterval-by-subinterval since DG solution is discontinuous. In summary, we
can get $u_{h}^{1,n+1}$ on $I_{j}$ by (3.31). Similarly, we can get
$u_{h}^{2,n+1}$ on $I_{j}$.

0 |  |
---|---|---
1 | 1 |
| $\frac{1}{2}$ | $\frac{1}{2}$
Table 3.1: SSPRK2 Butcher Tableau 0 |  |
---|---|---
$\frac{1}{2}$ | $\frac{1}{2}$ |
| 0 | 1
Table 3.2: RK2 Butcher Tableau

    3. (b2)

For SSPRK2 method with Butcher tableau: 2, we get $u^{1}_{h}(t^{(1)})$ from
(3.31), then compute

$\displaystyle\int_{I_{j}}u^{1,n+1}_{h}\Psi(x)dx=\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}\psi^{(1),n}dx+0.5\Delta
tL_{11}(U_{h}(t^{(1)}),t^{(1)},I_{j})$ (3.32) $\displaystyle+$
$\displaystyle\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}\psi^{(2),n}dx+0.5\Delta
tL_{12}(U_{h}(t^{(1)}),t^{(1)},I_{j}),$

where $t^{(1)}=t^{n+1}$, $u^{1}_{h}(t^{(1)})$ and $u^{2}_{h}(t^{(1)})$ are
defined on background mesh $I_{j}$.

    4. (b3)

For general RK methods with intermediate stages, we will update intermediate
RK solutions on background mesh as in [11]. For example, for a 2nd order RK
method with Butcher tableau 2. It has an intermediate stage at
$t^{(1)}=t^{n+\frac{\Delta t}{2}}$, we propose the following steps as [11],
see Figure 3.3.

      1. i.

We denote the dynamic domain tracking $I_{j}$ from $t^{(1)}$ to $t_{n}$ with
speed $\nu^{(1)}_{j\pm\frac{1}{2}}$ at mesh point $x_{j\pm\frac{1}{2}}$ as
$\tilde{I}^{(1)}_{j,(1)}(t)$, see the green domain in Figure 3.3.
$\tilde{I}^{(2)}_{j,(1)}(t)$ related to second characteristic is defined
similarly. Then we can update $U^{(1)}_{h}$ on $\tilde{I}^{(1)}_{j,(1)}(t)$
and $\tilde{I}^{(2)}_{j,(1)}(t)$ from $t^{n}$ to $t^{(1)}$ as in a forward-
Euler method.

      2. ii.

We update $U_{h}^{n+1}$ on dynamic domain $\tilde{I}^{(1)}_{j}(t)$ and
$\tilde{I}^{(2)}_{j}(t)$ from $U_{h}^{n}$ with projection onto $I_{j}^{*,(1)}$
and $U_{h}(t^{(1)})$ with projection onto $\tilde{I}^{(1)}_{j,(1)}(t^{(1)})$.

${\tilde{I}^{(1)}_{j,(1)}}(t)$
$t^{n}$$t^{n+1}$$x_{j-\frac{1}{2}}$$x_{j+\frac{1}{2}}$$x^{\star,(1)}_{j-\frac{1}{2}}$
$x^{\star,(1)}_{j+\frac{1}{2}}$$\tilde{I}^{(1)}_{j}(t^{n})=I^{\star,(1)}_{j}$$\tilde{I}^{(1)}_{j}(t^{n+1})=I_{j}$$\tilde{I}^{(1)}_{j}(t^{(1)})$$\nu^{(1)}_{j-\frac{1}{2}}$
$\nu^{(1)}_{j+\frac{1}{2}}$ Figure 3.3: Update RK intermediate solution at
background mesh (red line) from the first characteristic family of a
hyperbolic system.

######  Theorem 3.3.

(Mass conservation) The proposed fully discrete ELDG scheme is locally mass
conservative. In particular, given a DG solution $u_{h}(x,t^{n})\in V_{h}^{k}$
with a periodic boundary condition, we have

$\sum_{i=1}^{N}\int_{I_{j}}U_{h}(x,t^{n+1})dx=\sum_{i=1}^{N}\int_{I_{j}}U_{h}(x,t^{n})dx.$

Proof. We firstly consider the forward Euler time discretization. Taking
$\Psi=1$ and $F=0$ in the scheme (3.31), we have

$\displaystyle\int_{\Omega}u^{1,n+1}_{h}dx=\sum_{j}\int_{I_{j}}u^{1,n+1}_{h}dx$
(3.33)
$\displaystyle=\sum_{j}\left[\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}dx+\Delta
t\left(-\widehat{f^{11}_{j+\frac{1}{2}}}(U_{h}^{(1)}(t^{n}))+\widehat{f^{11}_{j-\frac{1}{2}}}(U_{h}^{(1)}(t^{n}))+\int_{I_{j}^{*,(1)}}(r^{p}_{11}l_{p}^{(1)})_{x}AU_{h}^{n}dx\right)\right]$
$\displaystyle+\sum_{j}\left[\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}dx+\Delta
t\left(-\widehat{f^{12}_{j+\frac{1}{2}}}(U_{h}^{(2)}(t^{n}))+\widehat{f^{12}_{j-\frac{1}{2}}}(U_{h}^{(2)}(t^{n}))+\int_{I_{j}^{*,(2)}}(r^{p}_{12}l_{p}^{(2)})_{x}AU^{n}_{h}dx\right)\right]$
$\displaystyle=\sum_{j}\left(\int_{I_{j}^{*,(1)}}r^{p}_{11}l_{p}^{(1)}U^{n}_{h}dx+\int_{I_{j}^{*,(2)}}r^{p}_{12}l_{p}^{(2)}U^{n}_{h}dx\right)$
$\displaystyle+\Delta
t\sum_{j}\left(\int_{I_{j}^{*,(1)}}(r^{p}_{11}l_{p}^{(1)})_{x}AU_{h}^{n}dx+\int_{I_{j}^{*,(2)}}(r^{p}_{12}l_{p}^{(2)})_{x}AU^{n}_{h}dx\right)$
$\displaystyle=\int_{\Omega}(r^{p}_{11}l_{p}^{(1)}+r^{p}_{12}l_{p}^{(2)})U^{n}_{h}dx+\Delta
t\int_{\Omega}(r^{p}_{11}l_{p}^{(1)}+r^{p}_{12}l_{p}^{(2)})_{x}AU_{h}^{n}dx$
$\displaystyle=\int_{\Omega}u^{1,n}_{h}dx,$

which follows from the cancellation of unique fluxes at cell boundaries,
$r^{p}_{11}l_{p}^{(1)}+r^{p}_{12}l_{p}^{(2)}=[1,0]$ and (3.19) with
integration in a subinterval-by-subinterval fashion. The mass conservation for
the fully discrete ELDG scheme can be proved in a similar fashion.

######  Remark 3.4.

To maintain the mass-conservative property, the choice of eigenvectors
$R_{p}(x)$ is not necessarily exact for ELDG scheme, as long as $R_{p}(x)$ and
$R_{p}^{-1}(x)$ are a consistent pair throughout the domain.

###  3.2 2D linear wave equation

The solution for high-dimensional hyperbolic systems is given by means of a
characteristic cone, rather than individual characteristic lines. Numerically,
characteristic Galerkin [30] or evolution Galerkin [3, 22, 23] methods have
been proposed and developed to solve high dimensional hyperbolic system. This
method is constructed by taking into account information propagated in all
bicharacteristic directions and involving integrals around the characteristic
cone. However, the awkward integrals over the mantle, involving intermediate
times, limit both the accuracy and the stability of the resulting schemes.
Thus the finite volume evolution Galerkin (FVEG) schemes are introduced, which
is in a predictor-corrector plus finite volume framework to get higher
accuracy. FVEG method has been developed in [28] and widely applied in [21,
28, 24], and the stability and accuracy have been investigated in [25, 28, 27,
26]. Even though the FVEG method can achieve high-order accuracy and stability
with extra large step, the algorithm implementation is still very complex for
high-dimensional problems. In this paper, we use the splitting method for
higher dimensional problem.

Consider a linear 2D wave equation

$u_{tt}=(a(x,y,t)u_{x})_{x}+(b(x,y,t)u_{y})_{y},(x,y)\in\Omega.$ (3.34)

It can be written as a first order linear hyperbolic system

$U_{t}+(A(x,y,t)U)_{x}+(B(x,y,t)U)_{y}=0,(x,y)\in\Omega,$ (3.35)

where

$U=\begin{bmatrix}u_{t}\\\ u_{x}\\\ u_{y}\end{bmatrix},\qquad
A=\begin{bmatrix}0&-a(x,y,t)&0\\\ -1&0&0\\\ 0&0&0\end{bmatrix},\qquad
B=\begin{bmatrix}0&0&-b(x,y,t)\\\ 0&0&0\\\ -1&0&0\end{bmatrix}.$

The domain $\Omega$ is partitioned into rectangular meshes with each
computational cell
$\Omega_{ij}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j+\frac{1}{2}},y_{j+\frac{1}{2}}]$,
where we use the piecewise $Q^{k}$ tensor-product polynomial spaces. Then we
extend ELDG algorithm to 2D problems via dimensional splitting [31].

Figure 3.4: Illustration of the 2D ELDG scheme via Strang splitting. k = 3.

  1. 1.

We first locate $(k+1)^{2}$ tensor-product Gaussian nodes on cell
$\Omega_{ij}:(x_{i,p},y_{j,q}),\ p,q=0,...,k$. For example, see Figure 3.4
(left) for the case of $k=3$.

  2. 2.

Then, the equation (3.35) is split into two 1D hyperbolic problems based on
the quadrature nodes in both $x-$ and $y-$ directions:

$\displaystyle U_{t}+(A(x,y,t)U)_{x}$ $\displaystyle=0,$ (3.36) $\displaystyle
U_{t}+(B(x,y,t)U)_{y}$ $\displaystyle=0.$ (3.37)

Based on a 1D ELDG formulation, the split equations (3.36) and (3.37) are
evolved via Strang splitting over a time step $\Delta t$ as follows.

     * •

Evolve 1D equation (3.36) at different $y_{j,q}^{\prime}s$ for a half time-
step $\Delta t/2$, see Figure 3.4 (middle). For each $y_{j,q}$, the $(k+1)$
point values are mapped to a $P^{k}$ polynomial per cell, then the 1D system
(3.36) is evolved by the proposed ELDG scheme. Finally, we can map the evolved
$P^{k}$ polynomial back to the $(k+1)$ point values to update the solution.

     * •

Evolve 1D system (3.37) at different $x_{i,p}^{\prime}s$ for a full time-step
$\Delta t$ as above, see Figure 3.4 (right).

     * •

Evolve 1D system (3.36) at different $y_{j,q}^{\prime}s$ for another half
time-step $\Delta t/2$.

The splitting 2D ELDG formulation maintains high order accuracy in space,
extra large time stepping size with stability and mass conservation; and has a
second order splitting error.

##  4 Numerical results

In this section, we show numerical results of the proposed scheme for wave
equations. We set the time stepping size as $\Delta t=\frac{CFL}{a}\Delta x$
for 1D and $\Delta t=\frac{CFL}{\frac{a}{\Delta x}+\frac{b}{\Delta y}}$, where
a and b are maximum eigenvalues of coefficient matrixes in $x$\- and
$y$-directions respectively. We use RK4 [18] for time discretization. We study
the following aspects: the spatial order of convergence by using small enough
time stepping sizes, the spatial super-convergence of the post-processed
solutions [2, 9] produced by convolving the ELDG solution with a suitable
kernel consisting of B-splines for the purpose of improving spatial
convergence rate, the temporal order of convergence and numerical stability
under a large time stepping size by varying CFL for a fixed spatial mesh.

###  4.1 1D wave equations

######  Example 4.1.

(1D wave equation with constant coefficient.) We consider the 1D wave equation
eq. (3.1) with constant coefficient $a(x)=1$ and the source term $f(x,t)=0$.
The initial data is $u(x,0)=\sin(x),x\in[0,2\pi]$ with periodic boundary
condition. The exact solution is $u(x,t)=\sin(x+t)$. For the constant
coefficient problem, if using exact characteristic velocity fields for space-
time partition and exact eigenvectors, the proposed ELDG method is the same as
SL DG, then it is unconditionally stable. Here we perturb the characteristic
velocity $\nu_{j+\frac{1}{2}}^{(1)}$ in (3.6) at cell boundaries and/or
$a_{p}(x)$ in (3.17) related to approximating eigenvectors to get ELDG, ELDG1,
ELDG2 and ELDG3 schemes respectively. Similarly we implement the non mass-
conservative ELDG methods denoted as NMC ELDG, NMC ELDG1, NMC ELDG2 and NMC
ELDG3. Related parameters of these ELDG methods are given in Table 4.3.

Table 4.3: The numerical parameters of ELDG, ELDG1, ELDG2, ELDG3 methods and NMC ELDG, NMC ELDG1, NMC ELDG2, NMC ELDG3 methods for $u_{tt}=u_{xx}$.  setting | ELDG | ELDG1 | ELDG2 | ELDG3
---|---|---|---|---
$\nu^{(1)}_{j+\frac{1}{2}}$ in (3.6), $\nu^{(2)}_{j+\frac{1}{2}}=-\nu^{(1)}_{j+\frac{1}{2}}$ | $1$ | $1+\sin(x_{j})\Delta x$ | $1$ | $1+\sin(x_{j})\Delta x$
$a_{p}(x)$ in (3.17) | $1$ | $1$ | $1+\sin(x)\Delta x$ | $1+\sin(x)\Delta x$
setting | NMC ELDG | NMC ELDG1 | NMC ELDG2 | NMC ELDG3
$\nu^{(1)}_{j+\frac{1}{2}}$ in (3.6), $\nu^{(2)}_{j+\frac{1}{2}}=-\nu^{(1)}_{j+\frac{1}{2}}$ | $1$ | $1+\sin(x_{j})\Delta x$ | $1$ | $1+\sin(x_{j})\Delta x$
$a_{j}$ in (3.10) | $1$ | $1$ | $1+\sin(x_{j})\Delta x$ | $1+\sin(x_{j})\Delta x$

Table 4.4 and 4.5 report spatial accuracies of the ELDG, ELDG1, ELDG2 and
ELDG3 methods for this example under the same time stepping size without and
with post-processing technique [2, 9]. We can observe the optimal convergence
rate $k+1$ and $2k+1$. We vary time stepping sizes, with fixed well-resolved
spatial meshes, and plot error vs. $CFL$ in Figure 4.5 and 4.6 for ELDG,
ELDG1, ELDG2 and ELDG3 schemes without and with post-processed technique
respectively, after a long time $T=100$. The plots from post-processed ELDG
schemes better show the fourth order temporal convergence. ELDG2 and ELDG3
perform comparably; they have a more restricted time step constraint than
ELDG1. It indicates that, stability is affected by approximations of
characteristic via the space-time partition and approximation of eigenvectors.
We also note that, in both Figure 4.5 and 4.6, the CFL allowed with stability
is much larger than that of the RK DG method which is $\frac{1}{2k+1}$.
Further, we verify the mass conservative property of the ELDG schemes are
around machine precision and the non-mass conservative property of the NMC
ELDG schemes is presented in Figure 4.7.

Table 4.4: 1D wave equation with constant coefficient. $u_{tt}=u_{xx}$ with initial condition $u(x,0)=\sin(x)$. $T=1$. We use $CFL=0.3$ and $CFL=0.18$ with RK4 time discretization for all $P^{1}$ and $P^{2}$ respectively. The error for only $u^{1}=u_{t}$ was shown in this table. Mesh |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
---|---|---|---|---|---|---|---|---
|  $P^{1}$ ELDG |  $P^{1}$ ELDG1 |  $P^{1}$ ELDG2 |  $P^{1}$ ELDG3
20 | 2.54E-03 | – | 2.42E-03 | – | 2.55E-03 | – | 2.49E-03 | –
40 | 6.18E-04 | 2.03 | 5.97E-04 | 2.02 | 6.18E-04 | 2.04 | 5.99E-04 | 2.06
80 | 1.58E-04 | 1.96 | 1.55E-04 | 1.94 | 1.58E-04 | 1.96 | 1.55E-04 | 1.95
160 | 3.66E-05 | 2.11 | 3.62E-05 | 2.10 | 3.66E-05 | 2.11 | 3.62E-05 | 2.10
|  $P^{2}$ ELDG |  $P^{2}$ ELDG1 |  $P^{2}$ ELDG2 |  $P^{2}$ ELDG3
20 | 5.92E-05 | – | 6.91E-05 | – | 6.01E-05 | – | 7.02E-05 | –
40 | 7.48E-06 | – | 7.83E-06 | 3.14 | 7.49E-06 | 3.00 | 7.81E-06 | 3.17
80 | 9.17E-07 | 3.03 | 9.29E-07 | 3.08 | 9.17E-07 | 3.03 | 9.29E-07 | 3.07
160 | 1.17E-07 | 2.97 | 1.18E-07 | 2.98 | 1.17E-07 | 2.97 | 1.18E-07 | 2.98
Table 4.5: 1D wave equation with constant coefficient. $u_{tt}=u_{xx}$ with initial condition $u(x,0)=\sin(x)$. $T=1$. We use $CFL=0.3$ and $CFL=0.18$ with RK4 time discretization and post-processed technique for all $P^{1}$ and $P^{2}$ respectively. The error with post-processed technique for only $u^{1}=u_{t}$ was shown in this table. Mesh |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order |  $L^{1}$ error | Order
---|---|---|---|---|---|---|---|---
|  $P^{1}$ ELDG |  $P^{1}$ ELDG1 |  $P^{1}$ ELDG2 |  $P^{1}$ ELDG3
20 | 2.26E-04 | – | 2.49E-04 | – | 2.38E-04 | – | 2.39E-04 | –
40 | 2.36E-05 | 3.26 | 2.40E-05 | 3.38 | 2.40E-05 | 3.31 | 2.35E-05 | 3.35
80 | 2.66E-06 | 3.15 | 2.64E-06 | 3.18 | 2.67E-06 | 3.16 | 2.62E-06 | 3.16
160 | 3.15E-07 | 3.08 | 3.11E-07 | 3.08 | 3.15E-07 | 3.09 | 3.11E-07 | 3.08
|  $P^{2}$ ELDG |  $P^{2}$ ELDG1 |  $P^{2}$ ELDG2 |  $P^{2}$ ELDG3
20 | 2.15E-06 | – | 2.27E-06 | – | 2.19E-06 | – | 2.28E-06 | –
40 | 3.63E-08 | 5.89 | 3.86E-08 | 5.87 | 3.67E-08 | 5.90 | 3.89E-08 | 5.87
80 | 6.40E-10 | 5.83 | 6.79E-10 | 5.83 | 6.46E-10 | 5.83 | 6.84E-10 | 5.83
160 | 1.27E-11 | 5.66 | 1.33E-11 | 5.68 | 1.28E-11 | 5.66 | 1.34E-11 | 5.68

Figure 4.5: The $L^{\infty}$ error versus $CFL$ of ELDG methods, ELDG1, ELDG2
and ELDG3 methods for 1D wave equation with constant coefficient:
$u_{tt}=u_{xx}$ with initial condition $u(x,0)=\sin(x)$. A long time
simulation is performed with $T=100$ and mesh size $N=160$.

Figure 4.6: The $L^{\infty}$ error versus $CFL$ of ELDG methods, ELDG1, ELDG2
and ELDG3 methods with post-processed technique for 1D wave equation with
constant coefficient: $u_{tt}=u_{xx}$ with initial condition $u(x,0)=\sin(x)$.
A long time simulation is performed with $T=100$ and mesh size $N=160$.

Figure 4.7: The error of mass versus time of ELDG2, ELDG3, NMC ELDG2 and NMC
ELDG3 methods for 1D wave equation with constant coefficient: $u_{tt}=u_{xx}$
with initial condition $u(x,0)=\sin(x)$. A long time simulation $T=100$ is
performed with meshes $N=160$, $CFL=0.1$ and RK4 time discretization.

######  Example 4.2.

Then, we test the ELDG schemes for eq. (3.1) with a Gaussian initial condition

$u^{1}=exp(-\frac{x^{2}}{0.005}),u^{2}=0.$

The computational domain is $[-1,1]$ with the periodic boundary conditions.
The exact solutions
$u^{1}=0.5(exp(-\frac{(x+t)^{2}}{0.005})+exp(-\frac{(x-t)^{2}}{0.005}))$ and
$u^{2}=0.5(exp(-\frac{(x+t)^{2}}{0.005})-exp(-\frac{(x-t)^{2}}{0.005}))$ are
the superposition of two Gaussian functions with a periodic extensions. We
plot $u_{1}$ ELDG3 with $P^{1}$ and $P^{2}$ numerical solutions at time
$T=50.5$ for this system in Figure 4.8. We can observed that there is no
significant phase difference with a long time simulation, meanwhile the
dissipation can be improved by the mesh refinement and higher order spatial
approximation.

Figure 4.8: Plots of the exact and numerical solutions $u_{1}$ at time
$T=50.5$ of ELDG3 scheme for solving $u_{tt}=u_{xx}$ with Gaussian function
initial condition. The mesh size of $N=80$ and $N=320$ are used. Left: $k=1$
ELDG3 with $CFL=1.5$. Right: $k=2$ ELDG3 with $CFL=0.9$.

######  Example 4.3.

Next, we test the ELDG schemes for eq. (3.1) on $[0,2\pi]$ with the periodic
boundary conditions and the following discontinuous initial condition

$\displaystyle u^{1}_{0}(x)=\left\\{\begin{aligned} &1.,&if\ 0.95\pi\leq x\leq
1.05\pi,\\\ &0.5,\ &otherwise,\end{aligned}\right.$ (4.1) $\displaystyle
u^{2}_{0}(x)=1.$

The exact solutions $u^{1}$ and $u^{2}$ are discontinuous piecewise constants
with moving discontinuities. It is a challenging test for controlling
oscillations around discontinuities. We adopt a simple TVD limiter on
background mesh at each RK stages with M = 0 in [10] for all schemes. As shown
in Figure 4.5, the $CFL$ constraint with stability is slightly less than 1 for
ELDG3 scheme. We plot the numerical solutions $u_{1}$ of ELDG3 scheme with
$P^{1}$ and $P^{2}$, $CFL=0.9$ in Figure 4.9. It is found that oscillations
are well controlled with the TVD limiter and ELDG method performs well for
large time stepping size. Moreover, we track the mass conservation of ELDG
methods VS NMC ELDG methods for eq. (3.1) and present results in Figure 4.10.
It shows that the ELDG schemes maintain the mass conservation at the level of
machine error, while the NMC ELDG schemes do not.

Figure 4.9: Plots of the numerical solutions $u_{1}$ of ELDG3 scheme with TVD
limiter for solving $u_{tt}=u_{xx}$ with step function initial condition. The
final integration time T is 2.85. The mesh of $160$ is used. Left: $k=1$
ELDG3+TVDlimiter with $CFL=1.5$. Right: $k=2$ ELDG3+TVDlimiter with $CFL=0.9$.

Figure 4.10: The error of mass versus time of ELDG2, ELDG3, NMC ELDG2 and NMC
ELDG3 methods with TVD limiter for 1D wave equation with constant coefficient:
$u_{tt}=u_{xx}$ with initial condition step function. $T=2.85$, $N=160$,
$CFL=0.9$ and RK4 time discretization are performed for the simulation.

######  Example 4.4.

(1D wave equation with variable coefficient and source term.) We consider the
1D wave equation eq. (3.1) with variable coefficient $a(x)=2+\sin(x)$ and
exact solution $u(x,t)=\sin(x-2t)$ is periodic on $[0,2\pi]$. The source term
is
$f(x,t)=-4\sin(x-2t)+\sin(x-2t)(2+\sin(x))^{2}-2(2+\sin(x))\cos(x)\cos(x-2t)$.
For computation, we choose mesh velocity
$\nu^{(1)}_{j+\frac{1}{2}}=a(x_{j+\frac{1}{2}}),\nu^{(2)}_{j+\frac{1}{2}}=-\nu^{(1)}_{j+\frac{1}{2}}$
and exact eigenvectors with $a_{p}(x)=2+\sin(x)$ in (3.17).

Table 4.6: 1D wave equation with variable coefficient and source term. $u_{tt}=(a^{2}(x)u_{x})_{x}+f(x,t)$ with initial condition $u(x,0)=\cos(x)$. $T=1$. We use $CFL=0.1$ for $P^{1}$ and $P^{2}$ with RK4 time discretization. The error for only $u^{1}=u_{t}$ was shown in this table. Mesh |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
| $P^{1}$
20 | 5.20E-03 | – | 6.70E-03 | – | 2.39E-02 | –
40 | 1.32E-03 | 1.98 | 1.74E-03 | 1.94 | 6.54E-03 | 1.87
80 | 3.28E-04 | 2.00 | 4.40E-04 | 1.99 | 1.68E-03 | 1.96
160 | 8.11E-05 | 2.02 | 1.09E-04 | 2.01 | 4.16E-04 | 2.01
| $P^{2}$
20 | 1.16E-04 | – | 1.61E-04 | – | 5.48E-04 | –
40 | 1.48E-05 | 2.97 | 2.01E-05 | 3.00 | 6.70E-05 | 3.03
80 | 1.88E-06 | 2.98 | 2.46E-06 | 3.03 | 7.84E-06 | 3.10
160 | 2.31E-07 | 3.02 | 3.13E-07 | 2.98 | 1.04E-06 | 2.91
Table 4.7: 1D wave equation with variable coefficient and source term. $u_{tt}=(a^{2}(x)u_{x})_{x}+f(x,t)$ with initial condition $u(x,0)=\cos(x)$. $T=1$. We use $CFL=0.1$ for $P^{1}$ and $P^{2}$ with post-processed technique and RK4. The error for only $u^{1}=u_{t}$ was shown in this table. Mesh |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
| $P^{1}$
20 | 9.10E-04 | – | 1.07E-03 | – | 1.91E-03 | –
40 | 1.07E-04 | 3.09 | 1.26E-04 | 3.09 | 2.26E-04 | 3.08
80 | 1.29E-05 | 3.05 | 1.53E-05 | 3.05 | 2.75E-05 | 3.04
160 | 1.58E-06 | 3.02 | 1.88E-06 | 3.02 | 3.39E-06 | 3.02
| $P^{2}$
20 | 5.34E-06 | – | 6.39E-06 | – | 1.56E-05 | –
40 | 8.95E-08 | 5.90 | 1.03E-07 | 5.96 | 2.82E-07 | 5.79
80 | 1.73E-09 | 5.69 | 1.94E-09 | 5.73 | 3.32E-09 | 6.41
160 | 6.62E-11 | 4.71 | 7.59E-11 | 4.67 | 1.46E-10 | 4.51

Figure 4.11: The $L^{\infty}$ error versus $CFL$ of ELDG method for
$u_{tt}=(a^{2}(x)u_{x})_{x}+f(x,t)$ with initial condition $u(x,0)=\cos(x)$.
$T=1$. $\Delta t=CFL\Delta x$.

Figure 4.12: The $L^{\infty}$ error versus $CFL$ of ELDG method with post-
processed technique for $u_{tt}=(a^{2}(x)u_{x})_{x}+f(x,t)$ with initial
condition $u(x,0)=\cos(x)$. $T=1$. $\Delta t=CFL\Delta x$.

The expected optimal spatial accuracies of the ELDG methods without and with
post-processing technique are shown in Table 4.6 and Table 4.7 respectively.
In Figure 4.11 and 4.12, we plot the $L^{\infty}$ error versus $CFL$ of ELDG
methods without and with post-processing technique respectively. The following
observations are made: (1) The high order accuracy of the RK method reduce the
error magnitude when large time stepping size is used; (2) The ELDG methods
with RK4 time discretization perform well around and before $CFL=1$, which is
well above the stability constraint of the RK DG method $1/(2k+1)$ for $P^{k}$
approximations. (3) After $CFL=1$ and before stability constraint of the
method, the temporal convergence order is observed to be consistent with the
order of RK discretization; (4) The ELDG methods with post-processing
technique have smaller error magnitute than those without post-processing.

######  Example 4.5.

(Two-dimensional linear system with constant coefficient matrices.) The second
order wave equation $u_{tt}=u_{xx}+u_{yy}$ can be written as the following
first order linear hyperbolic system:

$\begin{cases}\left(\begin{matrix}u\\\
v\end{matrix}\right)_{t}+\left(\begin{matrix}-1&0\\\
0&1\end{matrix}\right)\left(\begin{matrix}u\\\
v\end{matrix}\right)_{x}+\left(\begin{matrix}0&-1\\\
-1&0\end{matrix}\right)\left(\begin{matrix}u\\\
v\end{matrix}\right)_{y}=\left(\begin{matrix}0\\\ 0\end{matrix}\right),\\\
u(x,y,0)=\frac{1}{2\sqrt{2}}\sin(x+y)-\frac{1}{2\sqrt{2}}\cos(x+y),\\\
v(x,y,0)=\frac{\sqrt{2}-1}{2\sqrt{2}}\sin(x+y)+\frac{\sqrt{2}+1}{2\sqrt{2}}\cos(x+y)\end{cases}$
(4.2)

with period boundary conditions in both $x$ and $y$ directions. The exact
solution is

$\begin{cases}u(x,y,t)=\frac{1}{2\sqrt{2}}\sin(x+y+\sqrt{2}t)-\frac{1}{2\sqrt{2}}\cos(x+y-\sqrt{2}t),\\\
v(x,y,t)=\frac{\sqrt{2}-1}{2\sqrt{2}}\sin(x+y+\sqrt{2}t)+\frac{\sqrt{2}+1}{2\sqrt{2}}\cos(x+y-\sqrt{2}t).\end{cases}$
(4.3)

We notice that the two matrices in equation (4.2) don’t commute, thus the
linear system can not be reduced to 2-D scalar problems. We test accuracy for
$Q^{k}$ ELDG methods with RK4 and 4th splitting method [5] at $T=1$ for
$k=1,2$ with $CFL=0.1$ in Table 4.8. As expected, the $(k+1)$th order
convergence is observed for these methods. We plot the $L^{\infty}$ error
versus $CFL$ of ELDG methods with $Q^{1}$ (left) and $Q^{2}$ (right)
polynomial spaces for this case with Strang splitting and 4th order splitting
in Figure 4.13, which shows that second and forth order splitting errors are
dominant when time-stepping sizes are large enough. The CFL constraint with
stability for ELDG method is larger than that for general RK DG method when
high order time discretization is applied.

Table 4.8: Two-dimensional linear system with constant coefficient matrices. $Q^{k}$ EL DG methods (k = 1, 2) with RK4 and 4th splitting time discretization methods for (4.2) with the smooth initial condition at T = 1. $CFL=0.1$.  Mesh |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
| $Q^{1}$
$20^{2}$ | 8.03E-04 | – | 9.47E-04 | – | 1.85E-03 | –
$40^{2}$ | 2.16E-04 | 1.89 | 2.50E-04 | 1.92 | 4.56E-04 | 2.02
$80^{2}$ | 5.57E-05 | 1.96 | 6.40E-05 | 1.97 | 1.13E-04 | 2.02
$160^{2}$ | 1.43E-05 | 1.96 | 1.64E-05 | 1.97 | 2.84E-05 | 1.99
| $Q^{2}$
$20^{2}$ | 1.70E-04 | – | 1.90E-04 | – | 3.12E-04 | –
$40^{2}$ | 2.21E-05 | 2.95 | 2.47E-05 | 2.94 | 4.14E-05 | 2.91
$80^{2}$ | 2.75E-06 | 3.00 | 3.08E-06 | 3.00 | 5.21E-06 | 2.99
$160^{2}$ | 3.38E-07 | 3.02 | 3.80E-07 | 3.02 | 6.45E-07 | 3.01

Figure 4.13: The $L^{\infty}$ error versus $CFL$ of ELDG method with Strang
splitting and 4th splitting, RK4 time discretization for (4.2). $T=1$.

######  Example 4.6.

(Two-dimensional linear system with variable coefficient matrices.) The second
order wave equation $u_{tt}=(a^{2}(x,y)u_{x})_{x}+(b^{2}(x,y)u_{y})_{y}$ can
be written as the following first order linear hyperbolic system:

$\begin{cases}\left(\begin{matrix}u_{1}\\\ u_{2}\\\
u_{3}\end{matrix}\right)_{t}+\left[\left(\begin{matrix}0&-(a(x,y))^{2}&0\\\
-1&0&0\\\ 0&0&0\end{matrix}\right)\left(\begin{matrix}u_{1}\\\ u_{2}\\\
u_{3}\end{matrix}\right)\right]_{x}+\left[\left(\begin{matrix}0&0&-(b(x,y))^{2})\\\
0&0&0\\\ -1&0&0\end{matrix}\right)\left(\begin{matrix}u_{1}\\\ u_{2}\\\
u_{3}\end{matrix}\right)\right]_{y}=\left(\begin{matrix}0\\\ 0\\\
0\end{matrix}\right).\end{cases}$ (4.4)

The simplest case is to choose $a(x,y)=b(x,y)=1$, computational domain
$[-1,1]\times[-1,1]$ with periodic boundary condition and the Gaussian
function initial condition:

$\begin{cases}u_{1}(x,y,0)=exp(-\frac{x^{2}+y^{2}}{0.005}),\\\
u_{2}(x,y,0)=0,\\\ u_{3}(x,y,0)=0.\end{cases}$ (4.5)

For this example, we show the numerical ELDG $Q^{2}$ solution $u_{1}$ at times
$T=0.5;1;1.5;2$ in Figure 4.14.

Figure 4.14: Plots of the ELDG numerical solutions $u_{1}$ and their contour
plots at $T=0.5;1;1.5;2$ for 2D wave equation $u_{tt}=u_{xx}+u_{yy}$ with
Gaussian function initial condition. The mesh of $80\times 80$ is used with
4th splitting method and RK4 time discretization, $CFL=1$.

Next we consider the system (4.4) with the initial condition

$\begin{cases}u_{1}(x,y,0)=2\cos(x+y),\\\ u_{2}(x,y,0)=\cos(x+y),\\\
u_{3}(x,y,0)=\cos(x+y),\end{cases}$ (4.6)

where $a(x,y)=1+0.5\sin(x+y)$, $b(x,y)=\sqrt{(4-(1+0.5\sin(x+y))^{2})}$ and
the boundary condition is periodic in both $x$ and $y$ directions. The exact
solution is

$\begin{cases}u_{1}(x,y,t)=2\cos(x+y+2t),\\\ u_{2}(x,y,t)=\cos(x+y+2t),\\\
u_{3}(x,y,t)=\cos(x+y+2t).\end{cases}$ (4.7)

We report the spatial accuracy of $Q^{k}$ ELDG methods in Table 4.9. The
expected optimal convergence is observed. We plot the $L^{\infty}$ error
versus $CFL$ of ELDG methods in Figure 4.15. The ELDG methods perform as well
as that for the linear system with constant coefficient matrices, and the CFL
allowed with stability is much larger than that of the RK DG method.

Table 4.9: Two-dimensional linear system with variable matrices. $Q^{k}$ ELDG methods (k = 1, 2) with RK4 and 4th splitting time discretization methods for (4.4) with the smooth initial condition at T = 0.1, $CFL=0.1$.  Mesh |  $L^{1}$ error | Order |  $L^{2}$ error | Order |  $L^{\infty}$ error | Order
---|---|---|---|---|---|---
| $Q^{1}$
$20^{2}$ | 1.89E-03 | – | 2.46E-03 | – | 5.29E-03 | –
$40^{2}$ | 4.58E-04 | 2.05 | 6.01E-04 | 2.03 | 1.26E-03 | 2.07
$80^{2}$ | 1.13E-04 | 2.02 | 1.50E-04 | 2.00 | 3.10E-04 | 2.02
$160^{2}$ | 2.81E-05 | 2.01 | 3.73E-05 | 2.01 | 7.67E-05 | 2.01
| $Q^{2}$
$20^{2}$ | 2.95E-04 | – | 3.66E-04 | – | 8.43E-04 | –
$40^{2}$ | 4.06E-05 | 2.86 | 4.91E-05 | 2.90 | 1.01E-04 | 3.07
$80^{2}$ | 5.15E-06 | 2.98 | 6.20E-06 | 2.99 | 1.30E-05 | 2.95
$160^{2}$ | 6.49E-07 | 2.99 | 7.80E-07 | 2.99 | 1.63E-06 | 2.99

Figure 4.15: The $L^{\infty}$ error versus $CFL$ of ELDG method with Strang
splitting and 4th splitting, RK4 time discretization for (4.4). $T=1$, mesh
size$40^{2}$.

##  5 Conclusion

In this paper, we have developed a mass-conservative Eulerian-Lagrangian
discontinuous Galerkin (ELDG) method for wave equations written in the form of
hyperbolic systems. The new framework track the information of each
characteristic family by the corresponding characteristic region, and combine
in a mass-conservative fashion; The method inherits advantages in stability
under large time stepping sizes, and in mass conservation, compactness and
high order accuracy. These advantages are numerically verified by extensive
numerical tests for 1D and 2D linear wave equations. Future works include
further theoretic development and application of ELDG methods for nonlinear
hyperbolic problems.

##  6 Appendix (The notation of test function)

We give some notations which is used in the implementation of the fully
discrete ELDG scheme with RK time discretization. For convenience, we only
give the definitions related to $\Omega_{j}^{(1)}$ below because we can
similarly get the definitions related to $\Omega_{j}^{(2)}$. As scalar case,
we take test function $\psi^{(1)}(x,t)$ as
$\psi^{(1)}_{j,m}(x,t)=\Psi_{j,m}(x-\alpha^{(1)}(t-t^{n+1}))$ in the adjoint
problem which is a set of basis of $P^{k}(\tilde{I}^{(1)}_{j}(t))$, where
$\tilde{I}^{(1)}_{j}(t)$ is donated by a domain related $\lambda^{(1)}$ as
$\tilde{I}_{j}(t)$ in scalar problem. Here, we also take $\Psi_{j,m}$ as
orthogonal basis on $I_{j}$, and let

$\displaystyle
u_{h}^{1,(1)}(x,t)=\sum_{l=0}^{k}\hat{u}^{1,(1);(l)}_{j}(t)\psi^{(1)}_{j,l}(x,t),\quad\mbox{on}\quad\tilde{I}^{(1)}_{j}(t),$
(6.1)

where $\hat{u}^{1,(1);(l)}$ are coefficients for the basis. Let
$\hat{U}^{1,(1)}_{j}(t)=(\hat{u}_{j}^{1,(1);(0)}(t),\cdots,\hat{u}_{j}^{1,(1);(k)}(t))^{T}$
be the coefficient vector of size $(k+1)\times 1$. Then we have

$[\int_{\tilde{I}^{(1)}_{j}(t)}u_{h}^{1,(1)}(x,t)\psi^{(1)}_{j,0}(x,t)dx,...,\int_{\tilde{I}^{(1)}_{j}(t)}u_{h}^{1,(1)}(x,t)\psi^{(1)}_{j,k}(x,t)dx]^{T}=\hat{U}^{1,(1)}_{j}(t),\forall
t\in[t^{n},t^{n+1}].$

$\hat{U}^{2,(1)}_{j}(t)$ can be similarly defined. Similar definition can be
made to $\Omega^{(2)}_{j}$ and $\psi^{(2)}(x,t)$ for the second
characteristics family.

## References

  * [1] M. Baccouch.  A local discontinuous Galerkin method for the second-order wave equation.  Computer Methods in Applied Mechanics and Engineering, 209:129–143, 2012.
  * [2] J. H. Bramble and A. H. Schatz.  Higher order local accuracy by averaging in the finite element method.  Mathematics of Computation, 31(137):94–111, 1977.
  * [3] D. S. Butler.  The numerical solution of hyperbolic systems of partial differential equations in three independent variables.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 255(1281):232–252, 1960.
  * [4] X. Cai, W. Guo, and J.-M. Qiu.  A high order conservative semi-Lagrangian discontinuous Galerkin method for two-dimensional transport simulations.  Journal of Scientific Computing, 73(2-3):514–542, 2017.
  * [5] X. Cai, J. Qiu, and J.-M. Qiu.  A conservative semi-Lagrangian HWENO method for the Vlasov equation.  Journal of Computational Physics, 323:95–114, 2016.
  * [6] X. Cai, J.-M. Qiu, and Y. Yang.  An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics.  Journal of Computational Physics, 439:110392, 2021.
  * [7] M. Celia, T. Russell, I. Herrera, and R. Ewing.  An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation.  Advances in water resources, 13(4):187–206, 1990.
  * [8] C.-S. Chou, C.-W. Shu, and Y. Xing.  Optimal energy conserving local discontinuous Galerkin methods for second-order wave equation in heterogeneous media.  Journal of Computational Physics, 272:88–107, 2014.
  * [9] B. Cockburn, M. Luskin, C.-W. Shu, and E. Suli.  Enhanced accuracy by post-processing for finite element methods for hyperbolic equations.  Mathematics of Computation, 72(242):577–606, 2003.
  * [10] B. Cockburn and C.-W. Shu.  TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework.  Mathematics of computation, 52(186):411–435, 1989.
  * [11] M. Ding, X. Cai, W. Guo, and J.-M. Qiu.  A semi-Lagrangian discontinuous Galerkin (DG)–local DG method for solving convection-diffusion equations.  Journal of Computational Physics, 409:109295, 2020.
  * [12] D. R. Durran.  Numerical methods for wave equations in geophysical fluid dynamics, volume 32.  Springer Science & Business Media, 2013.
  * [13] R. S. Falk and G. R. Richter.  Explicit finite element methods for symmetric hyperbolic equations.  SIAM Journal on Numerical Analysis, 36(3):935–952, 1999.
  * [14] M. J. Grote, A. Schneebeli, and D. Schötzau.  Discontinuous Galerkin finite element method for the wave equation.  SIAM Journal on Numerical Analysis, 44(6):2408–2431, 2006.
  * [15] M. J. Grote and D. Schötzau.  Optimal error estimates for the fully discrete interior penalty DG method for the wave equation.  Journal of Scientific Computing, 40(1):257–272, 2009.
  * [16] W. Guo, R. Nair, and J.-M. Qiu.  A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere.  Monthly Weather Review, 142(1):457–475, 2013.
  * [17] X. Hong and Y. Xia.  Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for hyperbolic equations involving $\delta$-singularities.  SIAM Journal on Numerical Analysis, 58(1):125–152, 2020.
  * [18] G.-S. Jiang and C.-W. Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [19] N. A. Kampanis, V. Dougalis, and J. A. Ekaterinaris.  Effective computational methods for wave propagation.  CRC Press, 2008.
  * [20] C. Klingenberg, G. Schnücke, and Y. Xia.  Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for conservation laws: analysis and application in one dimension.  Mathematics of Computation, 86(305):1203–1232, 2017.
  * [21] T. Kröger and M. Lukáčová-Medvidová.  An evolution Galerkin scheme for the shallow water magnetohydrodynamic equations in two space dimensions.  Journal of Computational Physics, 206(1):122–149, 2005.
  * [22] P. Lax.  Systems of conservation laws.  Technical report, LOS ALAMOS NATIONAL LAB NM, 1959.
  * [23] M. Lukáčová-Medvidová, K. Morton, and G. Warnecke.  Evolution Galerkin methods for hyperbolic systems in two space dimensions.  Mathematics of Computation, 69(232):1355–1384, 2000.
  * [24] M. Lukáčová-Medvidová, S. Noelle, and M. Kraft.  Well-balanced finite volume evolution Galerkin methods for the shallow water equations.  Journal of computational physics, 221(1):122–147, 2007.
  * [25] M. Lukáčová-Medvid’ová, J. Saibertová, and G. Warnecke.  Finite volume evolution Galerkin methods for nonlinear hyperbolic systems.  Journal of Computational Physics, 183(2):533–562, 2002.
  * [26] M. Lukácová-Medvidová and E. Tadmor.  On the entropy stability of the Roe-type finite volume methods.  In Proceedings of the twelfth international conference on hyperbolic problems, American Mathematical Society, eds. J.-G. Liu et al, volume 67, 2009.
  * [27] M. Lukáčová-Medviďová, G. Warnecke, and Y. Zahaykah.  On the stability of evolution Galerkin schemes applied to a two-dimensional wave equation system.  SIAM journal on numerical analysis, 44(4):1556–1583, 2006.
  * [28] M. Lukáčová-Medvid’ová, G. Warnecke, and Y. Zahaykah.  Finite volume evolution Galerkin (FVEG) methods for three-dimensional wave equation system.  Applied numerical mathematics, 57(9):1050–1064, 2007.
  * [29] P. Monk and G. R. Richter.  A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media.  Journal of Scientific Computing, 22(1):443–477, 2005.
  * [30] S. Ostkamp.  Multidimensional characteristic Galerkin methods for hyperbolic systems.  Mathematical methods in the applied sciences, 20(13):1111–1125, 1997.
  * [31] J.-M. Qiu and C.-W. Shu.  Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov–Poisson system.  Journal of Computational Physics, 230(23):8386–8409, 2011.
  * [32] B. Riviere and M. F. Wheeler.  Discontinuous finite element methods for acoustic and elastic wave problems.  Contemporary Mathematics, 329(271-282):4–6, 2003.
  * [33] T. Russell and M. Celia.  An overview of research on Eulerian-Lagrangian localized adjoint methods (ELLAM).  Advances in water resources, 25(8):1215–1231, 2002.
  * [34] H. Wang, R. Ewing, G. Qin, S. Lyons, M. Al-Lawatia, and S. Man.  A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations.  Journal of Computational Physics, 152(1):120–163, 1999.
  * [35] K. Wang, H. Wang, and M. Al-Lawatia.  An Eulerian-Lagrangian discontinuous Galerkin method for transient advection-diffusion equations.  Numerical Methods for Partial Differential Equations: An International Journal, 23(6):1343–1367, 2007.
  * [36] Y. Xing, C.-S. Chou, and C.-W. Shu.  Energy conserving local discontinuous Galerkin methods for wave propagation problems.  Inverse Problems & Imaging, 7(3):967, 2013.