# Stabilizing Discontinuous Galerkin Methods Using Dafermos’ Entropy Rate
Criterion: II – Systems of Conservation Laws and Entropy Inequality Predictors

Simon-Christian<EMAIL_ADDRESS>
###### Abstract

A novel approach for the stabilization of the Discontinuous Galerkin method
based on the Dafermos entropy rate crition is presented. First, estimates for
the maximal possible entropy dissipation rate of a weak solution are derived.
Second, families of conservative Hilbert-Schmidt operators are identified to
dissipate entropy. Steering these operators using the bounds on the entropy
dissipation results in high-order accurate shock-capturing DG schemes for the
Euler equations, satisfying the entropy rate criterion and an entropy
inequality.

##  1 Introduction

Discontinuous Galerkin methods [5] are a popular tool to design numerical
schemes for hyperbolic systems of conservation laws [8]

${\frac{\partial{f(u)}}{\partial{x}}}+{\frac{\partial{u}}{\partial{t}}}=0\quad\text{for}\quad
u(x,t):\mathbb{R}\times\mathbb{R}\to\mathbb{R}^{m},\quad
f:\mathbb{R}^{m}\to\mathbb{R}^{m}.$ (1)

A cell of the subdivision $\mathcal{T}$ of the domain $\Omega$ | $T$
---|---
The left and right boundaries of cell $T$ |  $T_{l}$, $T_{r}$
The set of test functions | $\mathcal{D}$
The space of ansatz functions for an cell $T$ | $V^{T}=\operatorname{span}\\{\varphi_{1},\varphi_{2},\dots,\varphi_{N}\\}$
Polynomials of degree of $p$ in cell $T$ | $V^{T,p}$
$\mathrm{L}^{2}$ projection of $u$ onto $V$ |  $\operatorname{\mathbb{P}}_{V}u$.
Interpolation of $u$ on $V$ w.r.t. the collocation points $(\xi_{k})_{k=1}^{N}$ | $\operatorname{\mathbb{I}}_{V}u$
Vector of nodal values in cell $T$ at time $t$ | $u^{T}(t)$
ansatz function in cell $T$ at position $x$ and time $t$ | $u^{T}(x,t)$
Inner product on cell $T$ | $\left\langle{u},{v}\right\rangle_{T}=\int_{T}\left(u(x)\middle)\cdot\middle(v(x)\right)\,\mathrm{d}x$
Surface inner product on cell $T$ | $\left[{v},{f}\right]_{T}=\int_{\partial T}\left(v\middle)\cdot\middle(f\right)\,\mathrm{d}O$
Gramian Matrices on cell $T$ | $M_{k,l}^{T}=\left\langle{\varphi_{k}},{\varphi_{l}}\right\rangle_{T},\,S_{k,l}^{T}=\left\langle{{\frac{\partial{\varphi_{k}}}{\partial{x}}}},{\varphi_{l}}\right\rangle_{T}$
The total entropy in cell $T$ | $E_{u,T}(t)$
The discrete total entropy in cell $T$ | $E^{T}(t)$
The inner product on $T$ discretised using $\omega_{k}$ | $\left\langle{u},{v}\right\rangle_{T,\omega}=\sum_{k}\omega_{k}u_{k}v_{k}$
Entropy variables in cell $T$ | ${\frac{\partial{U}}{\partial{u}}}(u^{T}(x,t))$
Vector of nodal values of the entropy variables in cell $T$ at $t$ | $\operatorname{{\frac{\partial{U}}{\partial{u}}}}^{T}(t)$
Interpolation of the entropy variables in cell $T$ on $V$ | $\operatorname{{\frac{\partial{U}}{\partial{u}}}}^{T}(x,t)$
The canonical inner product between $a,b\in\mathbb{R}^{n}$ |  $a\cdot b$ or $\left(a\middle)\cdot\middle(b\right)$
The inner product between $u$ and $v$ on cell $T$ | $\left\langle{u},{v}\right\rangle_{T}$
The $p$ norm of $u$ in cell $T$ | $\left\|u\right\|_{T,p}$
Exact solution to the initial condition $u(x,t_{0})$ after $t-t_{0}$ |  $H(u(\cdot,t_{0}),t-t_{0})$.
Mean value of subcell $k$ of $N$, $u^{T}$ as initial condition | $u^{T,N}_{k}$
The convex hull of a set $A$ | $\operatorname{ch}A$
Table 1: Notation used. As a general rule, quantities with only $t$ as an
argument are vectors of nodal values at a certain time. Values with $x$ and
$t$ in their argument list are functions that were evaluated at these values.
Objects with $T$ added as exponent are approximations of the quantity in the
cell $T$.

An intriguing feature of DG methods is their ability to transfer the
definition of a weak solution to a hyperbolic conservation law [15]

$\displaystyle\forall\varphi\in\mathcal{D}:\int_{0}^{\infty}\int_{\mathbb{R}}\left(u(x,t)\middle)\cdot\middle({\frac{\partial{\varphi(x,t)}}{\partial{t}}}\right)$
$\displaystyle+\left(f\circ
u(x,t)\middle)\cdot\middle({\frac{\partial{\varphi(x,t)}}{\partial{x}}}\right)\,\mathrm{d}x\,\mathrm{d}t$
(2)
$\displaystyle+\int_{\mathbb{R}}\left(u(x,0)\middle)\cdot\middle(\varphi(x,0)\right)\,\mathrm{d}x=0.$

to the semidiscrete level [2, 4, 3]. Using a method of lines approach this
leads to the set of equations

$\forall
T\in\mathcal{T},\varphi\in\mathcal{D}:\left\langle{{\frac{\partial{\varphi}}{\partial{x}}}},{f}\right\rangle_{T}-\left[{\varphi},{f}\right]_{T}-\left\langle{\varphi},{{\frac{\partial{u^{T}}}{\partial{t}}}}\right\rangle_{T}=0$

for every cell $T\in\mathcal{T}$ of a subdivision $\mathcal{T}$ of the domain
into cells. The solution $u(x,t)$ is approximated in every cell by
$u^{T}(x,t)\in V^{T}$ out of a finite dimensional space of ansatz functions
$V^{T}$. Using an approximation of the inner products as point evaluations
results in the matrix vector form

$M^{T}{\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}=S^{T}f(u^{T}(t))-\begin{pmatrix}\varphi^{T}_{1}(x_{r})f^{*}_{r}-\varphi^{T}_{1}(x_{l})f^{*}_{l}\\\
\vdots\\\
\varphi^{T}_{N}(x_{r})f^{*}_{r}-\varphi^{T}_{N}(x_{l})f^{*}_{r}\end{pmatrix}.$
(3)

Sadly, the constructed schemes lack robustness and stability in the high order
case and some stabilization measures and robustness enhancements are needed.
Popular are overintegration, flux-differencing, modal filtering, sub-cells and
(W)ENO recoveries [14, 29, 12, 3, 42, 27, 34]. In this publication the
procedure first presented in [21] will be refined, connections to some other
stabilization techniques will be shown, and the technique will be tested on a
catalog of problems for the Euler system of conservation laws. The method in
[21] is based on the entropy rate admissibility criterion [6, 11]. An entropy
[24] is a convex functional $U:\mathbb{R}^{m}\to\mathbb{R}$ satisfying

${\frac{\mathrm{d}f}{\mathrm{d}u}}{\frac{\mathrm{d}U}{\mathrm{d}u}}={\frac{\mathrm{d}F}{\mathrm{d}u}}$

in conjunction with a entropy flux function $F:\mathbb{R}^{m}\to\mathbb{R}$.
One can show that for this pair $(F,U)$ holds

${\frac{\partial{U(u(x,t))}}{\partial{t}}}+{\frac{\partial{F}}{\partial{x}}}\leq
0$ (4)

in the sense of distributions [24]. If the solution is smooth one can even
show

${\frac{\partial{U(u(x,t))}}{\partial{t}}}+{\frac{\partial{F}}{\partial{x}}}=0.$

The entropy rate criterion states that the total entropy

$E_{u}(t)=\int U(u(x,t))\,\mathrm{d}x$

of the selected weak solution $u$ should reduce faster than the entropy of any
other existing weak solution $\tilde{u}$

$\forall
t>0:\quad{\frac{\mathrm{d}E_{u}}{\mathrm{d}t}}\leq{\frac{\mathrm{d}E_{\tilde{u}}}{\mathrm{d}t}}.$

A numerical approximation of this total entropy can be defined as

$E_{u,T}(t)=\int_{T}U(u^{T}(x,t))\,\mathrm{d}x\approx\sum_{k}\omega^{T}_{k}U(u^{T}(x_{k},t)),\quad
E_{u}=\sum_{T\in\mathcal{T}}E_{u,T}$ (5)

via a (positive) quadrature rule $\omega^{T}_{k}$ on each cell
$T\in\mathcal{T}$. The numerical enforcement of the criterion with respect to
such a definition of the discrete entropy happened in [21] in three steps

  * •

Calculate the time derivative of the ansatz function
${\frac{\mathrm{d}\tilde{u}^{T}}{\mathrm{d}t}}$ using a DG scheme

  * •

Calculate an error prediction $\delta^{T}$ for
${\frac{\mathrm{d}\tilde{u}^{T}}{\mathrm{d}t}}$ on $T$, i.e.
$\left\|{\frac{\mathrm{d}\tilde{u}^{T}}{\mathrm{d}t}}-{\frac{\partial{u}}{\partial{t}}}\right\|_{T}\leq\delta^{T}$.

  * •

Correct the time derivative into the direction of the steepest entropy descent

${\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}={\frac{\mathrm{d}\tilde{u}^{T}}{\mathrm{d}t}}-\frac{\delta^{T}}{\left\|h^{T}\right\|_{T}}h^{T},$

where $h$ shall be the steepest descent direction that does not change the
average value in cell $T$.

While the approach above is successful for scalar conservation laws [21]
significant improvements can be made by introducing two refinements. The first
one concerns the usage of an error indicator to estimate the entropy
correction needed. We will instead show that it is possible to directly give
bounds on how dissipative a weak solution can be. This will eliminate the need
for the error indicator while allowing a faster convergence, because the
derived bounds converge to zero significantly faster in the smooth case. A
second refinement concerns the direction used for the entropy correction. DG
methods can make use of modal filtering to remove unwanted high frequency
modes from the solution [29]. These filters can be sometimes expressed as
viscosity, and we will devise correction directions that at the same time
dissipate entropy and filter the solution from unwanted oscillations and
thereby combine the dissipation and filtering.

Our schemes will therefore follow the slightly different general layout of

  * •

Calculate a time derivative for the ansatz function
${\frac{\mathrm{d}\tilde{u}^{T}}{\mathrm{d}t}}$

  * •

Estimate the highest possible entropy dissipation speed $\sigma^{T}$ in cell
$T$

  * •

Calculate the correction direction $\upsilon^{T}$

  * •

Calculate the size $\lambda^{T}$ of the correction needed to achieve that

${\frac{\partial{u^{T}}}{\partial{t}}}={\frac{\partial{\tilde{u}^{T}}}{\partial{t}}}+\lambda^{T}\upsilon^{T}$

satisfies

${\frac{\mathrm{d}E_{u,T}}{\mathrm{d}t}}=\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{{\frac{\partial{\tilde{u}^{T}}}{\partial{t}}}+\lambda^{T}\upsilon^{T}}\right\rangle_{T}\leq\sigma^{T}+F^{*}_{l}-F^{*}_{r}.$
(6)

the dissipation mandated by the estimate. Here $F^{*}$ shall be a numerical
entropy flux [35, 36, 37].

The procedure makes only use of the fact that in our cells there exist local
ansatz functions and is therefore also applicable to similar schemes like the
spectral volume (SV) method [40]. The only difference would lie in the
evaluation of a different scheme for the uncorrected derivative
${\frac{\mathrm{d}\tilde{u}}{\mathrm{d}t}}$. One complication is brought in by
the fact that entropy dissipation implies the non-smoothness of the solution,
as otherwise the entropy equality applies. Therefore, dissipation can’t happen
in cells in the continuous setting, as polynomials are smooth. Instead,
dissipation is a process taking place at the cell edges were our different
ansatz functions transition. As we are not correcting the numerical fluxes
used between cells dissipation will be centered in cells and not at cell
edges, and we will show in section 3.1 how to work around this problem.

##  2 Entropy inequality predictors

###  2.1 Bounds for entropy and entropy dissipation

Our main tool to approximate the most dissipative weak solution using a DG
method will be a bound on the derivative of the total entropy. We will derive
a lower bound for the entropy dissipation

$s^{\theta}_{u}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}\int_{\theta}{\frac{\partial{U}}{\partial{t}}}+{\frac{\partial{F}}{\partial{x}}}\,\mathrm{d}x\,\mathrm{d}t\leq
0.$

Here $\theta\subset\Omega$ shall be an arbitrary open subdomain of the
complete domain. This value has to be smaller than zero for a solution that is
admissible with respect to the classical entropy inequality (4). Further, we
are interested in the entropy dissipation speed

$\sigma^{\theta}_{u}(t)={\frac{\partial{s^{\theta}_{u}(t_{0},t)}}{\partial{t}}}=\int_{\theta}{\frac{\partial{U}}{\partial{t}}}+{\frac{\partial{F}}{\partial{x}}}\,\mathrm{d}x\leq
0.$

If this value is known one can estimate the total entropy’s $E_{u}(t)$
derivative as

$E_{u}(t)\geq\sum_{\theta\in\Theta}s^{\theta}_{u}(0,t),\quad{\frac{\mathrm{d}E_{u}}{\mathrm{d}t}}\geq\sum_{\theta\in\Theta}\sigma^{\theta}_{u},$

when $\Theta=\\{\theta_{1},\theta_{2},\theta_{3},\dots,\theta_{L}\\}$ is
overlapping $\Omega$ in the sense of

$\Omega\subset\bigcup_{\theta\in\Theta}\theta.$

To achieve our goal of estimating $s^{\theta}$ we will view the problem in the
setting of classical Finite-Volueme schemes [33] and go over to the limit
$\Delta x\to 0$. In [7] it was shown that for scalar conservation laws the
flux $f$ of the solution to the Riemann problem
$u_{\text{R}}(u_{l},u_{r};x,t)$ is given by

$f\left(\operatorname{argmin}_{u\in\operatorname{ch}(u_{l},u_{r})}\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}(u_{r})-{\frac{\mathrm{d}U}{\mathrm{d}u}}(u_{l})},{f(u)}\right\rangle\right),$

i.e. by entering the value of $u$ into $f$ that when entered into the flux
yields the fastest entropy dissipation. In [22] it was shown that some
approximate Riemann solvers, for example the local Lax-Friedrichs flux, can be
also interpreted as approximate solutions to such variational descriptions of
two-point fluxes. While the aforementioned results hold for semidiscrete
schemes the new results below are new and aim at three point first order
Finite-Difference/Finite-Volume schemes for systems of conservation laws. As
one assumes piecewise constant functions in those first order methods any
quadrature exact for constants will yield the same result in equation (5). As
we only look at discrete time values in this part of the publication we will
write $E^{n}_{u}=E_{u}(t_{n})$ for the discrete total entropy at time level
$n$.

######  Lemma 1.

Let a system of hyperbolic conservation laws in conservation form and a
strictly convex entropy pair $(U,F)$ be given that is approximated by a
Finite-Volume scheme with grid constant $\lambda=\frac{\Delta t}{\Delta x}$.
Then the original Lax-Friedrichs scheme has the fastest dissipation of the
total entropy

$E^{n}_{u}=\sum_{k}U(u^{n}_{k})\Delta x$

under all consistent and conservative three-point numerical schemes.

###### Proof.

Assume $f(u_{k},u_{k+1})$ is a consistent numerical two point flux minimizing
the total entropy with maximal rate and let $u_{l},u_{r}$ be arbitrary in the
domain of admissible values for the conserved variables. We apply a scheme
using this flux to a Riemann problem, i.e. the initial data

$u^{0}_{k}=\begin{cases}u_{l}&k\leq 0\\\ u_{r}&k>0\\\ \end{cases}$

to query the flux value $f(u_{l},u_{r})=f(u_{0},u_{1})=f_{\frac{1}{2}}$ by
analyzing the solution. As the flux is consistent it holds

$\forall
k<0:f\left(u_{k},u_{k+1}\right)=f_{k+\frac{1}{2}}=f(u_{l}),\quad\forall
k>0:f\left(u_{k},u_{k+1}\right)=f_{k+\frac{1}{2}}=f(u_{l}).$

The scheme

$u^{1}_{k}=u^{0}_{k}+\lambda\left(f_{k-\frac{1}{2}}-f_{k+\frac{1}{2}}\right)$

therefore implies that $u^{1}_{k}=u^{0}_{k}$ for all $k\not\in\\{0,1\\}$. The
total entropy

$E^{1}_{u}=E^{0}_{u}-\Delta
x\left(U\left(u^{0}_{0}\right)+U\left(u^{0}_{1}\right)\right)+\Delta
x\left(U\left(u^{1}_{0}\right)+U\left(u^{1}_{1}\right)\right)$

is minimized by $u^{1}_{0}=u^{1}_{1}$, as $U$ is strictly convex. Entering
this into the scheme’s definition with $u^{0}_{0}=u_{l}$ and $u^{0}_{1}=u_{r}$
implies

$u_{l}+\lambda(f(u_{l})-f(u_{l},u_{r}))=u_{r}+\lambda(f(u_{l},u_{r})-f(u_{r})).$

Rearranging for $f(u_{l},u_{r})$ shows

$f(u_{l},u_{r})=\frac{f(u_{l})+f(u_{r})}{2}+\frac{u_{l}-u_{r}}{2\lambda},$

and this is the classical Lax-Friedrichs flux and therefore uniquely
determined by demanding maximal entropy rate. ∎∎

This result shows that the classical LF scheme is the most direct realisation
of a scheme satisfying Dafermos’ entropy rate criterion and therefore
justifies the use of the LF scheme in [20] as the most dissipative scheme
possible for systems of conservation laws. Similar results are also known for
scalar conservation laws. Tadmor showed in [35, 36] that every monotonicity
preserving scheme satisfying classical numerical entropy inequalities for a
scalar conservation law has a viscosity coefficient less or equal to that of
the LF scheme, and higher or equal than the viscosity coefficient of Godunov’s
scheme. Our result can be seen as a generalization of the LF part of this
result to systems of conservation laws, as it states that the LF flux is the
most dissipative flux for a selected time-step size. Using the scheme above
one can derive estimates for the highest possible entropy dissipation in a
time-step and using finite-differencing of this result, approximations for the
lowest possible derivative of the total entropy with respect to time.

######  corollary 1.

The biggest possible entropy dissipation during a discrete time-step of a
Finite-Volume scheme with grid constant $\lambda=\frac{\Delta t}{\Delta x}$ is
given by the difference

$E^{n+1}-E^{n}=\sum_{k}U\left(\frac{u_{k-1}^{n}+u_{k+1}^{n}}{2}+\lambda\left(f\left(u_{k-1}^{n}\right)+f\left(u_{k+1}^{n}\right)\right)\right)\Delta
x-U\left(u_{k}^{n}\right)\Delta x$

and an approximation to the total entropy’s minimal derivative by

${\frac{\mathrm{d}E}{\mathrm{d}t}}\approx\sum_{k}\frac{U\left(\frac{u_{k-1}+u_{k+1}}{2}+\lambda(f(u_{k-1})+f(u_{k+1}))\right)-U(u_{k})}{\lambda}.$

The second estimate above degenerates for $\lambda\to 0$ as the difference in
entropy is in general finite between cells $u_{k-1},u_{k},u_{k+1}$. A second,
more refined, estimate is given by the following lemma based on the ideas from
[18] and does not have these deficiencies.

######  Lemma 2.

Given bounds on the fastest signal speed to the left $a_{l}$ and the highest
signal speed to the right $a_{r}$ let
$M\geq\max(\left|a_{l}\right|,\left|a_{r}\right|)$. The maximum entropy
dissipation of a Riemann problem solution on the interval $\theta=(-M,M)$ is
bounded from below by

$E^{\theta}_{u}(t)-E^{\theta}_{u}(0)\geq
t\left((a_{r}-a_{l})U\left(u_{lr}\right)+a_{l}U(u_{l})-a_{r}U(u_{r})\right)$

with

$u_{lr}=\frac{a_{r}u_{r}-a_{l}u_{l}+f(u_{l})-f(u_{r})}{a_{r}-a_{l}}.$

The rate is bounded from below by

$\begin{aligned}
{\frac{\mathrm{d}E^{\theta}}{\mathrm{d}t}}|_{t=0}\geq(a_{r}-a_{l})U(u_{lr})+a_{l}U(u_{l})-a_{r}U(u_{r})\end{aligned}.$

The entropy dissipation is bounded by

$s^{\theta}\geq
t\left((a_{r}-a_{l})U(u_{lr})+a_{l}U(u_{l})-a_{r}U(u_{r})+F(u_{l})-F(U_{r})\right),$

and its rate by

$\sigma^{\theta}\geq(a_{r}-a_{l})U(u_{lr})+a_{l}U(u_{l})-a_{r}U(u_{r})+F(u_{l})-F(U_{r}).$

###### Proof.

$0$$t$$x$$a_{l}t=x$$a_{r}t=x$$(a_{r},1)$$1.0$$(a_{l},1)$ (a) The case
$a_{l}<0<a_{r}$

$0$$t$$x$$a_{l}t=x$$a_{r}t=x$$(a_{r},1)$$(a_{l},1)$ (b) The case
$0<a_{l}<a_{r}$ and vice-versa.

Figure 1: Layout of the integration areas in the proof (blue). As originally
used in [18]. To the left and right of the lines
${\frac{\mathrm{d}x}{\mathrm{d}t}}=a_{l}$ and
${\frac{\mathrm{d}x}{\mathrm{d}t}}=a_{r}$ the initial condition is unaltered.

The entropy of the initial condition in the interval $[-M,M]$ is given by

$\int_{-M}^{M}U(u(x,0))\,\mathrm{d}x=MU(u_{l})+MU(u_{r})$

for any $M>0$. Integrating over the triangle
$T=\operatorname{ch}\\{(0,0),(a_{l},1),(a_{r},1)\\}$ in spacetime and using
the conservation law yields

$\displaystyle 0$
$\displaystyle=\int_{T}{\frac{\partial{u}}{\partial{t}}}+{\frac{\partial{f}}{\partial{x}}}\,\mathrm{d}V(x,t)=\int_{\partial
T}\begin{pmatrix}f\\\ u\end{pmatrix}\cdot n\,\mathrm{d}O(x,t)$
$\displaystyle=\int_{a_{l}}^{a_{r}}\begin{pmatrix}f(u)\\\
u(x,1)\end{pmatrix}\cdot\begin{pmatrix}0\\\
1\end{pmatrix}\,\mathrm{d}x+\int_{0}^{1}\begin{pmatrix}f(u)\\\
u(ta_{l},t)\end{pmatrix}\cdot\begin{pmatrix}-1\\\
a_{l}\end{pmatrix}\,\mathrm{d}t+\int_{0}^{1}\begin{pmatrix}f(u)\\\
u(ta_{r},t)\end{pmatrix}\cdot\begin{pmatrix}1\\\
-a_{r}\end{pmatrix}\,\mathrm{d}t$
$\displaystyle=(a_{r}-a_{l})u_{lr}+a_{l}u_{l}-a_{r}u_{r}+f(u_{r})-f(u_{l}),$

in conjunction with the Gauß divergence theorem, cf. figure 1. Here $u_{lr}$
shall denote the mean value of $u(x,1)$ on $[a_{l},a_{r}]$ and is

$u_{lr}=\frac{a_{r}a_{l}-a_{l}u_{l}+f(u_{l})-f(u_{r})}{a_{r}-a_{l}}$

as apparent from the calculation above. Jensens inequality implies

$\displaystyle t(a_{r}-a_{l})U(u_{lr})$
$\displaystyle=t(a_{r}-a_{l})U\left(\frac{1}{t(a_{r}-a_{l})}\int_{ta_{l}}^{ta_{r}}u(x,t)\,\mathrm{d}x\right)$
(7)
$\displaystyle\leq\frac{t(a_{r}-a_{l})}{t(a_{r}-a_{l})}\int_{ta_{l}}^{ta_{r}}U(u(x,t))=E^{(ta_{l},ta_{r})}_{u}(t).$

Therefore it follows

$\displaystyle E^{\theta}(1)-E^{\theta}(0)\geq$
$\displaystyle(a_{r}-a_{l})U(u_{lr})+(M-a_{r})U(u_{r})+(M+a_{l})U(u_{l})$ (8)
$\displaystyle-M(U(u_{l})+U(u_{r}))$ $\displaystyle=$
$\displaystyle(a_{r}-a_{l})U(u_{lr})+a_{l}U(u_{l})-a_{r}U(u_{r})$

for the entropy dissipation between $t=0$ and $t=1$ and using the invariance
under transformations $(x,t)\mapsto(\mu x,\mu t)$ for $\mu>0$ yields

$\displaystyle{\frac{\mathrm{d}E}{\mathrm{d}t}}|_{t=0}\geq(a_{r}-a_{l})U(u_{lr})+a_{l}U(u_{l})-a_{r}U(u_{r})$
(9)

for the rate. To calculate the entropy dissipation $s^{\theta}$ and its speed
$\sigma^{\theta}$ we just have to account for the entropy flowing in and out
of the intervall $\theta$ using the entropy flux $F$. This is possible as $u$
is constant to the left of $(ta_{l},t)$ and to the right of $(ta_{r},t)$. ∎∎

The estimate above does not depend on any grid constant, and reduces to the
previous one for $-a_{l}=c_{\text{max}}=a_{r}$, $\lambda c_{\text{max}}=1$,
and this is the CFL condition for the classical Lax-Friedrich scheme, i.e.
both estimates are compatible. A Godunov type scheme using the HLL approximate
Riemann solver is also compatible with the estimate above. The discrete total
entropy after one time-step is still less or equal than the bound given above.

$x$$t$$x$$u(x,t_{n+1})$ Figure 2: A set of noninteracting HLL approximate
Riemann solutions.

Let $\lambda c_{\text{max}}\leq\frac{1}{2}$ hold implying that the Riemann
problems do not interact and $u^{\mathrm{HLL}}(x,t_{n+1})$ be the picewise
constant solution of the HLL solver as in figure 2, but not averaged over the
cells, while $u^{n+1}_{k}$ shall be the corresponding cell averages. In this
case the total discrete entropy at the next time-step is given by

$\displaystyle E^{n+1}_{\mathrm{FV}}=$ $\displaystyle\sum_{k}\Delta
xU(u_{k}^{n+1})\leq\sum_{k}\frac{\Delta x}{\Delta
x}\int_{x_{k-\frac{1}{2}}}^{x_{k+\frac{1}{2}}}U(u^{HLL}(x,t))\,\mathrm{d}x$
$\displaystyle=$ $\displaystyle\int_{\Omega}U(u^{HLL}(x,t))\,\mathrm{d}x\leq
E^{n+1}_{u^{\mathrm{HLL}}}.$

Therefore the discrete entropy of the approximate solution is lower than the
entropy of any exact weak solution. The next subsection will move beyond first
order schemes by generalizing this lower bound to one that also allows smooth
solutions instead of piecewise constant ones.

###  2.2 Asymptotic analysis based entropy inequality predictor

$x$$u(x,0)$ (a) The problem for the generalized entropy inequality predictor.
Two piecewise smooth solutions are spliced together and we are interested in
the local residual of the entropy equality around the interface.

$x$$t$$a_{l}$$a_{r}$ (b) Solutions of generalized Riemann problems lack
scaling invariance. Still we assume that there exist bounds $a_{l}$ and
$a_{r}$ that the waves from the interaction of $u_{l}(x)$ and $u_{r}(x)$ do
not leave the cone $[ta_{l},ta_{r}]$ for small t.

$x$$u(x,t)$$ta_{l}$$ta_{r}$ (c) The assumed solution used in the generalized
entropy inequality predictor. Note that this solution follows the HLL idea of
assuming a constant function in the wedge formed by $ta_{l}$ and $ta_{r}$.

$x$$u(x,t)$$x_{k-\frac{5}{2}}$$x_{k-\frac{3}{2}}$$x_{k-\frac{1}{2}}$$x_{k+\frac{1}{2}}$$x_{k+\frac{3}{2}}$$x_{k+\frac{5}{2}}$$x_{k-2}$$x_{k-1}$$x_{k}$$x_{k-1}$$x_{k-2}$$\theta_{k-\frac{3}{2}}$$\theta_{k-\frac{1}{2}}$$\theta_{k+\frac{1}{2}}$$\theta_{k+\frac{3}{2}}$
(d) Application of the entropy inequality predictor to an open overlap
$\theta_{k+\frac{1}{2}}$ of the domain - centered on the discontinuities.

Figure 3: Construction and application of the generalized HLL dissipation
estimate.

The entropy inequality predictor in this section will be based on an
asymptotic analysis of the problem described in figure 3(a), i.e. two smooth
solutions splined together at an interface. An obstacle lies in the missing
self-similarity. This is a difference to the previous part where the self-
similarity of the initial condition and assumed self-similarity of the
solution induced the existence of a self-similar, i.e. constant, speed of the
entropy dissipation. We will therefore try to approximate

$s^{\theta}(t_{1},t_{2})=\int_{t_{1}}^{t_{2}}\int_{\theta}{\frac{\partial{F}}{\partial{x}}}+{\frac{\partial{U}}{\partial{t}}}\,\mathrm{d}x\,\mathrm{d}t$

for reasonably small $\left|t_{2}-t_{1}\right|$ and the discontinuity at the
interface in the interior of $\theta$. The schemes in which we will use these
entropy inequality predictors should converge with high orders for smooth
solutions, necessitating a convergence of the predictor to zero with a high
order for smooth solutions. This convergence is also dictated by the entropy
equality for smooth solutions. If $u_{l}(x)$ and $u_{r}(x)$ are piecewise
constant this problem is already solved by the methods described in the last
subsection. We will therefore now reiterate through the proof of lemma 2
assuming that $u_{l}(x)$ and $u_{r}(x)$ are smooth functions. The missing
self-similarity of Generalized Riemann problems [1], cf. figure 3(b), defies
the existence of the speed estimates $a_{l}$ and $a_{r}$, and we therefore
just assume that these speed estimates exist for small times. Further we
assume that for small times the solutions left of $(ta_{l},t)$ and right of
$(ta_{r},t)$ remain smooth, as no waves from the interaction arrive there and
$u_{l}(x),u_{r}(x)$ have bounded derivatives.

The average value $u_{lr}$ shall be determined by applying the conservation
law to the triangle $T=\operatorname{ch}\\{(0,0),(ta_{l},t),(ta_{r},t)\\}$

$\displaystyle 0=$
$\displaystyle\int_{T}{\frac{\partial{u}}{\partial{t}}}+{\frac{\partial{f}}{\partial{x}}}\,\mathrm{d}V(x,t)=\int_{\partial
T}\begin{pmatrix}f\\\ u\end{pmatrix}\cdot n\,\mathrm{d}O(x,t)$
$\displaystyle=$ $\displaystyle\int_{ta_{l}}^{ta_{r}}\begin{pmatrix}f(u)\\\
u(x,t)\end{pmatrix}\cdot\begin{pmatrix}0\\\
1\end{pmatrix}\,\mathrm{d}x+\int_{0}^{t}\begin{pmatrix}f(u)\\\ u(\tau
a_{l},\tau)\end{pmatrix}\cdot\begin{pmatrix}-1\\\
a_{l}\end{pmatrix}\,\mathrm{d}\tau$
$\displaystyle+\int_{0}^{t}\begin{pmatrix}f(u)\\\ u(\tau
a_{r},\tau)\end{pmatrix}\cdot\begin{pmatrix}1\\\
-a_{r}\end{pmatrix}\,\mathrm{d}\tau$ $\displaystyle=$
$\displaystyle\underbrace{\int_{ta_{l}}^{ta_{r}}u(x,t)\,\mathrm{d}x}_{t(a_{r}-a_{l})u_{lr}}+\int_{0}^{t}f(u_{r}(\tau
a_{r},\tau))-f(u_{l}(\tau a_{l},\tau))\,\mathrm{d}\tau$
$\displaystyle-\int_{ta_{l}}^{0}U(u_{l}(x,0))\,\mathrm{d}x-\int_{0}^{ta_{r}}U(u_{r}(x,0))\,\mathrm{d}x.$

Dividing this equation by $t$ and going over to the limit $t\to 0$ results in

$\displaystyle\frac{\int_{0}^{t}f(u_{r}(\tau a_{r},\tau))-f(u_{l}(\tau
a_{l},\tau))\,\mathrm{d}\tau}{t}$ $\displaystyle\xrightarrow{t\to
0}f(u_{r}(0,0))-f(u_{l}(0,0)),$
$\displaystyle\frac{\int_{ta_{l}}^{0}U(u_{l}(x,0))\,\mathrm{d}x-\int_{0}^{ta_{r}}U(u_{r}(x,0))\,\mathrm{d}x}{t}$
$\displaystyle\xrightarrow{t\to 0}a_{l}U(u_{l}(0,0))-a_{r}U(u_{r}(0,0))$

using the continuity of the integrands and the mean value theorem of
integration [26]. Therefore it follows

$u_{lr}(t)\xrightarrow{t\to
0}\frac{a_{r}u_{r}(0)-a_{l}u_{l}(0)+f(u_{l}(0))-f(u_{r}(0))}{a_{r}-a_{l}}$

for vanishing $t$. Equation (7) stays also valid in the case of piecewise
polynomial functions as initial conditions and for small $t>0$. We can
therefore conclude that a generalization of equation (8) holds in the form

$\displaystyle E^{\theta}(t)-E^{\theta}(0)\geq$ $\displaystyle
t(a_{r}-a_{l})U(u_{lr})+\int_{-M}^{ta_{l}}U(u(x,t))\,\mathrm{d}x+\int_{ta_{r}}^{M}U(u(x,t))\,\mathrm{d}x$
$\displaystyle-\int_{-M}^{M}U(u(x,0))\,\mathrm{d}x.$

Accounting for the entropy flowing in and out of $[-M,M]$ yields

$\displaystyle s^{\theta}(0,t)\geq$ $\displaystyle
t(a_{r}-a_{l})U(u_{lr})+\int_{-M}^{ta_{l}}U(u(x,t))\,\mathrm{d}x+\int_{ta_{r}}^{M}U(u(x,t))\,\mathrm{d}x$
$\displaystyle-\int_{-M}^{M}U(u(x,t))\,\mathrm{d}x+\int_{0}^{t}F(u(-M,\tau))-F(u(M,\tau))\,\mathrm{d}\tau.$

Applying the entropy equality to the subdomains $[-M,ta_{l}]\times[0,t]$ and
$[ta_{r},M]\times[0,t]$

$\int_{-M}^{ta_{l}}U(u(x,t))\,\mathrm{d}x-\int_{-M}^{ta_{l}}U(u(x,0))\,\mathrm{d}x=\int_{0}^{t}F(u(-M,\tau))-F(u(ta_{l},\tau))\,\mathrm{d}\tau,$

that holds for small $t>0$ because the solution stays smooth in the
subdomains, allows us to restate this as

$\displaystyle s^{\theta}(0,t)\geq
t(a_{r}-a_{l})U(u_{lr})-\int_{ta_{l}}^{ta_{r}}U(u(x,t))\,\mathrm{d}x-\int_{0}^{t}F(u(ta_{l},\tau))-F(u(ta_{r},\tau))\,\mathrm{d}\tau.$

Dividing by $t$ and going over to the limit, using the limit of $u_{lr}$ and
once more the mean value theorem, shows in this case also

$\sigma^{\theta}\geq(a_{r}-a_{l})U(u_{lr})-a_{r}U(u_{r}(0))+a_{l}U(u_{l}(0))-F(u_{l}(0))-F(u_{r}(0)).$
(10)

A significant problem of the derivation above lies in the fact that one can
only estimate the entropy dissipation speed in the interval $\theta=(-M,M)$,
but not in $(-M,0)$ as the true dissipation can be located anywhere in the
cone $[ta_{l},ta_{r}]$. As the cells in our numerical tests will be layed out
as in figure 3(d)

$\mathcal{T}=\left\\{T_{k}=\left[x_{k-\frac{1}{2}},x_{k+\frac{1}{2}}\right]\leavevmode\nobreak\
\middle|\leavevmode\nobreak\ k\in\mathbb{Z}\right\\},\quad
x_{k-\frac{1}{2}}<x_{k+\frac{1}{2}}$

is a suitable set of overlapping open intervals

$\Theta=\left\\{\theta_{k+\frac{1}{2}}=(x_{k}-\varepsilon,x_{k+1}+\varepsilon)\leavevmode\nobreak\
\middle|\leavevmode\nobreak\ k\in\mathbb{Z}\right\\},\quad
x_{k}=\frac{x_{k-\frac{1}{2}}+x_{k+\frac{1}{2}}}{2}.$

We are therefore left with the problem of how to split this dissipation onto
the two neighboring cells that have overlap with $\theta_{k+\frac{1}{2}}$.
This problem will be handled below in section 3.1.

###  2.3 Accounting for aliasing errors

In [21] one of the findings in the numerical tests section was that the
entropy dissipation of the numerical solutions started already shortly prior a
real entropy dissipating discontinuity formed. This was attributed to the fact
that while the entropy of the exact solution is still constant as long as the
solution is smooth this exact solution will in general not be representable in
our ansatz space. It is therefore wise to dissipate entropy to arrive at a
function that still lies in our space, and certainly better than selecting an
ansatz function that has more entropy than the true solution. A similar issue
could be the fact that in the $\mathrm{L}^{p}$ norms, for $p<\infty$, near
each piecewise continuous solution $u^{T}$ lies a $\mathcal{C}^{\infty}$
function that can be constructed via mollification. Therefore an infinitely
small perturbation of $u^{T}$ in the usual norms leads to a vanishing entropy
dissipation. Or, put differently, the dissipation bound as a functional is
discontinuous in the $L^{p}$ spaces. While unsatisfactory let us remark that
the functional is better behaved with respect to the $BV$ semi norms. The
discontinuity of the entropy dissipation bound is problematic with under-
resolved solutions where a lucky, or in this case better to be considered
unlucky, too smooth approximation of the solution in our piecewise polynomial
spaces induces wrong, i.e. too conservative entropy dissipation predictions.

We are therefore interested in allowing our entropy inequality predictor to be
also greedy, or one could say pessimistic, with respect to an under-resolved
solution. The key to this strengthening is the following lemma.

######  Lemma 3 (Order of the entropy dissipation bound).

The maximal entropy dissipation prediction (10) of a Riemann problem for a
smooth flux function with smooth entropy-entropy flux pair vanishes
quadratically with the jump of $u$ at the interface

$\left|\sigma^{\theta}\right|\in\operatorname{\mathcal{O}}(\left\|u_{l}-u_{r}\right\|^{2}).$

###### Proof.

As the entropy inequality holds it is clear that the entropy dissipation is
non-positive in the sense of distributions. As we only allow entropy
dissipative solutions the entropy dissipation on $\theta$ is a non-positive
constant for a fixed jump. On the contrary (10) has to be zero for
$u_{l}=u_{r}$ and is smooth, implying that the line $u_{l}=u_{r}$ consists of
local maxima. Therefore, a power expansion of (10) around $u_{l}$ in $u_{r}$
has to be of the form

$\sigma^{\theta}(u_{l},u_{r})=\left(u_{l}-u_{r}\middle)\cdot\middle(H(u_{l}-u_{r})\right)+\operatorname{\mathcal{O}}{\left\|u_{l}-u_{r}\right\|^{3}}$

with a negative semi-definite Hessian $H\in\mathbb{R}^{m\times m}$. This
proofs the claim. ∎∎

We are therefore in the relaxing position that even if our approximations of
$u_{l},u_{r}$ only satisfy
$\left\|u_{l}-u_{r}\right\|\in\operatorname{\mathcal{O}}((\Delta x)^{p})$ the
corresponding estimate will converge significantly faster with
$\sigma^{\theta}(u_{l},u_{r})\in\operatorname{\mathcal{O}}((\Delta x)^{2p})$.
Our basic DG method predicts values for our solution $u^{T}$ in a Hilbertspace
that is spanned by polynomials on every cell. In this case a suitable
orthonormal basis is spanned by Legendre polynomials and the limits of these
basis representations are $\mathrm{L}^{2}$ functions. But as explained before
our functional is not continuous on $\mathrm{L}^{2}$ and our ansatz $u^{T}$ is
only an approximation of a projection of the true solution onto our ansatz
space. We can therefore try to exploit different projections of our ansatz,
especially projections that assume less regularity of $u^{T}$, and estimate
our entropy dissipation with the strongest one encountered in all of these
different approximations of $u_{l}$ and $u_{r}$. A natural choice for
projections on spaces assuming less regularity are projections on lower order
polynomials. As the Legendre polynomials on each cell when truncated up to
polynomial $p$ are an orthogonal basis of the polynomials with degree less
than or equal to $p$, is the orthogonal projection onto these spaces given by
discarding the higher order coefficients in the Legendre expansion of $u^{T}$.
We can truncate down to order $p-1$ by discarding the highest coefficient and
still achieve a convergence order of at least $q=2p-2>p$ of our entropy
inequality predictor for $p>2$. This can be summed up in the following
procedure used above order $p=2$.

  * •

Assign
$\sigma^{\theta_{k+\frac{1}{2}}}_{p}=\sigma^{\theta_{k+\frac{1}{2}}}_{u(x,t)}$.

  * •

Project the ansatz $u^{T}$ in every cell onto $V^{T,p-1}$ using an orthonormal
projection

$u^{T,p-1}=\operatorname{\mathbb{P}}_{V^{T,p-1}}u(\cdot,t).$

  * •

Assign
$\sigma^{\theta_{k+\frac{1}{2}}}_{p-1}=\sigma^{\theta_{k+\frac{1}{2}}}_{u^{p-1}(\cdot,t)}$.

  * •

Use
$\min\left(\sigma^{\theta_{k+\frac{1}{2}}}_{p},\sigma^{\theta_{k+\frac{1}{2}}}_{p-1}\right)$
as entropy inequality prediction.

##  3 Suitable dissipation directions and filtering

$u^{T}$$-{\frac{\mathrm{d}U}{\mathrm{d}u}}$$E_{u}=\text{const.}$$-{\frac{\partial{f}}{\partial{x}}}$${\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}$$\nabla^{2}u^{T}$$\nabla^{8}u^{T}$
Figure 4: Use of alternative dissipation directions. The direction
$-{\frac{\partial{f}}{\partial{x}}}$ shall denote the $\mathrm{L}^{2}$
projection of the exact solutions’ derivative
${\frac{\partial{u}}{\partial{t}}}$ onto $V$. Our approximation
${\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}$ is not entropy dissipative in this
example and should be corrected into the entropy dissipative half-space
characterized by the normal $-{\frac{\mathrm{d}U}{\mathrm{d}u}}$. The
direction $\nabla^{2}u$ that stems from a discretization of the heat kernel is
suitable for this correction and has additional benefits compared to
$-{\frac{\mathrm{d}U}{\mathrm{d}u}}$. While the diffusion also has a smoothing
effect the addition of $-{\frac{\mathrm{d}U}{\mathrm{d}u}}$ can even result in
a sharpening effect. Higher even derivatives like $\nabla^{8}u$ will smooth
the solution but will not result in a dissipation for all entropies.

After deriving approximations for the entropy dissipation needed we will now
determine how to correct the time derivative of the DG scheme to dissipate the
amount of entropy needed. At the same time the resulting scheme is hopefully
still high order accurate for entropy conservative solutions. In the scalar
case the direction of the steepest descent of the entropy, corrected for
conservation, was used for this purpose. This approach incurs several
problems:

  * •

The direction of steepest entropy descent has in general no
smoothing/filtering effect.

  * •

Proving in the previous publication that a correction in the steepest descent
direction with the length taken from the error indicator results in enough
entropy dissipation was possible but resulted in highly technical
arguments[21].

  * •

Dissipation stems from the viscous and parabolic history of hyperbolic
conservation laws. A viscous flux

$f_{\varepsilon}\left(u,{\frac{\partial{u}}{\partial{x}}}\right)=f(u)-\varepsilon{\frac{\partial{u}}{\partial{x}}}$

associated with a viscous regularization of a hyperbolic conservation law is
proportional to the gradient of the solution for fixed viscosity and
proportional to the gradient of the solution for constant viscosity. If a
component of $u$ is smooth with a low magnitude of the first and second
derivative the viscous flux of this component will also only differ from the
hyperbolic flux by a small margin. If our scheme is corrected with the
steepest entropy descent direction one can ask if this correction can be
expressed using some viscosity distribution $\varepsilon(x)$ in the domain.
This will be false in general. Even worse, the steepest gradient descend of
the entropy can’t be bounded using the first derivative of the respective
component of the vector valued function $u(x,t)$, incurring an infinitely
large viscosity.

All of the above reasons motivate us to devise alternative directions for the
entropy correction. These alternative directions should have the following
properties

  * •

The dissipation direction should have a filtering effect, i.e. when the
direction only is used high order modes should be dissipated.

  * •

The direction should dissipate entropy.

  * •

The dissipation should stem from a viscosity added to the hyperbolic flux.

Our new correction directions will be based on the construction of filters,
i.e. operators that can regularize a solution $u$. A filter will in our case
be a special Hilbert-Schmidt operator $K$ [25].

######  Definition 1 (Filter).

An operator $K:\mathrm{L}^{2}(\Omega)\to\mathrm{L}^{2}(\Omega)$ is said to be
a filter if it is an integral operator whose pointwise evaluation results in a
weighted average, i.e.

$[Ku](x)=\int_{\Omega}k(x,y)u(y)\,\mathrm{d}y,\quad\text{with }\forall
x\in\Omega:\int_{\Omega}k(x,y)\,\mathrm{d}y=1.$

is satisfied and the kernel $k$ is of bounded Hilbert-Schmidt norm.

We are especially interested in conservative filters as they do not destroy
the conservation of our basic schemes when they are applied on a per cell
basis.

######  Lemma 4 (Conservative filter).

A filter $K:\mathrm{L}^{2}(\Omega)\mapsto\mathrm{L}^{2}(\Omega)$ is
conservative

$\int_{\Omega}[Ku](x)\,\mathrm{d}x=\int_{\Omega}u(x)\,\mathrm{d}x$

if it can be written as an integral operator with a kernel with mass one, i.e.

$[Ku](x)=\int_{\Omega}u(y)k(x,y)\,\mathrm{d}y,\quad\text{with }\forall
y\in\Omega:\int_{\Omega}k(x,y)\,\mathrm{d}x=1$

###### Proof.

Using Fubini’s theorem shows

$\int_{\Omega}[Ku](x)\,\mathrm{d}x=\int_{\Omega}\int_{\Omega}k(x,y)u(y)\,\mathrm{d}y\,\mathrm{d}x=\int_{\Omega}\underbrace{\int_{\Omega}k(x,y)\,\mathrm{d}x}_{=1}\,u(y)\,\mathrm{d}y=\int_{\Omega}u(y)\,\mathrm{d}y$

in this case. ∎∎

Please note that the weighted average property is stated using the integration
w.r.t. the second variable while the conservation results from the unit
measure in the first variable. Obviously a convolution with a convolution
kernel satisfying

$\int_{\mathbb{R}}k(y)\,\mathrm{d}y=1$

satisfies both as $k(x,y)=k(x-y)$ holds in this case, but not every operator
satisfying these properties is a convolution. Especially when one is
interested in bounded domains convolutions are not an option, but there still
exist suitable smoothing operators.

######  Theorem 1 (Universally dissipative filters).

A conservative filter $K$ is dissipative for all convex entropies $U$,

$E_{Ku}=\int_{\Omega}U([Ku](x))\,\mathrm{d}x\leq\int_{\Omega}U(u(x))\,\mathrm{d}x=E_{u},$

if it can be written as a conservative filter with a positive kernel, i.e.

$[Ku](x)=\int_{\Omega}k(x,y)u(y)\,\mathrm{d}y$

with $\forall x,y:k(x,y)\geq 0$.

###### Proof.

Using Jensens inequality [30], the positivity and conservation of the filter
allows us to show

$\begin{aligned}
\int_{\Omega}U([Ku](x))\,\mathrm{d}x=&\int_{\Omega}U\left(\int_{\Omega}k(x,y)u(y)\,\mathrm{d}y\right)\,\mathrm{d}x\leq\int_{\Omega}\int_{\Omega}k(x,y)U(u(y))\,\mathrm{d}y\,\mathrm{d}x\\\
=&\int_{\Omega}\underbrace{\int_{\Omega}k(x,y)\,\mathrm{d}x}_{=1}U(u(y))\,\mathrm{d}y=\int_{\Omega}U(u(y))\,\mathrm{d}y\end{aligned}.$

∎∎

These theorem shows that the first and second bullet above can be satisfied by
an integral operator with a suitable kernel. An example of a dissipation that
can be identified with a positive conservative filter is the filtering by the
time evolution of

${\frac{\partial{u}}{\partial{t}}}=\varepsilon\nabla^{2}_{x}u$

on the entire domain as the assorted filter has the heat kernel as kernel
function [10],

$k^{t}(x,y)=h(x-y,t),\quad
h(x,t)=\frac{\mathrm{e}^{-\frac{\left|x\right|^{2}}{4t}}}{(4\pi
t)^{n/2}},\quad t>0.$

Further, this filtering obviously stems from viscosity and has therefore a
direct physical interpretation. It is known that while a positive integral
operator always dissipates entropy a high order finite-difference
implementation will not dissipate all entropies [28] and similar theorems hold
for higher even derivatives even in the analytic case. We will therefore
outline how to construct a filter that is dissipative in the semidiscrete and
fully discrete setting and can therefore be used as a descent direction. We
begin by stating some discrete equivalents of the theorems above and will
analyze if usual dissipations/filters satisfy this property. We will assume
that $\omega_{k}\geq 0$ is a positive quadrature rule on the cell $T$ for the
rest of the chapter and all notions of conservation for our filters will be
centered around being conservative with respect to this quadrature rule. For a
general DG method with dense mass matrix a quadrature can be calculated via
$\sum_{l}M_{lk}=\omega_{k}$, i.e. by entering the constant one into the
discretised inner product, but positivity is not guaranteed in general. A
general view of our plan could be to not discretise the second derivative, but
its action as the generator of a Hilbert-Schmidt operator. We will therefore,
when given a discrete filter, consider also its (discrete) generator.

######  Definition 2 (Conservative and positive filter generator).

Let $G\in\mathbb{R}^{(p+1)\times(p+1)}$ be a square matrix. We call this
matrix a filter generator if

$\forall k\in\\{1,\dots,p+1\\}:\quad\sum_{l=1}^{p+1}G_{kl}=0$

holds. It will be conservative if

$\forall l\in\\{1,\dots,p+1\\}:\quad\sum_{k=1}^{p+1}\omega_{k}G_{kl}=0$

is satisfied. Further, we call it positive, if

$\forall l\in\\{1,\dots,p+1\\},\quad\forall
k\in\\{1,\dots,l-1,l+1,\dots,p+1\\}:\quad G_{kl}\geq 0$

holds.

######  Definition 3 (Discrete conservative and positive filter).

We call a matrix $\Upsilon\in\mathbb{R}^{(p+1)\times(p+1)}$ a filter, if

$\forall k\in\\{1,\dots,p+1\\}:\quad\sum_{l=1}^{p+1}\Upsilon_{kl}=1$

holds. It is termed conservative, if

$\forall
l\in\\{1,\dots,p+1\\}:\quad\sum_{k=1}^{p+1}\omega_{k}\Upsilon_{kl}=\omega_{l}$

is satisfied. Further, we call it positive, if

$\forall k,l\in\\{1,\dots,p+1\\}:\quad\Upsilon_{kl}\geq 0.$

Obviously, the definition of the conservative positive discrete filter mirrors
the definition of such a filter in the continuous case using the quadrature
rule. The definition of the averaging property on the other hand is not based
on the quadrature rule, as this rule is not used when applying the filter
pointwise

$\upsilon_{k}=\sum_{l=1}^{p+1}\Upsilon_{kl}u_{l}.$

Forward Euler steps connect the generators defined above with the filters, as
we will see in the lemma below.

######  Lemma 5 (Connecting generators and filters).

It holds

$G\text{ conservative as generator
}\implies\Upsilon=\operatorname{\mathrm{I}}+\Delta tG\text{ conservative as
filter.}$

Let further $\Delta t\max_{l}\left|G_{ll}\right|\leq 1$. Then it follows

$G\text{ positive as generator}\implies\Upsilon\text{ positive as filter}.$

###### Proof.

We begin by showing the conservativity and filter property. It holds

$\sum_{l=1}^{p+1}\operatorname{\mathrm{I}}_{kl}=1\implies\sum_{l=1}^{p+1}\Upsilon_{kl}=\sum_{l=1}^{p+1}(\operatorname{\mathrm{I}}+\Delta
tG)_{kl}=1.$

As the identity is conservative follows

$\sum_{k=1}^{p+1}\omega_{k}\operatorname{\mathrm{I}}_{kl}=\omega_{l}\implies\sum_{k=1}^{p+1}\omega_{k}\Upsilon_{kl}=\sum_{k=1}^{p+1}\omega_{k}(\operatorname{\mathrm{I}}+\Delta
tG)_{kl}=\omega_{l}.$

The positivity follows as for non-diagonal elements,

$k\neq l\implies(I+\Delta tG)_{kl}=\Delta tG_{kl}\geq 0$

is satisfied for any positive timestep size while the given restriction is
needed to enforce

$(I+\Delta tG)_{ll}\geq 1-\Delta t\left|\Upsilon_{ll}\right|\geq 0.$

∎∎

It is clear that a discrete filter that is positive and conservative is also
dissipative by reiterating through the arguments given above for the
continuous case. Sadly, it is also true that while in the continuous case the
filter which is generated by the second derivative, i.e. the heat kernel, is
positive, the second derivative discretised in our DG method is not a positive
generator and also does not generate a positive filter directly. We will
therefore show how to design a generator generating an approximation of the
heat kernel for forward Euler steps, thereby even allowing to prove the
dissipativity of the entropy dissipation operator for finite time steps. The
basis will be the heat equation with varying heat conductivity $\alpha(x)$
[19]

${\frac{\partial{u}}{\partial{t}}}=\sum_{k=1}^{n}{\frac{\partial{\leavevmode\nobreak\
}}{\partial{x_{k}}}}\alpha(x){\frac{\partial{u}}{\partial{x_{k}}}},\quad\alpha{\frac{\partial{u}}{\partial{n}}}\Bigg{|}_{\partial
T}=0$

on the (reference) element in conjunction with Neumann boundary conditions.
The Neumann boundary conditions enforce the conservation of the resulting
solution operator as any change of the cell mean values must happen through
the numerical flux of the basic DG method. Discretising this problem [19] with
the nodal basis of the basic DG method that is a continuous Galerkin method in
this case because a single element is considered, yields

${\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}=-M^{-1}Qu^{T},\quad
Q_{kl}=\int_{T}\left({\frac{\partial{\varphi_{k}}}{\partial{x}}}\middle)\cdot\middle(\alpha(x){\frac{\partial{\varphi_{l}}}{\partial{x}}}\right)\,\mathrm{d}x.$
(11)

As noted before, in general there exists no $\Delta t>0$ where
$\operatorname{\mathrm{I}}+\Delta t(-M^{-1}Q)$ is a positive operator because
the negative elements in $(-M^{-1}Q)$ prohibit it from being a positive
generator. Yet the following theorem shows that the exact ODE solution to this
problem for a $t>0$ big enough is in fact eligible as a filter.

######  Theorem 2.

If the quadrature $\omega$ is exact on $V^{T}$, the solution of (11) for a
positive initial condition $u_{0}\in\mathbb{R}^{p+1}$ satisfies for all $t>0$

  * •

$\sum_{k=1}^{p+1}\omega_{k}u_{k}(t)=\sum_{k=1}^{p+1}\omega_{k}u_{k}(0)$
(Conservation)

  * •

$u_{k}(t)=C_{kl}(t)u_{l}(0)$ with $\forall
k\in\\{1,\dots,p+1\\}:\quad\sum_{l=1}^{p+1}C_{kl}=1$ (averaging property)

Further, for a $t>0$ big enough it follows $\forall
k\in\\{1,\dots,p+1\\}:\quad u_{k}(t)\geq 0$.

###### Proof.

Entering $v=1$ into the weak form results in

$\int_{T}{\frac{\partial{u}}{\partial{t}}}\,\mathrm{d}x=-\int_{T}{\frac{\partial{1}}{\partial{x}}}{\frac{\partial{u}}{\partial{x}}}\,\mathrm{d}x=0.$

As the quadrature is exact for the basis functions the same follows for the
discretisation, and this shows the conservation. The matrix $C$ used to
describe the solution has the explicit form [25, sec. 34]

$u(t)=\underbrace{\mathrm{e}^{-tM^{-1}Q}}_{C(t)}u(0).$

Multiplying this matrix with the vector $v\in V$ representing the function $1$
from the right reveals

$M^{-1}Qv=0\implies\mathrm{e}^{-tM^{-1}Q}v=\operatorname{\mathrm{I}}v=v.$

This already shows the second result as the nodal representation of $C(t)$
must have unit row sum. The matrix $-Q$ is negative semi-definite, the $v$
vector is in its null space. If another linearly independent $u\in V$ would be
in its null space it would follow

$\left(u\middle)\cdot\middle(Qu\right)=\left\langle{{\frac{\partial{u}}{\partial{x}}}},{\alpha{\frac{\partial{u}}{\partial{x}}}}\right\rangle_{T}=\int_{T}\alpha(x)\left|{\frac{\partial{u}}{\partial{x}}}\right|^{2}\,\mathrm{d}x=0$

and this is a contradiction to ${\frac{\partial{u}}{\partial{x}}}\neq 0$, as
$u$ was assumed non-constant. Therefore, there exists an orthonormal
eigenvalue decomposition of the discretisation whose eigenvalues, apart from
the constant eigenfunction $\psi_{1}=v$ with eigenvalue $\lambda_{1}=0$, are
bounded away from zero,

$\forall k\in\\{1,\dots,p+1\\}:\quad-M^{-1}Q\psi_{k}=\lambda_{k}\psi_{k}.$

We assume that the eigenvectors are sorted by increasing absolute value of the
corresponding eigenvalues,

$0=\lambda_{1}<\left|\lambda_{2}\right|\leq\left|\lambda_{3}\right|\leq\dots\leq\left|\lambda_{p+1}\right|.$

The solution

$u(t)=\sum_{k=1}^{p+1}\mathrm{e}^{-\lambda_{k}t}\psi_{k}\left\langle{\psi_{k}},{u(0)}\right\rangle$

therefore converges to the average of $u(0)$, as

$\left\|u(t)-\psi_{0}\left\langle{\psi_{0}},{u(0)}\right\rangle\right\|^{2}=\left\|\sum_{k=2}^{p+1}\mathrm{e}^{-\lambda_{k}t}\psi_{k}\left\langle{\psi_{k}},{u(0)}\right\rangle\right\|^{2}\leq\mathrm{e}^{-2\lambda_{1}t}\left\|u_{0}\right\|^{2}$

holds. Because a positive initial condition has a positive average the
solution will converge to this positive average. ∎∎

Using the theorem above we can construct filters $\Upsilon$ simply by
calculating the matrix $C(t)=G$ used in the proof above. This matrix which
maps an initial state onto the solution at time $t$ is always a conservative
filter, and when $t$ is large enough also positive. In the implementation the
suitable $t$ was found using a bisection algorithm. Using
$\Upsilon=(C(t)-\operatorname{\mathrm{I}})/t$ the corresponding generator can
be found. We note in passing that numerous other possibilites exist to define
a positive conservative filter as defined above, but that the method given
above defines a filter than can be associated with viscosity.

######  Lemma 6.

Assume the null space of $G$ consists only of constants. Then for a non-
constant $u$ and a strictly convex entropy $U$ it holds

$\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{Gu}\right\rangle_{T,\omega}<0.$

If $U$ is just convex only »$\leq$« applies in the equation above.

###### Proof.

The discrete dissipativity

$\displaystyle E_{T}(u+\Delta t\lambda Gu)$
$\displaystyle=\sum_{k=1}^{p+1}\omega_{k}U\left(u_{k}+\Delta
t\lambda\sum_{l=1}^{p+1}G_{kl}u_{l}\right)$
$\displaystyle<\sum_{k=1}^{p+1}\omega_{k}\sum_{l=1}^{p+1}G_{kl}U(u_{l})=\sum_{l=1}^{p+1}\omega_{l}U(u_{l})=E_{u,T}$

follows from the positive conservative filter property of $G$ for
$\lambda\Delta t>0$ small enough as in lemma 5 in conjunction with the strict
convexity and Jensens inequality in the strict sense. Let now $\Delta t$ be
fixed and small enough for all $\lambda\in[0,1]$, and denote by
$\varepsilon=E^{T}(u+\Delta tGu)-E^{T}_{u}<0$ the entropy dissipation for
$\lambda=1$. The convexity of $U$ implies

$E_{T}(u+\lambda\Delta tGu)\leq E_{T}(u)+\lambda\left(E_{T}(u+\Delta
tGu)-E_{T,u}\right)=E_{T}(u)+\lambda\varepsilon$

Entering this into the definition of the derivative of $E^{T}$ with respect to
$\lambda$ shows

${\frac{\mathrm{d}E_{T}(u+\lambda\Delta
tGu)}{\mathrm{d}\lambda}}=\lim_{\lambda\to 0}\frac{E_{T}(u+\lambda\Delta
tGu)-E_{T}(u)}{\lambda}\leq\varepsilon,$

and therefore

$\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{Gu}\right\rangle_{T,\omega}=\sum_{k=1}^{p+1}\omega_{k}\left({\frac{\mathrm{d}U}{\mathrm{d}u}}(x_{k})\middle)\cdot\middle((Gu)_{k}\right)={\frac{\mathrm{d}E_{T}(u+\lambda\Delta
tGu)}{\mathrm{d}\lambda}}\frac{1}{\Delta t}=\frac{\varepsilon}{\Delta t}.$

If $U$ is not strictly convex the case $\varepsilon=0$ is possible, reducing
the result to »$\leq$«. ∎∎

The last step consists of selecting a suitable viscosity distribution
$\alpha$, i.e. one that is zero at the endpoints. The standard mollifier

$\alpha(x)=\begin{cases}\mathrm{e}^{1-\frac{1}{1-x^{2}}}&\left|x\right|\leq
1\\\ 0&\left|x\right|>1\end{cases}$

is smooth and zero at the ends of the reference element. Further, even its
derivatives vanish there. It was therefore selected.

###  3.1 Stable computation of the correction size required and timestep
restrictions

After we have calculated the entropy dissipation needed and a suitable
direction $\upsilon=Gu^{T}$ one would guess we only have to calculate
$\lambda$ as in (6) via

$\lambda\geq\frac{\sigma^{T}-\left(\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{{\frac{\mathrm{d}u}{\mathrm{d}t}}}\right\rangle_{T,\omega}-\left(F^{*}_{l}-F^{*}_{r}\right)\right)}{\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{\upsilon}\right\rangle_{T,\omega}}.$

It turns out that this process is significantly more intricate than one would
expect as this computation has to be stable with respect to roundoff errors.
Further, our estimates on the entropy dissipation can only estimate the
entropy dissipation that can take place at the interface between two adjacent
cells, but are not able to give an estimate of how this dissipation is split
between the two cells. Our method of calculating suitable values of
$\lambda^{T}$ therefore consists of two steps. First,

$\lambda_{\mathrm{ED}}^{T}=\max\left(0,-\frac{\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{{\frac{\mathrm{d}u}{\mathrm{d}t}}}\right\rangle_{T,\omega}-\left(F^{*}_{l}-F^{*}_{r}\right)}{\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{\upsilon}\right\rangle_{T,\omega}}\right)$
(12)

is calculated to enforce the per cell entropy dissipativity

$\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{{\frac{\mathrm{d}u}{\mathrm{d}t}}+\lambda_{\mathrm{ED}}^{T}\upsilon}\right\rangle_{T,\omega}\leq
F^{*}_{l}-F^{*}_{r}.$

In a second step a correction to enforce an entropy rate high enough

$\lambda_{\mathrm{ER}}^{\theta}=\max\left(0,\frac{\sigma^{\theta}-\sum_{T\cap\theta\neq\emptyset}\left(\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{{\frac{\mathrm{d}u}{\mathrm{d}t}}+\lambda_{\mathrm{ED}}^{T}\upsilon}\right\rangle_{T,\omega}-(F^{*}_{l,T}-F^{*}_{r,T})\right)}{\sum_{\theta\cap
T\neq\emptyset}\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{\upsilon}\right\rangle_{T,\omega}}\right)$
(13)

is determined for all $\theta\in\Theta$. Both corrections are then added
together

$\lambda^{T}_{\Sigma}=\lambda^{T}_{\mathrm{ED}}+\sum_{\theta\cap
T\neq\emptyset}\lambda^{\theta}_{\mathrm{ER}}$

for all cells $T\in\mathcal{T}$. Round-off errors tend to influence the
calculation out of two reasons. The division by
$\left\langle{{\frac{\mathrm{d}U}{\mathrm{d}u}}},{\upsilon}\right\rangle$ in
equation (12) and (13) can approach a division by zero for a solution
approaching a constant in the cell, as $\upsilon\to 0$ follows in this case.
Further, we saw in lemma (1) that the entropy inequality predictor can vanish
with a high order for smooth solutions, and an accurate DG scheme will also
have a vanishing entropy error vanishing with a high order. The difference of
these two values, i.e. the denominator of the fraction above, will in general
not vanish that fast because round-off in the difference becomes important.
Therefore $\lambda$ will, for highly resolved smooth solutions, be to big
because round-off errors propagate into the calculation. Our solution to this
problem is to calculate

$\lambda_{(\cdot)}^{T}=\max\left(\frac{ab}{b^{2}+c^{2}},0\right),\text{
instead of }\frac{a}{b}$

every time a $\lambda$ is calculated by a division in the procedure above.
Here, $a$ shall be the nominator, $b$ shall be the denominator and $c$ shall
be a suitable bound on the round-off error, a constant small with respect to
$a,b$ but large with respect to the machine precision. In our implementation
this is selected as $c=\sqrt{10^{-16}}$, i.e. the square root of the machine
precision for a solution scaled to be of unit magnitude. The addition of $c$
can be seen as the one-dimensional version of Tikhonov regularization [23].
Clipping the calculation of $\lambda$ at $0$ ensures that if $a$ or $b$ become
negative from rounding errors $\lambda$ will not become negative, i.e.
$\lambda\upsilon$ will not be antidissipative. In a last step,

$\lambda^{T}=\min\left(\lambda_{\mathrm{max}},\lambda^{T}_{\Sigma}\right),$

the upper limit $\lambda_{\text{max}}$ is introduced for stability reasons as
we want to enforce stability of

${\frac{\mathrm{d}u^{T}}{\mathrm{d}t}}=\upsilon=\lambda^{T}Gu^{T}.$ (14)

If a Runge-Kutta time integration method can be written as convex combination
of forward Euler steps, i.e. is Strong Stability Preserving (SSP) [16, 31, 32]
and the time-steps satisfy $\Delta t\lambda\leq 1$ during every Euler step,
the lemma 5 allows us to show that the solution is also entropy dissipative in
the discrete case. If the time integration method used is just a conditionally
stable Runge-Kutta method [9, 41] we are interested in limiting the operator
norm of $\Delta t\left\|\lambda G\right\|\leq R$ in order to at least avoid a
linear instability. The exact size depends on the time integration methods’
stability region as we would like to fit the half-circle

$C=\\{z\in\mathbb{C}|\left\|z\right\|\leq R\wedge\operatorname{im}z\leq 0\\}$

into the stability region of the method.

##  4 Numerical tests

Our tests will be carried out for the Euler equations of gas dynamics in
conservation form [17]

$u=(\rho,\rho v,E)\quad f(\rho,\rho v,E)=\begin{bmatrix}\rho v\\\ \rho
v^{2}+p\\\ v(E+p)\end{bmatrix}\quad p=(\gamma-1)\left(E-\frac{1}{2}\rho
v^{2}\right)$

in conjunction with the physical entropy [38, 17]

$U(\rho,\rho v,E)=-\rho S\quad F(\rho,\rho v,E)=-\rho vS\quad
S=\ln(p\rho^{-\gamma}).$

The tests below will focus on the cases $p=3$ and $p=7$ as the latter are
popular in applications because they amount to $4$ and $8$ nodes, suitable for
SIMD processor instructions. Results for values in between are essentially
interpolatory to the ones reported for $p=3$ and $p=7$ and the source code is
available to carry out tests for all values $p>0$. Time integration will be
carried out using the SSPRK(4, 3) method for most solutions, while the
convergence analysis for $p=7$ below will use the Hairer-Wanner DOPRI8 method,
to achieve the needed convergence speed of the time integration. In all images
below the ansatz functions of all cells are shown without any post-processing.

Property | Tested Solver | Reference
---|---|---
Type | DDG | first order FV
Intercell Flux | local Lax-Friedrichs | Lax-Friedrichs
CFL number | $0.1/(p^{2}+p)$ | $0.5$
Time Integration | SSPRK(4,3), DoPri8 | Forward Euler
Dissipation |  $G$ from sec. 3 | Built-in
$\lambda_{\mathrm{max}}$ | $\frac{1}{\Delta t}$ | Does not apply
Number of Cells | 13, 25, 50, 100, 200 | $3\cdot 10^{4}$
Table 2: Used schemes in the numerical tests

###  4.1 Shock tube tests

First, a series of shock tube tests was done to highlight the effectivity of
the entropy correction in shock calculations as this is the primary aim of
this publication. The first initial condition [39, 32, Problem I, Section
4.3.3 and Problem 6a]) is

$\displaystyle\rho_{0}(x,0)=\begin{cases}1,\\\ 0.125,\end{cases}\quad
v_{0}(x,0)=\begin{cases}0,\\\
0,\end{cases}p_{0}(x,0)=\begin{cases}1.0,&x<5,\\\ 0.1,&x\geq 5.\end{cases}$

Our second shock tube is the time-evolution of the following Riemann problem
[32, Problem 6b])

$\displaystyle\rho_{0}(x,0)=\begin{cases}0.445,\\\ 0.5,\end{cases}\quad
v_{0}(x,0)=\begin{cases}0.698,\\\
0,\end{cases}p_{0}(x,0)=\begin{cases}3.528,&x<5.0,\\\ 0.571,&x\geq
5.0.\end{cases}$

(a) Density, $25$ cells

(b) Density, $100$ cells

(c) Moment density, $25$ cells

(d) Moment density, $100$ cells

(e) Energy density, $25$ cells

(f) Energy density, $100$ cells

Figure 5: Shock tube 1 at $t=1.8$ with $25$ cells corresponding to $100$
degrees of freedom and $100$ cells corresponding to $400$ degrees of freedom
(p=3).

(a) Density, $13$ cells

(b) Density, $100$ cells

(c) Moment density, $13$ cells

(d) Moment density, $100$ cells

(e) Energy density, $13$ cells

(f) Energy density, $100$ cells

Figure 6: Shock tube 1 at $t=1.8$ with $13$ cells corresponding to $104$
degrees of freedom and $100$ cells corresponding to $800$ degrees of freedom
(p=7).

(a) Density, $25$ cells

(b) Density, $100$ cells

(c) Moment density, $25$ cells

(d) Moment density, $100$ cells

(e) Energy density, $25$ cells

(f) Energy density, $100$ cells

Figure 7: Shock tube 2 at $t=1.2$ with $25$ cells corresponding to $100$
degrees of freedom and $100$ cells corresponding to $400$ degrees of freedom
(p=3).

(a) Density, $13$ cells

(b) Density, $100$ cells

(c) Moment density, $13$ cells

(d) Moment density, $100$ cells

(e) Energy density, $13$ cells

(f) Energy density, $100$ cells

Figure 8: Shock tube 2 at $t=1.2$ with $25$ cells corresponding to $100$
degrees of freedom and $100$ cells corresponding to $400$ degrees of freedom
(p=7).

The shock tube tests were always carried out for two different numbers of
cells. First for $N=N_{\text{typ}}/(p+1)$ cells, where $N_{typ}=100$ is the
usual number of cells used in comparisons for Finite-Volume methods. This was
done so that the same number of degrees of freedom has to be saved. The
results look satisfactory and highlight the effectivity of the method in
figures 5, 6, 7, 8. All shocks are sharp and concentrated to less than one
cell width. Yet, only slight overshoots and oscillations are visible directly
around the shocks. These distortions are confined to the cell directly next to
the shock. Contact discontinuities are slightly smeared over one cell, but
after they have been smeared to this width no additional smearing takes place.
The computational complexity per timestep is still low as no recovery stencil
selection has to be carried out and only $1/(p+1)$ times the number of two-
point fluxes need to be evaluated. Because some other publications use 100
cells also for DG methods we carried out the tests once more for $N=100$
cells, amounting to $400$ and $800$ degrees of freedom for orders $p=3$ and
$p=7$.

###  4.2 Numerical validation of the entropy rate criterion

To verify the entropy rate criterion the total entropy of the solution to the
first shock tube above was compared to the solution calculated by a Lax-
Friedrichs scheme with $3\cdot 10^{4}$ cells. Similar comparisons were carried
out in [20, 21, 22]. Please note that the Godunov solver used previously was
swapped for a LF scheme to evade the need for an exact Riemann solver. This is
also supported by our finding in lemma 1 and corollary 1 as a Lax-Friedrichs
solution therefore has to comply with the entropy rate criterion. A scheme
should in these comparisons have the same entropy dissipation rate (in the
limit) as the Lax-Friedrichs scheme in the limit. Comparisons for orders $p=3$
and $p=7$ in figure 9 show that this seems to be the case. The DG scheme
always has an entropy that lies below the entropy of the LF scheme. As the
entropy inequality for vanishing viscosity solutions is also desirable it was
also verified on a per-cell basis. We just note that the small positive
violations in figure 9 are of the same magnitude as the precision achievable
during the calculation of $\lambda$ using our procedure with double precision
floats.

(a) Plot of the total entropy ($p=3$)

(b) Plot of the total entropy $(p=7)$

(c) Logarithm of the positive violation of the entropy equality ($p=3$)

(d) Logarithm of the positive violation of the entropy equality ($p=7$)

(e) Logarithm of the negative violation of the entropy equality ($p=3$)

(f) Logarithm of the negative violation of the entropy equality ($p=7$)

Figure 9: Entropy tests for the first shock tube.

###  4.3 Shu-Osher test

To showcase a combination of shocks and smooth areas the well established
shock-sine interaction problem from [32, Problem 8] was tested. The initial
conditions are given by

$\displaystyle\rho_{0}(x,0)=\begin{cases}3.857153\\\
1+\varepsilon\sin(5x)\end{cases}\quad v_{0}(x,0)=\begin{cases}2.629\\\
0\end{cases}p_{0}(x,0)=\begin{cases}10.333&x<1\\\ 1&x\geq 1\end{cases}.$

The parameter $\varepsilon$ was set to the canonical value of
$\varepsilon=0.2$.

(a) $p=3$, $50$ cells.

(b) $p=7$, $50$ cells.

(c) $p=3$, $100$ cells.

(d) $p=7$, $100$ cells.

(e) $p=3$, $200$ cells.

(f) $p=7$, $200$ cells.

Figure 10: Shu-Osher test for $50,100$ and $200$ cells and order $p=3$ and
$p=7$ and therefore $200,400,800$ and $1600$ degrees of freedom at $t=1.8$.

The results look satisfactory already when only $N=100$ cells are used in the
calculation. Yet, we note that this already corresponds to $400$ and $800$
degrees of freedom for the selected orders. When $N=200$ cells are used the
solution is nearly indistinguishable from the reference solution.

###  4.4 Convergence Analysis

While the main aim of our modification was to devise a new DG scheme usable
for shock-capturing calculations the scheme also converges with high order of
accuracy for smooth solutions in our experiments. As an example the solution
of

$\rho_{0}(x,0)=3.857153+\varepsilon(x)\sin(2x),\quad v_{0}(x,0)=2.0,\quad
p_{0}(x,0)=10.33333,$

with

$\varepsilon(x)=\mathrm{e}^{(x-3)^{2}}.$

and periodic boundary conditions was calculated using our modified DG method.
The analytical solution for this test problem is

$\rho(x,t)=3.857153+\varepsilon(x-2t)\sin(2x-4t),\quad v(x,t)=2.0,\quad
p(x,t)=10.33333,$

with suitable periodic boundary conditions.

(a) p = 3, $\mathrm{L}^{1}$ norm

(b) p = 7, $\mathrm{L}^{2}$ norm

(c) p = 7, $\mathrm{L}^{1}$ norm

(d) p = 7, $\mathrm{L}^{2}$ norm

Figure 11: Convergence Analysis for $p=3$ and $p=7$.

After the solution was calculated for
$N=\\{10,15,20,25,30,40,50,60,70,80,90,100\\}$ cells for $p=3$ and with the
same stepping up to $50$ cells for $p=7$ up to $t=5$ the $\mathrm{L}^{1}$ and
$\mathrm{L}^{2}$ errors were calculated. The convergence in figure 11 seems to
take place with too high an order for the ansatz polynomials used. The reason
for this could be that the accuracy of the basic scheme is significantly
higher for these solutions than the accuracy of the corrected scheme, because
the entropy dissipation estimate still falsely reports high amounts of entropy
dissipation. When the grid is refined the entropy dissipation estimate
converges with a higher speed than the basic scheme following lemma 3 and
because the error introduced to enforce the dissipation dominates a higher
convergence speed than expected is observed.

###  4.5 Timestep Analysis

An important result of any modification to a basic scheme can be an impact on
the allowed timestep size. In the first part of this publication [21] this
influence was tested by measuring the maximal timestep possible before a blow-
up occurs. This was done once more.

Figure 12: Maximal timestep sizes for $p=3$ and $p=7$ and the shock tube .

The maximal timestep possible for the first shock tube for orders $p=3$ and
$p=7$ is shown in figure 12. Obviously this timestep is acceptable and when
corrected for the larger maximal wave speed of the Riemann problem used for
testing, larger than the timestep reported in the previous part, highlighting
the superiority of the new dissipation direction.

##  5 Conclusion

The method described in [21] to enforce an entropy rate criterion for DG
methods was improved. By using a direct indicator for the entropy dissipation
the error indicator used before could be replaced, resulting in a lower
dissipation in situations like contact discontinuities. For smooth solutions
this new method to quantify the amount of dissipation needed converges
significantly faster to zero than the error estimate used before, and
therefore allows us to recover the convergence speed of the basic DG scheme
that was reduced by one degree before. Further, the direct quantification of
the entropy dissipation needed allowed us to consider different dissipation
directions, especially combining smoothing and dissipation and therefore
bridging into the field of modal filtering. The effectivity of the refined
method was demonstrated for the Euler system of gas dynamics. The method is
not only high order accurate but also able to handle shocks, contact
discontinuities, and rarefactionwaves. The next logical steps can be the
application to two-dimensional problems, the application of the designed
entropy inequality predictors to other schemes like continuous Galerkin and
Spectral Volume schemes, where several adjustments will be needed, and
revisiting the splitting into a fully discrete scheme already explored in
[21]. The presented method to estimate the entropy dissipation needed could
also be used with artificial viscosity shock-capturing as for example
described in [13].

##  6 Competing Interests

The author has no relevant financial or non-financial interests to disclose.

##  7 Data Availability

The commented implementation of the schemes is available under
https://github.com/simonius/dgdafermos.

##  8 Bibliography

## References

  * Ben-Artzi and Falcovitz [2011] Matania Ben-Artzi and Joseph Falcovitz.  _Generalized Riemann problems in computational fluid dynamics._ , volume 11 of _Camb. Monogr. Appl. Comput. Math._ Cambridge: Cambridge University Press, reprint of the 2003 hardback ed. edition, 2011.  ISBN 978-0-521-17327-8.  doi: 10.1017/CBO9780511546785.
  * Chavent and Cockburn [1989] Guy Chavent and Bernardo Cockburn.  The local projection $P^{0}-P^{1}$-discontinuous-Galerkin finite element method for scalar conservation laws.  _RAIRO, Modélisation Math. Anal. Numér._ , 23(4):565–592, 1989.  ISSN 0764-583X.  doi: 10.1051/m2an/1989230405651.
  * Chen and Shu [2017] Tianheng Chen and Chi-Wang Shu.  Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws.  _J. Comput. Phys._ , 345:427–461, 2017.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2017.05.025.
  * Cockburn and Shu [1989] Bernardo Cockburn and Chi-Wang Shu.  TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II: General framework.  _Math. Comput._ , 52(186):411–435, 1989.  ISSN 0025-5718.  doi: 10.2307/2008474.
  * Cockburn and Shu [2001] Bernardo Cockburn and Chi-Wang Shu.  Runge–Kutta discontinuous Galerkin methods for convection-dominated problems.  _J. Sci. Comput._ , 16(3):173–261, 2001.  ISSN 0885-7474.  doi: 10.1023/A:1012873910884.
  * Dafermos [1972] Constantine M. Dafermos.  The entropy rate admissibility criterion for solutions of hyperbolic conservation laws.  _Journal Of Differential Equations_ , pages 202–212, 1972.
  * Dafermos [2009] Constantine M. Dafermos.  A variational approach to the Riemann problem for hyperbolic conservation laws.  _Discrete Contin. Dyn. Syst._ , 23(1-2):185–195, 2009.  ISSN 1078-0947.  doi: 10.3934/dcds.2009.23.185.
  * Dafermos [2016] Constantine M. Dafermos.  _Hyperbolic conservation laws in continuum physics_ , volume 325\.  Berlin: Springer, 2016.  ISBN 978-3-662-49449-3; 978-3-662-49451-6.  doi: 10.1007/978-3-662-49451-6.
  * Dahlquist [1963] Germund G. Dahlquist.  A special stability problem for linear multistep methods.  _BIT, Nord. Tidskr. Inf.-behandl._ , 3:27–43, 1963.  ISSN 0006-3835.  doi: 10.1007/BF01963532.
  * Evans [2010] Lawrence C. Evans.  _Partial differential equations_ , volume 19 of _Grad. Stud. Math._ Providence, RI: American Mathematical Society (AMS), 2nd ed. edition, 2010\.  ISBN 978-0-8218-4974-3.
  * Feireisl [2014] Eduard Feireisl.  Maximal dissipation and well-posedness for the compressible Euler system.  _J. Math. Fluid Mech._ , 16(3):447–461, 2014\.  ISSN 1422-6928.  doi: 10.1007/s00021-014-0163-8.
  * Gassner [2013] Gregor J. Gassner.  A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods.  _SIAM J. Sci. Comput._ , 35(3):a1233–a1253, 2013\.  ISSN 1064-8275.  doi: 10.1137/120890144.
  * Glaubitz et al. [2019] J. Glaubitz, A. C. jun. Nogueira, J. L. S. Almeida, R. F. Cantão, and C. A. C. Silva.  Smooth and compactly supported viscous sub-cell shock capturing for discontinuous Galerkin methods.  _J. Sci. Comput._ , 79(1):249–272, 2019.  ISSN 0885-7474.  doi: 10.1007/s10915-018-0850-3.
  * Glaubitz et al. [2018] Jan Glaubitz, Philipp Öffner, and Thomas Sonar.  Application of modal filtering to a spectral difference method.  _Math. Comput._ , 87(309):175–207, 2018.  ISSN 0025-5718.  doi: 10.1090/mcom/3257.
  * Godlewski and Raviart [1991] Edwige Godlewski and Pierre-Arnaud Raviart.  _Hyperbolic Sytems of Conservation Laws_.  ellipses, 1991.
  * Gottlieb et al. [2001] Sigal Gottlieb, Chi-Wang Shu, and Eitan Tadmor.  Strong stability-preserving high-order time discretization methods.  _SIAM Review_ , 43, 05 2001.  doi: 10.1137/S003614450036757X.
  * Harten [1983] Amiram Harten.  On the symmetric form of systems of conservation laws with entropy.  _Journal of Computational Physics_ , 49:151–164, 1983.
  * Harten et al. [1983] Amiram Harten, Peter D. Lax, and Bram van Leer.  On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.  _SIAM Rev._ , 25:35–61, 1983.  ISSN 0036-1445.  doi: 10.1137/1025002.
  * Johnson [2009] Claes Johnson.  _Numerical solution of partial differential equations by finite element method._ Mineola, NY: Dover Publications, reprint of the 1987 English ed. edition, 2009.  ISBN 978-0-486-46900-3.
  * Klein [2022] Simon-Christian Klein.  Using the Dafermos entropy rate criterion in numerical schemes.  _BIT_ , 62(4):1673–1701, 2022.  ISSN 0006-3835.  doi: 10.1007/s10543-022-00927-x.
  * Klein [2023] Simon-Christian Klein.  Stabilizing discontinuous galerkin methods using dafermos’ entropy rate criterion: I—one-dimensional conservation laws.  _Journal of Scientific Computing_ , 95(2):55, 2023\.
  * Klein and Sonar [2023] Simon-Christian Klein and Thomas Sonar.  Entropy-aware non-oscillatory high-order finite volume methods using the dafermos entropy rate criterion, 2023.  URL https://arxiv.org/abs/2302.08971.
  * Kress [1998] Rainer Kress.  _Ill-Conditioned Linear Systems_ , pages 77–92.  Springer New York, New York, NY, 1998.  ISBN 978-1-4612-0599-9.  doi: 10.1007/978-1-4612-0599-9˙5.  URL https://doi.org/10.1007/978-1-4612-0599-9_5.
  * Lax [1971] Peter D. Lax.  Shock waves and entropy.  _Contributions to Nonlinear Functional Analysis_ , pages 603–634, 1971.
  * Lax [2002] Peter D. Lax.  _Functional Analysis_.  Wiley Interscience, 2002.
  * Lax et al. [1976] Peter D. Lax, Samuel Burstein, and Anneli Lax.  _Calculus with applications and computing. Vol. I_.  Undergraduate Texts Math. Springer, Cham, 1976.
  * Luo et al. [2007] Hong Luo, Joseph D. Baum, and Rainald Löhner.  A hermite weno-based limiter for discontinuous galerkin method on unstructured grids.  _Journal of Computational Physics_ , 225(1):686–713, 2007.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2006.12.017.  URL https://www.sciencedirect.com/science/article/pii/S0021999106006164.
  * Ranocha [2019] Hendrik Ranocha.  Mimetic properties of difference operators: Product and chain rules as for functions of bounded variation and entropy stability of second derivatives.  _BIT Numerical Mathematics_ , 59(2):547–563, 06 2019.  doi: 10.1007/s10543-018-0736-7.
  * Ranocha et al. [2018] Hendrik Ranocha, Jan Glaubitz, Philipp Öffner, and Thomas Sonar.  Stability of artificial dissipation and modal filtering for flux reconstruction schemes using summation-by-parts operators.  _Appl. Numer. Math._ , 128:1–23, 2018.  ISSN 0168-9274.  doi: 10.1016/j.apnum.2018.01.019.
  * Rudin [1966] Walter Rudin.  Real and complex analysis.  McGraw-Hill Series in Higher Mathematics. New York etc.: McGraw-Hill Book Company. xi, 412 p. (1966)., 1966.
  * Shu and Osher [1988] Chi-Wang Shu and Stanley Osher.  Efficient implementation of essentially non-oscillatory shock-capturing schemes.  _Journal of Computational Physics_ , 77:439–471, 1988.
  * Shu and Osher [1989] Chi-Wang Shu and Stanley Osher.  Efficient implementation of essentially non-oscillatory shock-capturing schemes II.  _Journal of Computational Physics_ , 83:439–471, 1989.
  * Sonar [2016] T. Sonar.  Chapter 3 - classical finite volume methods.  In Rémi Abgrall and Chi-Wang Shu, editors, _Handbook of Numerical Methods for Hyperbolic Problems_ , volume 17 of _Handbook of Numerical Analysis_ , pages 55–76. Elsevier, 2016.  doi: https://doi.org/10.1016/bs.hna.2016.09.005.  URL https://www.sciencedirect.com/science/article/pii/S157086591630014X.
  * Sonntag and Munz [2014] Matthias Sonntag and Claus-Dieter Munz.  Shock capturing for discontinuous Galerkin methods using finite volume subcells.  In Jürgen Fuhrmann, Mario Ohlberger, and Christian Rohde, editors, _Finite Volumes for Complex Applications VII-Elliptic, Parabolic and Hyperbolic Problems_ , pages 945–953, Cham, 2014. Springer International Publishing.  ISBN 978-3-319-05591-6.
  * Tadmor [1984a] Eitan Tadmor.  The large-time behavior of the scalar, genuinely nonlinear Lax-Friedrichs scheme.  _Math. Comput._ , 43:353–368, 1984a.  ISSN 0025-5718.  doi: 10.2307/2008281.
  * Tadmor [1984b] Eitan Tadmor.  Numerical viscosity and the entropy condition for conservative difference schemes.  _Math. Comput._ , 43:369–381, 1984b.  ISSN 0025-5718.  doi: 10.2307/2008282.
  * Tadmor [1987] Eitan Tadmor.  The numerical viscosity of entropy stable schemes for systems of conservation laws.  _Mathematics of Computation_ , 49:91–103, 1987.
  * Tadmor [2003] Eitan Tadmor.  Entropy stability theory for difference approximations of nonlinear conservation laws and related time dependent problems.  _Acta Numerica_ , pages 451–512, 2003.
  * Toro [2009] Eleuterio F. Toro.  _Riemann solvers and numerical methods for fluid dynamics. A practical introduction_.  Berlin: Springer, 2009.  ISBN 978-3-540-25202-3; 978-3-540-49834-6.  doi: 10.1007/b79761.
  * Wang [2002] Z. J. Wang.  Spectral (finite) volume method for conservation laws on unstructured grids. Basic formulation.  _J. Comput. Phys._ , 178(1):210–251, 2002.  ISSN 0021-9991.  doi: 10.1006/jcph.2002.7041.  URL semanticscholar.org/paper/4e9aef2784954f33505124a59c0e69d2248a2c96.
  * Wanner et al. [1978] Gerhard Wanner, Ernst Hairer, and Syvert P. Nørsett.  Order stars and stability theorems.  _BIT, Nord. Tidskr. Inf.-behandl._ , 18:475–489, 1978.  ISSN 0006-3835.  doi: 10.1007/BF01932026.
  * Zhu and Qiu [2011] Jun Zhu and Jianxian Qiu.  Local DG method using WENO type limiters for convection–diffusion problems.  _Journal of Computational Physics_ , 230(11):4353–4375, 2011.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2010.03.023.  URL https://www.sciencedirect.com/science/article/pii/S0021999110001336.  Special issue High Order Methods for CFD Problems.