∎

11institutetext: Ge Yan, Graduate Student 22institutetext: Rensselaer
Polytechnic Institute
22email<EMAIL_ADDRESS>33institutetext: Sharanjeet Kaur, Graduate
Student 44institutetext: Rensselaer Polytechnic Institute
44email<EMAIL_ADDRESS>55institutetext: Jeffrey W. Banks, Associate Professor
66institutetext: Rensselaer Polytechnic Institute
66email<EMAIL_ADDRESS>77institutetext: Jason E. Hicken, Associate Professor
88institutetext: Rensselaer Polytechnic Institute
88email<EMAIL_ADDRESS>
# Entropy-stable discontinuous Galerkin difference methods for hyperbolic
conservation laws

Ge Yan  Sharanjeet Kaur  Jeffery W. Banks  Jason E. Hicken

(Received: date / Accepted: date)

###### Abstract

The paper describes the construction of entropy-stable discontinuous Galerkin
difference (DGD) discretizations for hyperbolic conservation laws on
unstructured grids. The construction takes advantage of existing theory for
entropy-stable summation-by-parts (SBP) discretizations. In particular, the
paper shows how DGD discretizations — both linear and nonlinear — can be
constructed by defining the SBP trial and test functions in terms of
interpolated DGD degrees of freedom. In the case of entropy-stable
discretizations, the entropy variables rather than the conservative variables
must be interpolated to the SBP nodes. A fully-discrete entropy-stable scheme
is obtained by adopting the relaxation Runge-Kutta version of the midpoint
method. In addition, DGD matrix operators for the first derivative are shown
to be dense-norm SBP operators. Numerical results are presented to verify the
accuracy and entropy-stability of the DGD discretization in the context of the
Euler equations. The results suggest that DGD and SBP solution errors are
similar for the same number of degrees of freedom. Finally, an investigation
of the DGD spectra shows that spectral radius is relatively insensitive to
discretization order; however, the high-order methods do suffer from the
linear instability reported for other entropy-stable discretizations.

###### Keywords:

Galerkin difference summation-by-parts entropy stable unstructured grid

††journal: Journal of Scientific Computing

## Mathematics Subject Classification

65M60, 65M70, 65M12

##  1 Introduction

Many studies have demonstrated that high-order discretizations can simulate
hyperbolic conservation laws with greater efficiency than first- and second-
order discretizations. However, high-order methods tend to be less robust, and
this has hindered their widespread adoption. The need for robust and efficient
high-order discretizations motivates this work on entropy-stable discontinuous
Galerkin difference methods.

Our interest in the Galerkin difference (GD) family of methods stems from the
attractive properties of their underlying basis functions. For example, in the
original GD method proposed by Banks and Hagstrom Banks:GD2016 , the solution
is expressed in terms of piecewise continuous functions that extend over
multiple elements. Consequently, the number of GD degrees of freedom remains
constant as the polynomial degree increases, similar to conventional finite-
difference and finite-volume methods. Furthermore, like finite-difference
methods, GD schemes have time-step restrictions that are relatively
insensitive to their order of accuracy.

The GD method was recently generalized to discontinuous basis functions Banks-
DGD-paper , and this discontinuous Galkerin difference (DGD) method is the
starting point for the discretization considered in this work. Building on
Banks-DGD-paper , we present an entropy-stable formulation of DGD to address
robustness for discretizations of symmetrizable hyperbolic systems. In
addition, unlike Banks-DGD-paper , we consider unstructured grids and
construct basis function stencils following the approach presented by Li et
al. Li2019 .

To derive the entropy-stable DGD method, we leverage the existing theory for
entropy-stable summation-by-parts (SBP) methods. Over the past decade,
researchers have used SBP operators to construct high-order entropy-stable
finite-difference fisher:thesis ; fisher:2013 ; Fisher2013discretely , finite-
element Chan2018discretely , and spectral-element-type methods
Carpenter2014entropy ; parsani:2016 ; Gassner2016well ; Ranocha2018stability ;
Friedrich2018entropy ; Friedrich2019entropy ; Shadpey2020entropy ;
Rojas2021robustness . Entropy-stable discretizations provide a form of
nonlinear stability and, therefore, robustness by ensuring the entropy is non-
increasing111Entropy here refers to mathematical entropy..

In summary, the objectives of this work are to present a framework for
constructing entropy-stable DGD discretizations and to study the properties of
the resulting discretization. The principle contributions are listed below:

  * •

DGD discretizations of linear and nonlinear hyperbolic systems can be
constructed using diagonal-norm SBP operators.

  * •

DGD matrix operators are, themselves, dense-norm SBP operators.

  * •

When applied to entropy-stable DGD semi-discretizations, explicit time
marching schemes require the solution of a coupled nonlinear algebraic system
and, therefore, do not offer a significant computational advantage over
implicit methods.

The roadmap of the paper is as follows. Section 2 reviews the DGD method,
including the construction of the DGD basis functions and the semi-
discretization of the two-dimensional linear advection problem. Section 3
reviews multi-dimensional SBP operators and shows how they can be used in the
DGD semi-discretization of linear advection. This section also proves that DGD
operators are dense-norm SBP operators. In Section 4, we review generic
entropy-conservative/stable SBP discretizations, and then develop the entropy-
stable DGD discretization. Section 5 presents numerical experiments in order
to verify the accuracy and stability properties of the DGD discretization, as
well as characterize its efficiency. Finally, Section 6 concludes this study
with a summary and a discussion of potential future developments.

##  2 The discontinuous Galerkin difference method on unstructured grids

Galerkin difference methods are finite element methods that use piecewise
polynomial basis functions. One of our objectives in this section is to
familiarize readers with these somewhat unconventional basis functions. We
also highlight some important differences that arise when DGD basis functions
are constructed on unstructured versus tensor-product grids and when the
polynomial degree is not even. We conclude this section by illustrating how
the DGD basis functions are used to discretize the constant-coefficient linear
advection equation.

###  2.1 DGD basis functions

Let $T_{h}=\\{\Omega_{k}\\}_{k=1}^{K}$ denote a tesselation of a closed and
bounded domain $\Omega$ into $K$ elements, in which element $k$ has subdomain
$\Omega_{k}$, boundary $\Gamma_{k}$, and centroid $\tilde{\bm{x}}_{k}$. The
stencil
$N_{k}=\\{\nu_{1},\nu_{2},\ldots,\nu_{n_{k}}\\}\subset\\{1,2,3,...K\\}$
represents the ordered set of elements whose basis functions are nonzero on
element $k$ and, therefore, influence the discrete solution over $\Omega_{k}$.
We will sometimes refer to the stencil $N_{k}$ as the patch for element $k$.
The number of elements in the stencil is denoted by $n_{k}$.

Consider a bounded function $\mathcal{U}(\bm{x})$ and let
$\tilde{u}_{\nu}=\mathcal{U}(\tilde{\bm{x}}_{\nu})$ denote its value at the
centroid of element $\nu$. Then the DGD interpolation of $\mathcal{U}$, which
we denote by $\mathcal{U}_{h}(\bm{x})$, takes the following form on element
$k$:

$\mathcal{U}_{h}(\bm{x})=\sum_{\nu=1}^{K}\tilde{u}_{\nu}\phi_{\nu}(\bm{x})=\sum_{\nu\in
N_{k}}\tilde{u}_{\nu}\phi_{\nu}(\bm{x}),\quad\bm{x}\in\Omega_{k},$ (1)

where $\phi_{\nu}(\bm{x})$ is the piecewise polynomial basis associated with
element $\nu$ whose construction is the focus of this section. To help readers
gain some intuition about DGD basis functions, some one dimensional examples
of $\phi_{\nu}(\bm{x})$ are illustrated in Figure 1.

In general, the DGD basis functions can be expressed as a linear combination
of some (standard) polynomial basis functions on each element. Specifically,
the basis of the $i$th neighbour in the stencil $N_{k}$ can be written as

$\phi_{\nu_{i}}(\bm{x})=\sum_{j=1}^{n_{p}}c_{ji}\mathcal{V}_{j}(\bm{x}),\qquad\forall\bm{x}\in\Omega_{k},$
(2)

where $\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}$ is a given basis for
$\mathbb{P}_{p}(\Omega_{k})$, the space of total degree $p$ polynomials on
$\Omega_{k}$, and $n_{p}$ is the dimension of the basis.

Ideally, the coefficients $c_{ji}$ are chosen such that the DGD basis
$\phi_{\nu_{i}}(\bm{x})$ has a value of one at the centroid
$\tilde{\bm{x}}_{\nu_{i}}$ and a value of zero at the centroids of the other
elements in $N_{k}$. In other words, for all $\nu_{i},\nu_{m}\in N_{k}$, we
have

$\phi_{\nu_{i}}(\tilde{\bm{x}}_{\nu_{m}})=\sum_{j=1}^{n_{p}}c_{ji}\mathcal{V}_{j}(\tilde{\bm{x}}_{\nu_{m}})=\delta_{im},$
(3)

where $\delta_{im}$ denotes the Kronecker delta. Equation (3) can be expressed
in matrix form as

$\tilde{\mathsf{V}}_{k}\mathsf{C}_{k}=\mathsf{I}_{k},$ (4)

where $\mathsf{I}_{k}$ is the $n_{k}\times n_{k}$ identity matrix and
$\mathsf{C}_{k}\in\mathbb{R}^{n_{p}\times n_{k}}$ holds the to-be-determined
coefficients. $\tilde{\mathsf{V}}_{k}\in\mathbb{R}^{n_{k}\times n_{p}}$ is the
generalized Vandermonde matrix:

$\tilde{\mathsf{V}}_{k}\equiv\begin{bmatrix}\mathcal{V}_{1}(\tilde{\bm{x}}_{\nu_{1}})&\mathcal{V}_{2}(\tilde{\bm{x}}_{\nu_{1}})&\ldots&\mathcal{V}_{n_{p}}(\tilde{\bm{x}}_{\nu_{1}})\\\
\mathcal{V}_{1}(\tilde{\bm{x}}_{\nu_{2}})&\mathcal{V}_{2}(\tilde{\bm{x}}_{\nu_{2}})&\ldots&\mathcal{V}_{n_{p}}(\tilde{\bm{x}}_{\nu_{2}})\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{V}_{1}(\tilde{\bm{x}}_{\nu_{n_{k}}})&\mathcal{V}_{2}(\tilde{\bm{x}}_{\nu_{n_{k}}})&\ldots&\mathcal{V}_{n_{p}}(\tilde{\bm{x}}_{\nu_{n_{k}}})\end{bmatrix}.$

In order for (4) to have a unique solution for $\mathsf{C}_{k}$, the number of
elements in the stencil must be the same as the number of basis functions,
$n_{k}=n_{p}$, and the element centroids must form an unisolvent point set for
$\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}$. Assuming the degrees of freedom are
stored at element centers (i.e., on the dual grid), the constraint
$n_{k}=n_{p}$ is straightforward to enforce for tensor-product DGD elements of
even degree $p$, but it is problematic for odd $p$ and unstructured grids. For
instance, only three elements are required in the stencil to ensure
$n_{k}=n_{p}$ when constructing a piecewise linear ($p=1$) total degree DGD
basis in two dimensions; however, on a triangular mesh, restricting $|N_{k}|$
to three would require excluding either the element $k$ itself or one of its
adjacent neighbors. As $p$ increases, the choice of which elements to exclude
from the stencil $N_{k}$ becomes increasingly arbitrary.

In this work, we favor increasing $n_{k}$ rather than excluding particular
neighbors and potentially biasing the stencil. One consequence of this is that
we sacrifice interpolation when $n_{k}>n_{p}$, since we solve (4) for
$\mathsf{C}_{k}$ in a least squares sense when there are more elements in the
stencil than basis functions:

$\mathsf{C}_{k}=(\tilde{\mathsf{V}}_{k}^{T}\tilde{\mathsf{V}}_{k})^{-1}\tilde{\mathsf{V}}_{k}^{T}.$
(5)

Thus, unlike the tensor-product scheme in Banks-DGD-paper , the present DGD
solution does not necessarily interpolate data at the element centroids;
nevertheless, the basis is still capable of representing degree $p$
polynomials exactly. To prove this (see also Li2019 ), it is sufficient to
show that the DGD basis can represent any basis function
$\mathcal{V}_{m}\in\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}$. Indeed, if we choose
$\tilde{u}_{\nu_{i}}=[\tilde{\mathsf{V}}_{k}]_{im}=\mathcal{V}_{m}(\tilde{\bm{x}}_{\nu_{i}})$
in the DGD-basis expansion (1) then, for all $\bm{x}\in\Omega_{k}$, we find

$\begin{split}\mathcal{U}_{h}(\bm{x})=\sum_{\nu\in
N_{k}}\tilde{u}_{\nu}\phi_{\nu}(\bm{x})&=\sum_{i=1}^{n_{k}}[\tilde{\mathsf{V}}_{k}]_{im}\sum_{j=1}^{n_{p}}c_{ji}\mathcal{V}_{j}(\bm{x})\\\
&=\sum_{j=1}^{n_{p}}\mathcal{V}_{j}(\bm{x})\underbrace{\sum_{i=1}^{n_{k}}[(\tilde{\mathsf{V}}_{k}^{T}\tilde{\mathsf{V}}_{k})^{-1}\tilde{\mathsf{V}}_{k}^{T}]_{ji}[\tilde{\mathsf{V}}_{k}]_{im}}_{\displaystyle=\delta_{jm}}=\mathcal{V}_{m}(\bm{x}).\end{split}$
(6)

It follows that the DGD basis can represent any polynomial of total degree $p$
on the entire domain $\Omega$, since the element $k$ is arbitrary.

One-dimensional DGD basis functions of degree $p=1$ to $p=4$ are illustrated
in Figure 1. Notice that the even degrees $p=2$ and $p=4$ produce
interpolatory basis functions, with $\phi_{\nu}(\bm{x})$ equal to one at
$\tilde{\bm{x}}_{\nu}$ and zero at the other nodes. This is because
$n_{k}=n_{p}$ for these basis functions, so (4) has a unqiue solution. In
contrast, the basis functions $\phi_{\nu}(\bm{x})$ for the odd degrees do not
equal one at $\tilde{\bm{x}}_{\nu}$ since a symmetric stencil is enforced on
the interior elements, resulting in $n_{k}>n_{p}$.

(a) $p=1$

(b) $p=2$

(c) $p=3$

(d) $p=4$

Figure 1: Discontinuous Galerkin difference basis functions, for degrees
$p=1$, $p=2$, $p=3$, and $p=4$, on a uniform one dimensional grid. The
$x$-axis labels indicate the element indices. Only the basis functions on the
left half of the domain are shown, since the right-half basis functions are
symmetric.

###  2.2 Reconstruction operator

In practice, we only need to evaluate the DGD basis at quadrature points in
order to compute the integrals that arise in the finite-element method. Thus,
in this section, we show how the basis functions can be used to construct a
linear mapping from the centroids in $N_{k}$ to quadrature points. We will
refer to this linear mapping as the prolongation matrix; in the finite-volume
literature this matrix is often called the reconstruction operator.

Let $X_{k}\equiv\\{\bm{x}_{q}\\}_{q=1}^{n_{q}}$ denote the quadrature points
adopted for the element domain $\Omega_{k}$. Then the discrete solution at
$\bm{x}_{q}\in X_{k}$ is given by

$\mathcal{U}_{h}(\bm{x}_{q})=\sum_{\nu\in
N_{k}}\phi_{\nu}(\bm{x}_{q})\tilde{u}_{\nu}=\sum_{\nu=1}^{K}(\mathsf{P}_{k})_{q\nu}\tilde{u}_{\nu},\quad\forall
q=1,2,\dots,n_{q},$ (7)

where $(\mathsf{P}_{k})_{q\nu}\equiv\phi_{\nu}(\bm{x}_{q})$ is simply the
basis function corresponding to element $\nu$ evaluated at the quadrature
point $\bm{x}_{q}$. Equation (7) can be written succinctly as

$\mathsf{P}_{k}\tilde{\bm{u}}=\bm{u}_{k},$ (8)

where $\mathsf{P}_{k}\in\mathbb{R}^{n_{q}\times K}$ is the prolongation matrix
for element $k$, $\tilde{\bm{u}}\in\mathbb{R}^{K}$ is a vector that holds DGD
coefficients associated with the centroids of the elements, and $\bm{u}_{k}$
is $\mathcal{U}_{h}$ evaluated at, or “prolonged” to, the quadrature points.

To avoid constructing the DGD basis explicitly, we express the prolongation
matrix $\mathsf{P}_{k}$ in terms of the polynomial basis
$\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}$. The desired expression can be inferred
from (7), (5), and (2):

$\mathsf{P}_{k}=\mathsf{V}_{k}\left(\tilde{\mathsf{V}}_{k}^{T}\tilde{\mathsf{V}}_{k}\right)^{-1}\tilde{\mathsf{V}}_{k}^{T}\mathsf{Z}_{k},$
(9)

where $\mathsf{V}_{k}\in\mathbb{R}^{n_{q}\times n_{p}}$ and its entries
$[\mathsf{V}_{k}]_{qj}=\mathcal{V}_{j}(\bm{x}_{q})$ are the polynomial basis
evaluated at the quadrature points $X_{k}$. Furthermore, the matrix
$\mathsf{Z}_{k}\in\mathbb{R}^{n_{k}\times K}$ maps from the global indices
$\\{1,2,\ldots,K\\}$ to the local indices $\\{1,2,\ldots,n_{k}\\}$ of $N_{k}$:

$(\mathsf{Z}_{k})_{i\nu}=\begin{cases}1,&\nu=\nu_{i}\in N_{k},\\\
0,&\text{otherwise}.\end{cases}$

###### Remark 1

The prolongation matrix can be constructed directly by seeking an
interpolation operator from the centroids of the stencil $N_{k}$ to the
quadrature points that is exact for all degree $p$ polynomials. The conditions
for such an operator are given by

$\mathsf{P}_{k}\mathsf{Z}_{k}^{T}\tilde{\mathsf{V}}_{k}=\mathsf{V}_{k}.$ (10)

This equation is underdetermined when $n_{q}\geq n_{p}$ — assuming a
unisolvent set of quadrature points — so additional conditions are necessary
to fix a unique $\mathsf{P}_{k}$. Here we seek the prolongation matrix that
has the minimum Frobenius norm and satisfies (10). It is straightforward to
show that such a matrix is given by (9).

###  2.3 DGD discretization of the linear advection equation

We conclude this section by applying the DGD method to semi-discretize the
two-dimensional linear advection equation. This exposition is intended to
further familiarize readers with the DGD finite-element method, but it will
also be used later to relate the method to summation-by-parts discetizations
and construct entropy-conservative/stable DGD schemes.

Consider the two dimensional, constant-coefficient linear advection equation
on the domain $\Omega$:

$\frac{\partial\mathcal{U}}{\partial
t}+\lambda_{x}\frac{\partial\mathcal{U}}{\partial
x}+\lambda_{y}\frac{\partial\mathcal{U}}{\partial
y}=0,\qquad\forall\;(x,y)\in\Omega,$ (11)

where $\bm{\lambda}=[\lambda_{x},\lambda_{y}]^{T}$ is the advection velocity.
In practice, the PDE (11) requires boundary conditions and an initial
condition; however, given our present focus on the spatial operators, we will
ignore the boundary and initial conditions for the time being.

The DGD semi-discretization of the linear advection equation is obtained by
following the usual weighted-residual approach. Let
$W_{h}\equiv\operatorname{span}\\{\phi_{k}\\}_{k=1}^{K}$ denote the DGD
finite-dimensional function space. Then the DGD weak formulation seeks
$\mathcal{U}_{h}\in W_{h}$ such that

$\sum_{k=1}^{K}\Bigg{[}\int_{\Omega_{k}}\mathcal{V}_{h}\frac{\partial\mathcal{U}_{h}}{\partial
t}\,d\Omega-\int_{\Omega_{k}}\lambda_{x}\frac{\partial\mathcal{V}_{h}}{\partial
x}\mathcal{U}_{h}\,d\Omega-\int_{\Omega_{k}}\lambda_{y}\frac{\partial\mathcal{V}_{h}}{\partial
y}\mathcal{U}_{h}\,d\Omega\\\
+\int_{\Gamma_{k}}\mathcal{V}_{h}\mathcal{U}_{h}\lambda_{x}n_{x}\,d\Gamma+\int_{\Gamma_{k}}\mathcal{V}_{h}\mathcal{U}_{h}\lambda_{y}n_{y}\,d\Gamma\Bigg{]}=0,$
(12)

for all $\displaystyle\mathcal{V}_{h}\in W_{h}$. The integrals are taken over
the elements and their boundaries to accommodate the discontinuous basis
functions.

As usual with finite-element methods, the bilinear forms in the weak
formulation (12) can be represented as matrices. To find this equivalent
representation, we express the trial and test functions as

$\mathcal{U}_{h}(\bm{x})=\sum_{k=1}^{K}\tilde{u}_{k}\phi_{k}(\bm{x}),\qquad\text{and}\qquad\mathcal{V}_{h}(\bm{x})=\sum_{k=1}^{K}\tilde{v}_{k}\phi_{k}(\bm{x})$

and substitute these expansions into (12). Thus, the DGD weak formulation is
equivalent to

$\tilde{\bm{v}}^{T}\tilde{\mathsf{H}}\frac{d\tilde{\bm{u}}}{dt}-\lambda_{x}\tilde{\bm{v}}^{T}\tilde{\mathsf{Q}}_{x}^{T}\tilde{\bm{u}}-\lambda_{y}\tilde{\bm{v}}^{T}\tilde{\mathsf{Q}}_{y}^{T}\tilde{\bm{u}}+\lambda_{x}\tilde{\bm{v}}^{T}\tilde{\mathsf{E}}_{x}\tilde{\bm{u}}+\lambda_{y}\tilde{\bm{v}}^{T}\tilde{\mathsf{E}}_{y}\tilde{\bm{u}}=0,$
(13)

for all $\tilde{\bm{v}}\in\mathbb{R}^{K}$, where
$\tilde{\bm{u}}\in\mathbb{R}^{K}$ is the vector of solution coefficients, and
the matrices $\tilde{\mathsf{H}}$, $\tilde{\mathsf{Q}}_{x}$, and
$\tilde{\mathsf{E}}_{x}$ are defined by

$\begin{gathered}\tilde{\mathsf{H}}_{ij}=\sum_{k=1}^{K}\int_{\Omega_{k}}\phi_{i}\phi_{j}\,d\Omega,\qquad(\tilde{\mathsf{Q}}_{x})_{ij}=\sum_{k=1}^{K}\int_{\Omega_{k}}\phi_{i}\frac{\partial\phi_{j}}{\partial
x}\,d\Omega,\\\
\text{and}\qquad(\tilde{\mathsf{E}}_{x})_{ij}=\sum_{k=1}^{K}\int_{\Gamma_{k}}\phi_{i}\phi_{j}n_{x}\,d\Gamma.\end{gathered}$
(14)

The matrices $\tilde{\mathsf{Q}}_{y}$ and $\tilde{\mathsf{E}}_{y}$ are defined
analogously to $\tilde{\mathsf{Q}}_{x}$ and $\tilde{\mathsf{E}}_{x}$,
respectively.

###### Remark 2

The DGD mass matrix $\tilde{\mathsf{H}}$ is sparse, since
$\tilde{\mathsf{H}}_{ij}$ is non-zero only if both elements $i$ and $j$ are
included in a common stencil: that is, there exists $k\in\\{1,2,\ldots,K\\}$
such that $i,j\in N_{k}$. For the same reason, the operators
$\tilde{\mathsf{Q}}_{x}$, $\tilde{\mathsf{Q}}_{y}$, $\tilde{\mathsf{E}}_{x}$,
and $\tilde{\mathsf{E}}_{y}$ are also sparse.

###### Remark 3

Unlike DG mass matrices, the DGD mass matrix is not block diagonal, which has
implications for explicit time marching methods. However, on structured grids,
Galerkin difference methods can take advantage of the tensor-product structure
to invert the mass matrix rapidly in linear time banks:GD_highorder2019 .

##  3 DGD and summation-by-parts operators

The goal of this section is two-fold. First, we will show how DGD
discretizations can be implemented using summation-by-parts (SBP) operators
and the element prolongation matrices, $\mathsf{P}_{k}$. Second, we will prove
that DGD operators are themselves dense-norm SBP operators.

###  3.1 Multidimensional SBP operators and their properties

Consider one of the elements $\Omega_{k}\in T_{h}$, and let
$X_{k}=\\{\bm{x}_{q}\\}_{q=1}^{n_{q}}$ be a set of nodes that are in the
closure of the element subdomain, $\bm{x}_{q}\in\bar{\Omega}_{k}$,
$\forall\,q=1,2,\ldots,n_{q}$. The notation for the nodes $X_{k}$ is the same
used earlier for quadrature points; as we shall see, this choice is
deliberate.

######  Definition 1 (Summation-by-parts first-derivative operator)

The matrix $\mathsf{D}_{x,k}\in\mathbb{R}^{n_{q}\times n_{q}}$ is a degree $p$
summation-by-parts operator approximating the first-derivative with respect to
$x$ at the nodes $X_{k}$ if the following conditions are satisfied.

  1. 1.

The difference operator exactly differentiates polynomials of total degree $p$
at the nodes $X_{k}$:

$\sum_{q=1}^{n_{q}}(\mathsf{D}_{x,k})_{rq}\mathcal{V}_{i}(\bm{x}_{q})=\frac{\partial\mathcal{V}_{i}}{\partial
x}(\bm{x}_{r}),\qquad\forall\,\bm{x}_{r}\in X_{k},$ (15)

and for all polynomials $\mathcal{V}_{i}$ in the basis
$\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}\subset\mathbb{P}_{p}(\Omega_{k})$.

  2. 2.

$\mathsf{D}_{x,k}=\mathsf{H}_{k}^{-1}\mathsf{Q}_{x,k}$, where $\mathsf{H}_{k}$
is a symmetric positive-definite matrix.

  3. 3.

$\mathsf{Q}_{x,k}=\mathsf{S}_{x,k}+\frac{1}{2}\mathsf{E}_{x,k}$, where
$(\mathsf{S}_{x,k})^{T}=-\mathsf{S}_{x,k}$ is a skew-symmetric matrix, and the
symmetric matrix $\mathsf{E}_{x,k}=(\mathsf{E}_{x,k})^{T}$ satisfies

$\sum_{r=1}^{n_{q}}\sum_{q=1}^{n_{q}}\mathcal{V}_{i}(\bm{x}_{r})(\mathsf{E}_{x,k})_{rq}\mathcal{V}_{j}(\bm{x}_{q})=\int_{\Gamma_{k}}\mathcal{V}_{i}\mathcal{V}_{j}n_{x}\,d\Gamma,$
(16)

for all basis polynomials
$\mathcal{V}_{i},\mathcal{V}_{j}\in\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}\subset\mathbb{P}_{p}(\Omega_{k})$,
where $n_{x}$ is the $x$ component of the outward pointing unit normal on
$\Gamma_{k}$.

An analogous definition holds for the SBP difference operator
$\mathsf{D}_{y,k}$, which approximates the first-derivative in the $y$
direction.

An important subset of SBP operators, called diagonal-norm operators, are
worth highlighting and will be used later. Diagonal-norm operators have
diagonal $\mathsf{H}_{k}$, and one can show that the diagonal entries
$\\{(\mathsf{H}_{k})_{qq}\\}_{q=1}^{n_{q}}$ and nodes $X_{k}$ constitute a
quadrature rule that is at least $2p-1$ exact hicken:quad2013 ;
hicken:mdimsbp2016 . That is, we have

$\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\mathcal{V}_{i}(\bm{x}_{q})=\int_{\Omega_{k}}\mathcal{V}_{i}\,d\Omega,$

for all
$\mathcal{V}_{i}\in\\{\mathcal{V}_{j}\\}_{j=1}^{n_{2p-1}}\subset\mathbb{P}_{2p-1}(\Omega_{k}).$
Note that the accuracy of $2p-1$ exactness is a lower bound. For the
subsequent analysis, we will assume that we are using diagonal-norm SBP
operators whose quadrature accuracy is at least $2p$ exact.

###  3.2 Implementation of DGD with SBP operators

In this section, we use diagonal-norm SBP operators to discretize the
constant-coefficient linear-advection equation, and then show the relationship
between this SBP discretization and the DGD discretization.

Suppose we have SBP operators $\mathsf{D}_{x,k}$ and $\mathsf{D}_{y,k}$ for
each element $k\in\\{1,2,\ldots,K\\}$ in the tesselation $T_{h}$ of $\Omega$.
Then, the SBP discretization of the linear-advection equation (11) is given by

$\frac{d\bm{u}_{k}}{dt}+\lambda_{x}\mathsf{D}_{x,k}\bm{u}_{k}+\lambda_{y}\mathsf{D}_{y,k}\bm{u}_{k}=\bm{0},\qquad\forall\,k=1,2,\ldots,K,$
(17)

where $\bm{u}_{k}\in\mathbb{R}^{n_{q}}$ denotes the SBP solution at the nodes
$X_{k}$ of element $k$. As with the DGD semi-discretization in Section 2.3,
this SBP discretization is for illustrative purposes only, since it lacks
imposition of boundary conditions and inter-element coupling.

Equation (17) is the strong form of the SBP discretization. The equivalent
weak form can be obtained by left multiplying by
$\bm{v}_{k}^{T}\mathsf{H}_{k}$, where $\bm{v}_{k}\in\mathbb{R}^{n_{q}}$
denotes a test function at the nodes $X_{k}$, and summing over all elements:

$\sum_{k=1}^{K}\Bigg{[}\bm{v}_{k}^{T}\mathsf{H}_{k}\frac{d\bm{u}_{k}}{dt}-\lambda_{x}\bm{v}_{k}^{T}\mathsf{Q}_{x,k}^{T}\bm{u}_{k}-\lambda_{y}\bm{v}_{k}^{T}\mathsf{Q}_{y,k}^{T}\bm{u}_{k}\\\
+\lambda_{x}\bm{v}_{k}^{T}\mathsf{E}_{x,k}\bm{u}_{k}+\lambda_{y}\bm{v}_{k}^{T}\mathsf{E}_{y,k}\bm{u}_{k}\Bigg{]}=0,$
(18)

for all $\bm{v}_{k}\in\mathbb{R}^{n_{q}}$ and all $k=1,2,\ldots,K$. To arrive
at (18), we used the identities
$\mathsf{Q}_{x,k}+\mathsf{Q}_{x,k}^{T}=\mathsf{E}_{x,k}$ and
$\mathsf{Q}_{y,k}+\mathsf{Q}_{y,k}^{T}=\mathsf{E}_{y,k}$.

Our goal is to relate the SBP discretization (18) and the DGD discretization
(13). To that end, we will need the following lemma, which provides identities
relating the SBP and DGD matrix operators.

###### Lemma 1

For each element $k\in\\{1,2,\ldots,K\\}$, let $\mathsf{D}_{x,k}$ be a
diagonal-norm SBP operator whose corresponding quadrature rule is at least
$2p$ exact. In addition, let $\mathsf{P}_{k}$ be the DGD prolongation matrix
defined in (7). Then the DGD matrices in (14) can be evaluated as

$\begin{gathered}\tilde{\mathsf{H}}_{ij}=\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\mathsf{P}_{k}\right]_{ij},\qquad(\tilde{\mathsf{Q}}_{x})_{ij}=\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{Q}_{x,k}\mathsf{P}_{k}\right]_{ij},\\\
\text{and}\qquad(\tilde{\mathsf{E}}_{x})_{ij}=\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{E}_{x,k}\mathsf{P}_{k}\right]_{ij}.\end{gathered}$
(19)

The proof of Lemma 1 can be found in A. Note that an analogous result exists
for the $y$-coordinate operators, $\tilde{\mathsf{Q}}_{y}$ and
$\tilde{\mathsf{E}}_{y}$.

We can now state and prove our first result.

###### Theorem 1

Let $\mathsf{D}_{x,k}$ and $\mathsf{D}_{y,k}$ be diagonal-norm SBP operators
for element $k\in\\{1,2,\ldots,K\\}$, and assume the quadrature associated
with these SBP operators is $2p$ exact. If we define
$\bm{u}_{k}=\mathsf{P}_{k}\tilde{\bm{u}}$ and
$\bm{v}_{k}=\mathsf{P}_{k}\tilde{\bm{v}}$ on each element, where
$\mathsf{P}_{k}$ is the DGD prolongation matrix defined in (7), the SBP
discretization (18) of the constant-coefficient linear advection is equivalent
to the DGD discretization (13).

###### Proof

The result will follow if we can show that each term in (18) is equivalent to
the corresponding term in (13). We will establish this equivalence for the
temporal term.

Substituting $\bm{u}_{k}=\mathsf{P}_{k}\tilde{\bm{u}}$ and
$\bm{v}_{k}=\mathsf{P}_{k}\tilde{\bm{v}}$ into the temporal term we find

$\sum_{k=1}^{K}\bm{v}_{k}^{T}\mathsf{H}_{k}\frac{d\bm{u}_{k}}{dt}=\sum_{k=1}^{K}\tilde{\bm{v}}^{T}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\frac{d}{dt}\left(\mathsf{P}_{k}\tilde{\bm{u}}\right)=\tilde{\bm{v}}^{T}\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\mathsf{P}_{k}\right]\frac{d\tilde{\bm{u}}}{dt}=\tilde{\bm{v}}^{T}\tilde{\mathsf{H}}\tilde{\bm{u}},$

where we used Lemma 1 in the final step. This demonstrates that the temporal
terms are equivalent under the assumptions. The result follows after applying
a similar substitution in the remaining terms and applying Lemma 1.

Theorem 1 is useful because it provides a procedure for constructing DGD
discretizations from SBP discretizations. This relationship between the
methods can be exploited for proofs of entropy conservation and stability, as
we will show in Section 4.

###  3.3 DGD operators are dense-norm SBP operators

The previous section described how DGD discretizations can be implemented
using SBP operators. Here we explore a closely related connection between
these discretizations; namely, that DGD difference operators are, themselves,
a type of SBP operator.

###### Theorem 2

Let $T_{h}=\\{\Omega_{k}\\}_{k=1}^{K}$ denote a tesselation of the domain
$\Omega$ into non-overlapping elements. Assume that a degree $2p$ exact
quadrature rule with positive weights and nodes $X_{k}$ is available on each
element $k$, and that

$\operatorname{rank}\left(\left[\mathsf{P}_{1}^{T},\mathsf{P}_{2}^{T},\ldots,\mathsf{P}_{K}^{T}\right]\right)=K,$

where $\mathsf{P}_{k}$ is the prolongation matrix from the element centers to
$X_{k}$. Then the degree $p$ DGD operator
$\tilde{\mathsf{D}}_{x}=\tilde{\mathsf{H}}^{-1}\tilde{\mathsf{Q}}_{x}$ is a
degree $p$ multidimensional SBP operator for the partial derivative in the $x$
direction over the domain $\Omega$.

###### Proof

The proof relies on Lemma 1 and, therefore, the availability of degree $p$
diagonal-norm SBP operators, $\mathsf{D}_{x,k}$, with quadrature accuracy at
least $2p$ on each element. The existence of such operators is guaranteed by
the assumption that there is a sufficiently accurate quadrature rule with
positive weights for each element hicken:mdimsbp2016 .

Let $\tilde{\bm{v}}_{i}\in\mathbb{R}^{K}$ denote an arbitrary basis polynomial
$\mathcal{V}_{i}\in\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}\subset\mathbb{P}_{p}(\Omega)$,
evaluated at the centroids of the elements in $T_{h}$, and
$\tilde{\bm{v}}_{i}^{\prime}\in\mathbb{R}^{K}$ denote its derivative,
$\partial\mathcal{V}_{i}/\partial x$, also evaluated at the centroids of the
elements. Similarly, let $\bm{v}_{i,k}$ and $\bm{v}_{i,k}^{\prime}$ denote the
basis function and its derivative, respectively, evaluated at the SBP nodes
$X_{k}$ of element $k$.

We begin with the SBP accuracy property (15), which, after multiplying both
sides of the identity by $\tilde{\mathsf{H}}$, is equivalent to

$\tilde{\mathsf{Q}}_{x}\tilde{\bm{v}}_{i}=\tilde{\mathsf{H}}\tilde{\bm{v}}_{i}^{\prime},$

Consider the left-hand side of the above equation. Substituting the expression
for $\tilde{\mathsf{Q}}_{x}$ from (19), we find

$\displaystyle\tilde{\mathsf{Q}}_{x}\tilde{\bm{v}}_{i}$
$\displaystyle=\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{Q}_{x,k}\mathsf{P}_{k}\right]\tilde{\bm{v}}_{i}=\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{Q}_{x,k}\bm{v}_{i,k}=\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\bm{v}_{i,k}^{\prime}$
$\displaystyle=\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\mathsf{P}_{k}\right]\tilde{\bm{v}}_{i}^{\prime}=\tilde{\mathsf{H}}\tilde{\bm{v}}_{i}^{\prime}.$

In the steps leading to this result, we used the fact that the generic SBP
operator $\mathsf{D}_{x,k}$ exactly differentiates polynomials of total degree
$p$; specifically, we used the equivalent statement
$\mathsf{Q}_{x,k}\bm{v}_{i,k}=\mathsf{H}_{k}\bm{v}_{i,k}^{\prime}$. We also
used the exactness of the prolongation operators when applied to the
polynomials $\tilde{\bm{v}}_{i}$ and $\tilde{\bm{v}}_{i}^{\prime}$; see Remark
1.

Next, we need to show that the mass, or norm, matrix is symmetric positive
definite. Symmetry is obvious from the definition. The mass matrix is also
positive definite; if $\tilde{\bm{u}}\in\mathbb{R}^{K}$ is an arbitrary
vector, then

$\tilde{\bm{u}}^{T}\tilde{\mathsf{H}}\tilde{\bm{u}}=\tilde{\bm{u}}^{T}\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\mathsf{P}_{k}\right]\tilde{\bm{u}}=\sum_{k=1}^{K}\bm{u}_{k}^{T}\mathsf{H}_{k}\bm{u}_{k}=\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}(\bm{u}_{k})_{q}^{2}.$

This sum is clearly non-negative, since the diagonal entries
$(\mathsf{H}_{k})_{qq}$ are strictly positive for diagonal-norm SBP operators.
Furthermore, the sum is strictly positive if $\tilde{\bm{u}}$ is nonzero. To
see this, suppose otherwise; that is, suppose
$\tilde{\bm{u}}^{T}\tilde{\mathsf{H}}\tilde{\bm{u}}=0$ for some
$\tilde{\bm{u}}\neq\bm{0}$. Then we must have
$\sum_{k=1}^{K}\bm{u}_{k}^{T}\mathsf{H}_{k}\bm{u}_{k}=0$, which is only
possible if $\bm{u}_{k}=\bm{0},\,\forall\,k$; in other words

$\begin{bmatrix}\mathsf{P}_{1}\\\ \mathsf{P}_{2}\\\ \vdots\\\
\mathsf{P}_{K}\end{bmatrix}\tilde{\bm{u}}=\begin{bmatrix}\bm{u}_{1}\\\
\bm{u}_{2}\\\ \vdots\\\ \bm{u}_{K}\end{bmatrix}=\begin{bmatrix}\bm{0}\\\
\bm{0}\\\ \vdots\\\ \bm{0}\end{bmatrix}.$

This contradicts the assumption that the matrix on the left has full rank of
$K$. Thus, we have shown that $\tilde{\mathsf{H}}$ is symmetric positive
definite.

We also need to show that the symmetric part of $\tilde{\mathsf{Q}}_{x}$ is
equal to $\frac{1}{2}\tilde{\mathsf{E}}_{x}$. This follows easily from Lemma
1, since

$\tilde{\mathsf{Q}}_{x}+\tilde{\mathsf{Q}}_{x}^{T}=\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\underbrace{(\mathsf{Q}_{x,k}+\mathsf{Q}_{x,k}^{T})}_{\mathsf{E}_{x,k}}\mathsf{P}_{k}=\tilde{\mathsf{E}}_{x}.$

Finally, we need to show that $\tilde{\mathsf{E}}_{x}$ satisfies (16). As
before, let $\tilde{\bm{v}}_{i}$ and $\tilde{\bm{v}}_{j}$ denote the basis
functions $\mathcal{V}_{i}$ and
$\mathcal{V}_{j}\in\\{\mathcal{V}_{r}\\}_{r=1}^{n_{p}}$, respectively,
evaluated at the centroids of the elements. Then we have

$\displaystyle\tilde{\bm{v}}_{i}^{T}\tilde{\mathsf{E}}_{x}\tilde{\bm{v}}_{j}=\tilde{\bm{v}}_{i}^{T}\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{E}_{x,k}\mathsf{P}_{k}\right]\tilde{\bm{v}}_{j}$
$\displaystyle=\sum_{k=1}^{K}\bm{v}_{i,k}^{T}\mathsf{E}_{x,k}\bm{v}_{j,k}$
$\displaystyle=\sum_{k=1}^{K}\int_{\Gamma_{k}}\mathcal{V}_{i}\mathcal{V}_{j}n_{x}\,d\Gamma=\int_{\Gamma}\mathcal{V}_{i}\mathcal{V}_{j}n_{x}\,d\Gamma.$

We arrived at the final line by using the fact that diagonal-norm SBP
operators satisfy (16), and the fact that the surface integrals over interior
faces cancel for polynomials of degree $p$. ∎

##  4 Entropy-stable discontinuous Galerkin difference discretizations

Having established the connection between Galerkin difference and SBP
operators, we can exploit existing entropy-stable SBP theory to construct
entropy-stable DGD schemes. To that end, this section begins with a brief
review of conservative hyperbolic systems that admit a strongly convex entropy
function. We then define semi-discrete entropy conservation and stability in
the context of generic diagonal-norm SBP methods. Subsequently, we show how
these SBP methods can be used to construct entropy-conservative/stable DGD
spatial discretizations. We conclude by describing our entropy-stable temporal
discretization.

###  4.1 Hyperbolic conservation laws with convex entropy

Consider a generic, hyperbolic conservation law in two dimensions, given by
the PDE

$\frac{\partial\bm{\mathcal{U}}}{\partial
t}+\frac{\partial\bm{\mathcal{F}}_{x}}{\partial
x}+\frac{\partial\bm{\mathcal{F}}_{y}}{\partial
y}=\bm{0},\quad\forall\bm{x}\in\Omega,\quad\forall t\in[0,T],$ (20)

where $\mathcal{U}(\bm{x},t)\in\mathbb{R}^{s}$ is the vector of $s$
conservative variables, and
$\bm{\mathcal{F}}_{x}:\mathbb{R}^{s}\rightarrow\mathbb{R}^{s}$ and
$\bm{\mathcal{F}}_{y}:\mathbb{R}^{s}\rightarrow\mathbb{R}^{s}$ are smooth flux
functions in the $x$ and $y$ coordinate directions, respectively.

We narrow our focus to conservation laws that have an associated convex
entropy function, $\mathcal{S}:\mathbb{R}^{s}\rightarrow\mathbb{R}$, that
satisfies

$\frac{\partial^{2}\mathcal{S}}{\partial\bm{\mathcal{U}}^{2}}\frac{\partial\bm{\mathcal{F}}_{*}}{\partial\bm{\mathcal{U}}}=\left[\frac{\partial^{2}\mathcal{S}}{\partial\bm{\mathcal{U}}^{2}}\frac{\partial\bm{\mathcal{F}}_{*}}{\partial\bm{\mathcal{U}}}\right]^{T},\qquad\bm{\mathcal{F}}_{*}\in\\{\bm{\mathcal{F}}_{x},\bm{\mathcal{F}}_{y}\\},$
(21)

where
$[\partial\bm{\mathcal{F}}_{*}/\partial\bm{\mathcal{U}}]_{ij}=\partial\mathcal{F}_{*,i}/\partial\mathcal{U}_{j}$
is a flux Jacobian, and
$[\partial^{2}\mathcal{S}/\partial\bm{\mathcal{U}}]_{ij}=\partial\mathcal{S}/\partial\mathcal{U}_{i}\partial\mathcal{U}_{j}$
is the positive definite Hessian of the entropy. The relations (21) imply
tadmor:2003  the existence of entropy fluxes
$\mathcal{G}_{x}:\mathbb{R}^{s}\rightarrow\mathbb{R}$ and
$\mathcal{G}_{y}:\mathbb{R}^{s}\rightarrow\mathbb{R}$, whose gradients
obey222We follow the convention that gradients are row vectors

$\frac{\partial\mathcal{G}_{*}}{\partial\bm{\mathcal{U}}}=\bm{\mathcal{W}}^{T}\frac{\partial\bm{\mathcal{F}}_{*}}{\partial\bm{\mathcal{U}}},$
(22)

where we have introduced the entropy variables,
$\bm{\mathcal{W}}\equiv[\partial\mathcal{S}/\partial\bm{\mathcal{U}}]^{T}$,
for convenience. For smooth states, (20) and (22) imply that the entropy is
conserved:

$\int_{\Omega}\bm{\mathcal{W}}^{T}\left[\frac{\partial\bm{\mathcal{U}}}{\partial
t}+\frac{\partial\bm{\mathcal{F}}_{x}}{\partial
x}+\frac{\partial\bm{\mathcal{F}}_{y}}{\partial
y}\right]\,d\Omega=\int_{\Omega}\left[\frac{\partial\mathcal{S}}{\partial
t}+\frac{\partial\mathcal{G}_{x}}{\partial
x}+\frac{\partial\mathcal{G}_{y}}{\partial y}\right]\,d\Omega=0.$ (23)

More generally, the entropy for physically relevant weak solutions to (20)
satisfies the following inequality (refer to the review tadmor:2003 , and the
references therein):

$\int_{\Omega}\left[\frac{\partial\mathcal{S}}{\partial
t}+\frac{\partial\mathcal{G}_{x}}{\partial
x}+\frac{\partial\mathcal{G}_{y}}{\partial y}\right]\,d\Omega\leq 0.$ (24)

###  4.2 Entropy-conservative and entropy-stable SBP discretizations

We are interested in discretizations that mimic (23) and (24), because these
properties imply an $L^{2}$ bound on the state Dafermos2010hyperbolic , and,
consequently, they impart a form of nonlinear stability. What we mean by mimic
will be made precise for diagonal-norm SBP discretizations in Definition 2
below.

Consider the following SBP discretization of (20) defined on each element $k$
in the tesselation $T_{h}$ of $\Omega$.

$\frac{d\bm{u}_{k}}{dt}+\bm{r}_{k}(\bm{u})=\bm{0},\qquad\forall
k=1,2,\ldots,K,$ (25)

where $\bm{u}_{k}\in\mathbb{R}^{sn_{q}}$ holds the conservative variables at
the quadrature nodes $X_{k}=\\{\bm{x}_{q}\\}_{q=1}^{n_{q}}$ on element $k$,
and
$\bm{u}=[\bm{u}_{1}^{T},\bm{u}_{2}^{T},\ldots,\bm{u}_{K}^{T}]^{T}\in\mathbb{R}^{sn_{q}K}$
is the compound vector of conservative variables over all $n_{q}K$ nodes. We
assume the SBP operators used in the spatial discretization $\bm{r}_{k}$ are
associated with a diagonal norm matrix
$\mathsf{H}_{k}\in\mathbb{R}^{n_{q}\times n_{q}}$. Note that
$\bm{r}_{k}(\bm{u})$ depends on the neighbours of element $k$ through
interface fluxes, which is why $\bm{r}_{k}$ is not a function of $\bm{u}_{k}$
alone.

In order to define discrete entropy conservation and stability for (25), we
need to be able to relate the discrete conservative and entropy variables to
one another. To this end, we will assume that the unknowns in $\bm{u}_{k}$ are
ordered such that the conservative variables at an arbitrary node
$\bm{x}_{q}\in X_{k}$ are sequential:

$\bm{u}_{q}=[\bm{u}_{k}]_{s(q-1)+1:sq},\qquad\forall\,q=1,2,\ldots,n_{q}.$

Then, it follows that the discrete entropy variables at the nodes are given by

$\bm{w}_{q}=\bm{\mathcal{W}}(\bm{u}_{q}),\qquad\forall\,q=1,2,\ldots,n_{q}.$

We will use
$\bm{w}_{k}=[\bm{w}_{1}^{T},\bm{w}_{2}^{T},\ldots,\bm{w}_{n_{q}}^{T}]^{T}$ to
denote the set of all entropy variables on element $k$.

We can now define what we mean by entropy conservative and entropy stable SBP
discretizations.

######  Definition 2 (Entropy-conservative/stable SBP discretizations)

The spatial discretization (25) is entropy-conservative with respect to the
diagonal, positive definite matrices $\mathsf{H}_{k}\in\mathbb{R}^{n_{q}\times
n_{q}}$, $k=1,2,\ldots,K$, if

$\sum_{i=k}^{K}\bm{w}_{k}^{T}\mathsf{H}_{k}\bm{r}_{k}(\bm{u})=0.$ (26)

Similarly, (25) is entropy-stable with respect to the $\mathsf{H}_{k}$ if

$\sum_{i=1}^{n}\bm{w}_{k}^{T}\mathsf{H}_{k}\bm{r}_{k}(\bm{u})\geq 0.$ (27)

###### Remark 4

Definition 2 is general enough to accommodate boundary conditions, although
constructing entropy-conservative and entropy-stable boundary conditions is
nontrivial, in general.

From Definition 2, we can infer that an entropy-conservative SBP semi-
discretization satisfies

$\sum_{i=k}^{K}\bm{w}_{k}^{T}\mathsf{H}_{k}\left[\frac{d\bm{u}_{k}}{dt}+\bm{r}_{k}(\bm{u})\right]=\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\bm{w}_{q}^{T}\frac{d\bm{u}_{q}}{dt}=\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\frac{ds_{q}}{dt}=0,$
(28)

where $s_{q}=\mathcal{S}(\bm{u}_{q})$ is the entropy based on $\bm{u}_{q}$.
Similarly, an entropy-stable SBP discretization implies the following
condition:

$\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\frac{ds_{q}}{dt}\leq
0.$ (29)

On periodic domains, (28) and (29) are discrete analogs of (23) and (24),
respectively, since the divergence of $(\mathcal{G}_{x},\mathcal{G}_{y})$ will
vanish on a periodic domain. Furthermore, since the SBP nodes $X_{k}$ and
corresponding diagonal entries in $\mathsf{H}_{k}$ define a quadrature rule
that is degree $2p-1$ exact (at least), then

$\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\frac{ds_{q}}{dt}=\int_{\Omega}\frac{\partial\mathcal{S}}{\partial
t}\,d\Omega+\text{O}(h^{2p}).$

###  4.3 Construction of entropy-conservative/stable DGD discretizations from
diagonal-norm SBP discretizations

Entropy-conservative and entropy-stable DGD discretizations can be constructed
from SBP discretizations following an approach similar to Theorem 1 for
constant-coefficient linear advection. That is, we can define the solution at
the SBP nodes by prolonging the DGD solution. However, unlike the linear
advection case, the choice of _which_ variable to prolong is important in the
context of nonlinear conservation laws. We must prolong the entropy variables,
not the conservative variables, to ensure the DGD scheme is entropy
conservative/stable.

Let $\tilde{\bm{w}}\in\mathbb{R}^{sK}$ denote the entropy variables associated
with the DGD degrees of freedom at the center of each element. As with the SBP
discretization, suppose the variables are ordered such that entropy variables
at the $k$th degree of freedom are sequential:

$\tilde{\bm{w}}_{k}=[\tilde{\bm{w}}]_{s(k-1)+1:sk},\qquad\forall\,k=1,2,\ldots,K.$

Note that the DGD vector $\tilde{\bm{w}}_{k}$ should not be confused with the
SBP vector $\bm{w}_{k}$: the former has $s$ entries while the latter has
$sn_{q}$ entries.

With the above ordering, we can prolong the entropy variables to the SBP nodes
of element $k$ using the (scalar) prolongation operator $\mathsf{P}_{k}$, and
the Kronecker product $\otimes$:

$\bm{w}_{k}=(\mathsf{P}_{k}\otimes\mathsf{I}_{s})\tilde{\bm{w}}=\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}\tilde{\bm{w}},\qquad\forall\,k=1,2,\ldots,K,$

where $\mathsf{I}_{s}$ is the $s\times s$ identity matrix, and $\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}\equiv\mathsf{P}_{k}\otimes\mathsf{I}_{s}$. Thus, for the DGD
discretization, the entropy variables at the SBP nodes are given by

$\bm{w}_{q}=\bm{w}_{q}(\tilde{\bm{w}})=\sum_{\nu=1}^{K}(\mathsf{P}_{k})_{q\nu}\tilde{\bm{w}}_{\nu},\qquad\forall\,\bm{x}_{q}\in
X_{k},\quad k=1,2,\ldots,K.$

Finally, the conservative variables at $\bm{x}_{q}$ are obtained by
transforming the entropy variables $\bm{w}_{q}$, viz.

$\bm{u}_{q}=\bm{\mathcal{U}}(\bm{w}_{q}),\qquad\forall\,q=1,2,\ldots,n_{q},$

and we will use the compound vector
$\bm{u}_{k}(\tilde{\bm{w}})=[\bm{u}_{1}^{T},\bm{u}_{2}^{T},\ldots,\bm{u}_{n_{q}}^{T}]^{T}$
to denote the conservative variables at the SBP nodes of element $k$ computed
from the prolonged entropy variables.

We construct the DGD discretization by substituting
$\bm{u}_{k}(\tilde{\bm{w}})$ into (25), left multiplying each equation by
$\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}^{T}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}$, and then summing the result over all elements:

$\sum_{k=1}^{K}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}^{T}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\left[\frac{d}{dt}\bm{u}_{k}(\tilde{\bm{w}})+\bm{r}_{k}(\tilde{\bm{w}})\right]=\bm{0},$
(30)

where $\bm{r}_{k}(\tilde{\bm{w}})$ is the SBP spatial discretization on
element $k$ evaluated using $\bm{u}_{k}(\tilde{\bm{w}})$, $k=1,2,\ldots,K$.
The discretization (30) is a system of $sK$ nonlinear ODEs in $sK$ unknowns
(i.e. the entries in $\tilde{\bm{w}}$).

###### Theorem 3

If the SBP discretization (25) is entropy-stable, then the DGD discretization
(30) satisfies

$\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\frac{ds_{q}}{dt}\leq
0,$

where $s_{q}=\mathcal{S}(\bm{w}_{q}(\tilde{\bm{w}}))$ denotes the entropy at
SBP node $q$ of element $k$ based on the prolonged DGD entropy variables. The
inequality becomes an equality if the SBP spatial discretization is entropy
conservative.

###### Proof

Suppose (25) is entropy stable. Then, if we left multiply the DGD
discretization (30) by $\tilde{\bm{w}}^{T}$ and recall that $\bm{w}_{k}=\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}\tilde{\bm{w}}$, we find

$\sum_{k=1}^{K}\tilde{\bm{w}}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern 1.5mu_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\left[\frac{d}{dt}\bm{u}_{k}(\tilde{\bm{w}})+\bm{r}_{k}(\tilde{\bm{w}})\right]\\\
=\sum_{k=1}^{K}\bm{w}_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\frac{d}{dt}\bm{u}_{k}(\tilde{\bm{w}})+\sum_{k=1}^{K}\bm{w}_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\bm{r}_{k}(\tilde{\bm{w}})=0.$ (31)

Now, the spatial term $\bm{r}_{k}(\tilde{\bm{w}})$ is evaluated using the
conservative variables $\bm{u}_{k}(\tilde{\bm{w}})$, which are themselves
evaluated using the prolonged entropy variables $\bm{w}_{k}=\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}\tilde{\bm{w}}$. Consequently, since (25) is entropy stable by
assumption, we have $\sum_{k=1}^{K}\bm{w}_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\bm{r}_{k}(\tilde{\bm{w}})\geq 0$, so (31) implies

$\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\bm{w}_{q}^{T}\frac{d\bm{u}_{q}}{dt}=\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\frac{ds_{q}}{dt}\leq
0,$

as desired. The proof for the corresponding entropy-conservative spatial
discretizations is similar.

###### Remark 5

Note that the proof fails if we prolong the conservative variables, because
the prolongation operation does not commute with the nonlinear entropy-
variable transformation:

$\bm{w}_{q}=\sum_{\nu=1}^{K}(\mathsf{P}_{k})_{q\nu}\tilde{\bm{w}}_{\nu}\neq\bm{\mathcal{W}}\left({\textstyle\sum_{\nu=1}^{K}(\mathsf{P}_{k})_{q\nu}\tilde{\bm{u}}_{\nu}}\right),$

where $\tilde{\bm{u}}_{\nu}$ denotes the conservative variables at DGD degree
of freedom $\nu$.

###### Remark 6

While the DGD discretization must prolong the entropy variables to the SBP
nodes to ensure entropy conservation or stability, this does not preclude
using the conservative variables as the unknown state. If the conservative
variables are adopted as the state, it merely introduces an additional step,
namely, conversion to the entropy variables, before prolongation. Such an
approach is used in Chan2018discretely  in the context of discontinuous
Galerkin finite-element methods.

###  4.4 Entropy-stable temporal discretization

So far we have focused exclusively on semi-discrete entropy-conservation and
stability, but it is desirable to mimic (23) and (24) in a fully discrete
sense to ensure entropy remains bounded. Therefore, for our temporal
discretization, we use an entropy-stable relaxation Runge-Kutta (RRK) as
described by Ranocha et al. Ranocha2020relaxation . This section briefly
summarizes our RRK implementation.

We adopt the implicit midpoint method as our baseline temporal discretization,
which is a second-order accurate Gauss-Legendre method. Thus, to move from
time $t^{n}$ to $t^{n+1}$ we first solve the midpoint discretization of (30):

$\sum_{k=1}^{K}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}^{T}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\left[\bm{u}_{k}(\tilde{\bm{w}}^{*})-\bm{u}_{k}(\tilde{\bm{w}}^{n})+\Delta
t\,\bm{r}_{k}(\tilde{\bm{w}}^{n+\frac{1}{2}})\right]=\bm{0},$ (32)

where $\Delta t=t^{n+1}-t^{n}$ and

$\tilde{\bm{w}}^{n+\frac{1}{2}}\equiv\frac{1}{2}\tilde{\bm{w}}^{*}+\frac{1}{2}\tilde{\bm{w}}^{n}.$

In our implementation, we solve (32) for $\tilde{\bm{w}}^{*}$ using Newton’s
method and sparse direct solves for the linear subproblems. The tolerances
used for Newton’s method will be discussed in Section 5.

In the conventional midpoint method, we would set
$\tilde{\bm{w}}^{n+1}=\tilde{\bm{w}}^{*}$. However, even though the midpoint
method is a symplectic integrator, $\tilde{\bm{w}}^{*}$ does not necessarily
ensure entropy conservation or stability. This is because the fully discrete
version of Theorem 3 does not hold automatically, since

$\sum_{k=1}^{K}(\tilde{\bm{w}}^{n+\frac{1}{2}})^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern 1.5mu_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\left[\bm{u}_{k}(\tilde{\bm{w}}^{*})-\bm{u}_{k}(\tilde{\bm{w}}^{n})\right]\neq\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}\left[s_{q}^{*}-s_{q}^{n}\right],$

where $s_{q}^{*}$ and $s_{q}^{n}$ are the nodal entropies based on
$\tilde{\bm{w}}^{*}$ and $\tilde{\bm{w}}^{n}$, respectively.

Rather than setting $\tilde{\bm{w}}^{n+1}=\tilde{\bm{w}}^{*}$, the RRK version
of the midpoint method sets

$\tilde{\bm{w}}^{n+1}_{\beta^{n}}=\tilde{\bm{w}}^{n}+\beta^{n}(\tilde{\bm{w}}^{*}-\tilde{\bm{w}}^{n}),$
(33)

where $\beta^{n}\in\mathbb{R}$ is chosen such that entropy
conservation/stability is respected discretely. To be more precise, let

$S(\tilde{\bm{w}})\equiv\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}(\mathsf{H}_{k})_{qq}s_{q}(\tilde{\bm{w}})$
(34)

denote the global (integrated) entropy based on the DGD entropy variables
$\tilde{\bm{w}}$, and let

$R(\tilde{\bm{w}})\equiv\sum_{k=1}^{K}\tilde{\bm{w}}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern 1.5mu_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\bm{r}_{k}(\tilde{\bm{w}})$

denote the change in entropy due to the spatial discretization. Then we seek
$\beta^{n}$ such that

$S(\tilde{\bm{w}}^{n+1}_{\beta^{n}})-S(\tilde{\bm{w}}^{n})+\beta^{n}\Delta
tR(\tilde{\bm{w}}^{n+\frac{1}{2}})=0,$ (35)

where $\tilde{\bm{w}}^{n+1}_{\beta^{n}}$ is defined by (33). Note that (35) is
a scalar equation in the unknown $\beta^{n}$, which we solve to a tolerance of
$10^{-13}$ using the secant method. Once $\beta^{n}$ has been determined,
$\tilde{\bm{w}}^{n+1}$ is computed using (33).

###### Remark 7

The entropy-stable DGD discretization under consideration does not appear to
be well suited to explicit RRK schemes. For example, the forward Euler scheme
would require solving

$\sum_{k=1}^{K}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}^{T}\mkern 1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\left[\bm{u}_{k}(\tilde{\bm{w}}^{*})-\bm{u}_{k}(\tilde{\bm{w}}^{n})+\Delta
t\,\bm{r}_{k}(\tilde{\bm{w}}^{n})\right]=\bm{0}$

for $\tilde{\bm{w}}^{*}$. Since $\bm{u}_{k}(\tilde{\bm{w}}^{*})$ is a
nonlinear function of $\tilde{\bm{w}}^{*}$, the solution cannot be determined
without using Newton’s method.

##  5 Numerical experiments

This section presents numerical experiments to verify the accuracy and
stability of a particular the DGD implementation. We also use these
experiments to investigate the spectra of entropy-stable DGD discretizations.

The section begins with a description of our particular DGD implementation. We
then review the Euler equations of gas dynamics, since they are the hyperbolic
conservation law targeted by the numerical experiments. The numerical results
are presented last.

###  5.1 DGD implementation including stencil construction

Our DGD and SBP discretizations are implemented in the Modular Finite Element
Methods (MFEM) mfem-library  library, and we restrict our focus to one-
dimensional intervals and two-dimensional triangular grids. The element-level
SBP discretization for $\bm{r}_{k}(\tilde{\bm{w}})$, which also appears in the
DGD discretization (30), is based on the one presented in Reference
crean:entropy2017 . However, we use the diagonal-norm SBP simplex operators
described in Hicken2020entropy  for the present results.

The discussion of the DGD basis in section 2.1, while general, omitted details
regarding the construction of the stencil, or patch,
$N_{k}=\\{\nu_{1},\nu_{2},\ldots,\nu_{n_{k}}\\}$, i.e., the elements that
influence the solution on element $k$. To construct $N_{k}$ on the triangular
grids considered in this work, we follow the approach in Li2019 , which we
briefly review here.

The element stencil $N_{k}$ is generated using an iterative process. The set
is initiated with the element $k$ itself, and then neighboring elements are
added until the stencil size is at least as large as the number of basis
functions for $\mathbb{P}_{p}(\Omega_{k})$. To be more precise, we introduce
the following recursive set definition:

$\begin{split}N_{k}^{0}&=\\{k\\},\\\ N_{k}^{j}&=N_{k}^{j-1}\cup
F(N_{k}^{j-1}),\end{split}$

where $F(N_{k}^{j-1})$ are the face-adjacent elements of the elements in
$N_{k}^{j-1}$ in the sense of graph theory. Then $N_{k}\equiv N_{k}^{j}$ for
the smallest $j$ such that the number of elements in $N_{k}^{j}$ is greater
than or equal to $n_{p}$. Figure 2 illustrates stencils that result when
degree $p=0$, $p=1$, and $p=2$ basis functions are used.

###### Remark 8

Other methods can be used to define $N_{k}$. For example, a stencil/patch may
be defined by all elements whose centers are inside the circle of radius $r$
centered at the $k$th element’s centroid.

(a) Stencil for $p=0$ (b) Stencil for $p=1$ (c) Stencil for $p=2$

Figure 2: Stencil construction

###  5.2 The Euler equations

We verify and investigate our entropy-stable DGD framework using the Euler
equations. The two-dimensional Euler equations take the form of (20), with the
conservative variables given by

$\bm{\mathcal{U}}=\left[\rho,\rho u,\rho v,e\right]^{T}$

where $\rho$ is the density, $[\rho u,\rho v]^{T}$ is the momentum per unit
volume, and $e$ is the total energy per unit volume. The Euler fluxes are

$\bm{\mathcal{F}}_{x}=\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho uv\\\
(e+p)u\end{bmatrix},\qquad\text{and}\qquad\bm{\mathcal{F}}_{y}=\begin{bmatrix}\rho
v\\\ \rho vu\\\ \rho v^{2}+p\\\ (e+p)u\end{bmatrix}.$ (36)

The pressure, $p$, is evaluated using the calorically perfect gas law as
$p=(\gamma-1)\left[e-\rho(u^{2}+v^{2})/2\right]$, with $\gamma=1.4$.

Many discretizations of the Euler equations use the fluxes (36) directly, but
entropy-conservative and entropy-stable SBP discretizations rely on two-point
entropy-conservative flux functions. In this work we use the Ismail-Roe flux
function Ismail2009affordable  in our element-level discretization
$\bm{r}_{k}$. Thus, based on this choice, the entropy function associated with
our discretization is

$\mathcal{S}=-\rho s/(\gamma-1),$

where $s=\ln{(p/\rho^{\gamma})}$ is the thermodynamic entropy.

###  5.3 Spatial accuracy verification using the steady isentropic vortex

We begin with a verification of spatial accuracy. In order to isolate the
spatial errors from the temporal errors, we favor a steady problem such as the
isentropic vortex. The vortex flow consists of circular streamlines and a
radially varying density and pressure. While simple, the isentropic vortex has
the advantage of offering an exact solution; see, for example, hicken:dual2014
.

The isentropic vortex problem is defined on a quarter annulus domain
$\Omega=\\{(r,\theta),|1\leq r\leq 3,\;0\leq\theta\leq\frac{\pi}{2}\\}$. The
mesh for the domain is generated in the manner described in Hicken2020entropy
. First, a uniform triangular mesh is defined in polar-coordinate space by
splitting $N\times N$ quadrilaterals into $2N^{2}$ triangles. Subsequently,
the triangles are mapped to physical space using a degree $p+1$ mapping for a
degree $p$ discretization. Figure 3(a) shows a sample mesh for $N=10$, and
Figure 3(b) illustrates the numerical density prolonged to SBP quadrature
points obtained on this mesh using a $p=3$ SBP entropy-stable discretization.

As in Reference Hicken2020entropy , we apply a slip-wall boundary condition
along the inner radius $r=1$. This boundary condition is imposed by evaluating
the Euler flux with the state’s velocity projected to be perpendicular to the
surface normal. For the remaining three sides of the domain, the exact
solution is imposed weakly using the Roe numerical-flux function. The steady
versions of (25) and (30) are solved using Newton’s method with absolute and
relative tolerances both set to $10^{-14}$.

(a) Example 2D mesh for N = 10

(b) density for $p=3$, $N=10$

Figure 3: Example mesh and solution for the steady-vortex problem

Figure 4 plots the error in the density from the SBP and DGD discretizations
of degree $p=1,2,3,$ and $4$ as a function of element size, $h=1/N$. Here, the
error is an approximation to the integral $L^{2}$ error:

$L^{2}\;\text{Error}\approx\sqrt{\sum_{k=1}^{K}(\bm{\rho}_{k}-\bm{\rho}_{k}^{\text{exact}})^{T}\hat{\mathsf{H}}_{k}(\bm{\rho}_{k}-\bm{\rho}_{k}^{\text{exact}})},$

where $\bm{\rho}_{k}$ is the density at the SBP nodes of element $k$ obtained
from either a degree $p$ SBP or DGD discretization,
$\bm{\rho}_{k}^{\text{exact}}$ is the exact density at the SBP nodes, and
$\hat{\mathsf{H}}_{k}$ denotes the $2p$ exact SBP norm matrix scaled by the
determinant of the mapping Jacobian.

The results in Figure 4(a) show that the errors produced by the SBP schemes
approach the optimal convergence rate of $p+1$ asymptotically. The DGD rates
of convergence are similar, although Figure 4(b) shows the errors are somewhat
super-convergent for degrees $p=1,3,4$ and sub-optimal for degree $p=2$.

(a) SBP discretization

(b) DGD discretization

Figure 4: Density solution $L^{2}$ error for SBP and DGD discretizations

Comparing the SBP and DGD schemes on the basis of Figures 4(a) and 4(b) would
be unfair, since a degree $p$ SBP scheme has significantly more degrees of
freedom than a degree $p$ DGD scheme on the same mesh. Thus, to compare these
two schemes more fairly, we plot the $L^{2}$ error versus degrees of freedom
in Figure 5. In this figure, the number of degrees of freedom is defined as
the number of elements, $K$, for the DGD discretization and $K$ times the
number of nodes per element for the SBP discretization. Figure 5 shows that,
when measured in terms of degrees of freedom, a degree $p$ DGD discretization
generally outperforms the degree $p$ SBP discretization. The one exception is
$p=1$, where the DGD error only begins to overlap the SBP error on the finest
mesh.

###### Remark 9

Error versus degrees of freedom, while better than error versus element size,
remains an imperfect means of comparing the SBP and DGD schemes. In
particular, the number of degrees of freedom does not reflect the potential
computational saving of the DGD scheme due to its better spectral radius and
conditioning; see, for example, the DGD spectra presented in the next section.

Figure 5: $L^{2}$ error versus degree-of-freedoms for the DGD and SBP
discretizations

We conclude our investigation of spatial accuracy by showing that the DGD
discretization produces superconvergent functional estimates. This study is
motivated by the interest in integral outputs, such as lift and drag, in
engineering analysis and design. To this end, we estimate the drag force on
the inner radius of $\Omega$ and compare this with the exact value of the
drag. Figure 6(a) shows that for degree $p=1$ and $p=2$, the drag errors of
the SBP scheme appear to approach a super-convergent rate of $2p+1$; for
degrees $p=3$ and $p=4$, the convergence rates approach $2p$ and $2p-1$,
respectively. Note that the drag error for $p=4$ stagnates around $10^{-13}$
due to the convergence tolerance used for Newton’s method. Similarly, Figure
6(b) shows that the DGD drag errors are close to $2p+1$ super-convergent for
degrees $p=1$ and $p=2$, while the convergence rates of the $p=3$ and $p=4$
schemes are closer to $2p$ and $2p-1$, respectively.

(a) SBP discretization

(b) DGD discretization

Figure 6: Drag error of vortex exerted on inner radius

###  5.4 Verification of entropy conservation/dissipation

In this section, we apply the DGD discretization to the two-dimensional
unsteady isentropic vortex problem to verify that total entropy is conserved,
in the case of an entropy-conservative scheme, or dissipated, in the case of
an entropy-stable scheme. On an unbounded domain, the analytical solution to
the unsteady vortex is given by

$\begin{split}\rho&=\left(1-\frac{\epsilon^{2}(\gamma-1)M^{2}}{8\pi^{2}}\exp{(f(x,y,t))}\right)^{\frac{1}{\gamma-1}},\\\
u&=\rho
k\left(1-\frac{\epsilon}{2\pi}k(y-y_{0})\exp{\left(\frac{f(x,y,t)}{2}\right)}\right),\\\
v&=\rho
k^{2}\frac{\epsilon}{2\pi}(x-x_{0})\exp{\left(\frac{f(x,y,t)}{2}\right)},\\\
e&=\frac{k^{2}p}{\gamma-1}+\frac{1}{2\rho}(u^{2}+v^{2}),\end{split}$ (37)

where $f(x,y,t)=1-((x-x_{0}-t)^{2}+(y-y_{0})^{2})$, and the vortex is
initially centered at $[x_{0},y_{0}]^{T}=[0.5,0.5]^{T}$. The Mach number is
$M=0.5$, $\epsilon=1.0$ is the vortex strength, and $k=15$ is a scaling factor
that controls the vortex speed.

We cannot simulate the vortex on an unbounded domain, so we consider the
square domain $\Omega=\\{(x,y)\;|\;x\in[0,1],y\in[0,1]\\}$ with periodic
boundary conditions applied on all four edges. The analytical solution (37)
with $t=0$ defines the initial condition, but we do not make any claims
regarding the convergence of the numerical solution to the exact solution,
given the finite domain. Despite possible solution errors, the total entropy
should still behave as theoretically predicted, so this case remains useful
for verification.

(a) Example mesh

(b) Initial value of $\mathcal{W}_{1}$

Figure 7: Mesh and initial condition for the unsteady vortex study

All simulations are run for a period of time $T=1/k$ time units, during which
the vortex travels in the positive $x$ direction and returns to its initial
position. The solution is advanced in time using the relaxation Runge-Kutta
(RRK) Ranocha2020relaxation  variant of the implicit midpoint method described
in Section 4.4. We use a constant CFL number of 10 based on the fastest
acoustic wave speed, and the system (32) is solved to a tolerance of $1e-12$
using Newton’s method. Figure 7 shows the mesh used for this study and the
initial condition corresponding to the first entropy variable.

Figure 8(a) and Figure 8(b) show the time histories of the change in total
entropy for the entropy-conservative and entropy-stable DGD schemes,
respectively. In both figures, the change in total entropy at time $t^{n+1}$
is given by

$\Delta S^{n+1}=S(\tilde{\bm{w}}^{n+1})-S(\tilde{\bm{w}}^{n}),$

where $S(\tilde{\bm{w}})$ is defined by (34). The entropy is conserved to
order $10^{-12}$ for the entropy conservative scheme, which is consistent with
the tolerances used in the Newton solver and secant method used for conserving
the entropy among time steps. Note that the piecewise-constant behavior in
Figure 8(a) is caused by the entropy changing in the last two digits only.
Finally, as predicted by the theory, the entropy change is always negative for
the entropy-stable scheme. Furthermore, less entropy is dissipated as the
solution degree $p$ increases.

(a) Entropy conservative scheme

(b) Entropy stable scheme

Figure 8: Change in entropy for the unsteady vortex problem

###  5.5 Spectra of the DGD Jacobians

This section investigates the spectra of the DGD Jacobians for the unsteady
vortex problem. Although the eigenvalues of the Jacobian are not relevant to
the nonlinear entropy-stability of the DGD discretizations, the spectra reveal
the potential efficiency of the schemes for conditionally stable time-marching
methods and iterative solvers. In addition, recent studies
Gassner2020stability ; Ranocha2020preventing  indicate that some entropy-
stable discretizations suffer from linear instabilities that raise doubts
about their accuracy for long-time simulations. The spectrum of the Jacobian
will reveal if the entropy-stable DGD methods suffer from similar linear
instabilities.

Linearizing the DGD semi-discretization (30) about the state
$\tilde{\bm{w}}_{0}$, we obtain the following linear ODE for the perturbation
$\tilde{\bm{w}}^{\prime}$:

$\frac{d\tilde{\bm{w}}^{\prime}}{dt}=\tilde{\mathsf{A}}\tilde{\bm{w}}^{\prime}+\bm{b},$

where $\tilde{\mathsf{A}}\equiv-\tilde{\mathsf{M}}^{-1}\tilde{\mathsf{J}}$,
and the global Jacobians are given by

$\tilde{\mathsf{M}}=\sum_{k=1}^{K}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern 1.5mu_{k}^{T}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{H}\mkern-1.5mu}\mkern
1.5mu_{k}\frac{\partial\bm{u}_{k}}{\partial\bm{w}_{k}}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k},\qquad\text{and}\qquad\tilde{\mathsf{J}}=\sum_{k=1}^{K}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}^{T}\frac{\partial\bm{r}_{k}}{\partial\bm{w}_{k}}\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern 1.5mu_{k}.$

The element-level Jacobians $\partial\bm{u}_{k}/\partial\bm{w}_{k}$ and
$\partial\bm{r}_{k}/\partial\bm{w}_{k}$ are evaluated at the prolonged
baseline state, $\mkern
1.5mu\overline{\mkern-1.5mu\mathsf{P}\mkern-1.5mu}\mkern
1.5mu_{k}\tilde{\bm{w}}_{0}$. The vector $\bm{b}$ is constant and is not
relevant to the spectral analysis.

Figures 9(a) and 9(b) show the eigenvalues333We compute the eigenvalues of the
generalized eigenvalue problem
$\tilde{\mathsf{J}}\bm{v}=\lambda\tilde{\mathsf{M}}\bm{v}$ of
$\tilde{\mathsf{A}}$ for the entropy-conservative and -stable schemes,
respectively, run with 1152 triangular elements. To normalize the results, the
eigenvalues are scaled by the spectral radius of the corresponding $p=1$
operators. The initial condition is adopted for the baseline state
$\tilde{\bm{w}}_{0}$, which is used to evaluate the Jacobians
$\tilde{\mathsf{M}}$ and $\tilde{\mathsf{J}}$; however, for the unsteady
vortex, we do not expect the spectra to change significantly over time,
because the exact solution is stationary under a suitable transformation.

(a)

(b)

Figure 9: Normalized eigenvalues from the entropy-conservative (upper row) and
entropy stable (lower row) DGD discretizations

It is noteworthy that the DGD spectral radius does not increase significantly
as the polynomial degree increases. This suggests that higher-order DGD
schemes can take larger time steps relative to DG-type discretizations when
using conditionally stable time-marching methods. On the other hand, the
appeal of using explicit time-marching schemes with the entropy-stable DGD
schemes is diminished by the structure of their temporal terms; see Remark 7.

Finally, Table 1 lists the maximum real part of the spectrum for each
discretization on two mesh sizes, denoted by $h$ and $h/2$. The table shows
that all of the entropy-conservative discretizations are linearly unstable,
with positive real parts on the order of $10^{-4}$. The linear instability of
the entropy-conservative schemes is not a significant concern, since some
entropy dissipation is typically desirable for optimal convergence rates. What
is more concerning is that the entropy-stable schemes also appear to be
linearly unstable, with the possible exception of a neutrally stable $p=1$
scheme.

###### Remark 10

In the context of a (nonlinearly) entropy-stable scheme, linear instability
does not necessarily imply that the solution will grow unbounded.
Nevertheless, linear instability may have implications for robustness — since
the growth can lead to a nonphysical state — and for long-time accuracy. For
the DGD discretization of the unsteady isentropic vortex, the impact of the
linear instability will take thousands of time steps to manifest itself, given
the relatively small size of the positive real parts.

Table 1: Maximum real part over all normalized eigenvalues in $\tilde{\mathsf{A}}$.  |  | degree
---|---|---
| mesh size | $p=1$ | $p=2$ | $p=3$ | $p=4$
entropy-conservative | h | 3.12e-04 | 2.76e-04 | 3.67e-04 | 3.33e-04
| h/2 | 1.09e-04 | 2.81e-04 | 2.82e-04 | 2.77e-04
entropy-stable | h | 2.6909e-14 | 1.6337e-07 | 9.0194e-08 | 1.4249e-05
| h/2 | 6.2607e-15 | 7.1420e-06 | 4.7396e-05 | 6.3008e-05

###  5.6 Shock-tube problem

We conclude the numerical experiments with a classical Riemann problem to
investigate whether or not the DGD discretization correctly predicts shock
speeds. We are interested in this study because the conservative properties of
the DGD discretization are not obvious: it has a nonlinear temporal term whose
Jacobian is non-diagonal.

The Riemann problem is similar to the classical Sod’s shock-tube problem. The
governing equations are the one dimensional compressible Euler equations and
the space-time domain is given by $x\in[0,1],\;t\in[0,0.3]$. The initial
conditions are given by

$\rho=\begin{cases}5\;&x<0.5\\\ 0.5\;&x\geq 0.5\end{cases},\quad
p=\begin{cases}1\;&x<0.5\\\ 1/10\;&x\geq 0.5\end{cases},\quad u=0.$

The initial density differs from the classical Sod problem. This change was
necessary because the present scheme does not include a shock-capturing
method, and the oscillations at the shock produce non-physical states when the
classical initial condition is used.

Figure 10 shows the density, velocity and pressure profiles at $t=0.3$ units
based on the $p=2$ DGD discretization on a uniform mesh with 200 elements.
Oscillations are present at the shock and contact discontinuity, because, as
explained above, the underlying scheme has no shock-capturing method. However,
the locations of the discontinuities agree with those of the exact solution.

For a more quantitative assessment, Table 2 lists the $L^{1}$ errors in the
first entropy variable, $\mathcal{W}_{1}=\partial\mathcal{S}/\partial\rho$,
for different mesh and degree configurations used to solve the shock-tube
problem. We see that the $L^{1}$ error decreases with refinement, which
suggests that the DGD scheme converges to the weak solution in the $L^{1}$
norm, despite the scheme’s nonlinear and non-diagonal temporal discretization.
Note that the rates of convergence for the even-order schemes are consistent
with the results in banks08a .

Table 2: $L^{1}$ errors for the Sod shock-tube problem, and estimated $L^{1}$ convergence rates.  |  | $p=1$ |  | $p=2$
---|---|---|---|---
$K$ |  |  $L^{1}$ error |  $L^{1}$ rate |  |  $L^{1}$ error |  $L^{1}$ rate
100 |  | 0.06436 | — |  | 0.03333 | —
200 |  | 0.04116 | 0.6449 |  | 0.02009 | 0.7303
400 |  | 0.02771 | 0.5708 |  | 0.01266 | 0.6662
800 |  | 0.02071 | 0.4201 |  | 0.008556 | 0.5653
|  | $p=3$ |  | $p=4$
$K$ |  |  $L^{1}$ error |  $L^{1}$ rate |  |  $L^{1}$ error |  $L^{1}$ rate
100 |  | 0.03872 | — |  | 0.02661 | —
200 |  | 0.02713 | 0.5132 |  | 0.01599 | 0.7348
400 |  | 0.01968 | 0.4632 |  | 0.008951 | 0.8370
800 |  | 0.01380 | 0.5121 |  | 0.005142 | 0.7997
Figure 10: Exact and numerical solutions of Sod’s shock problem. The DGD
solution ($p=2$) is plotted at the element centers.

##  6 Summary and Conclusions

We have presented entropy-conservative and entropy-stable
discontinuousGalerkin difference (DGD) discretizations of Euler equations on
unstructured grids. The entropy DGD discretization was constructed using three
steps.

  1. 1.

Construct high-order prolongation/reconstruction operators to map DGD degrees
of freedom from element centers to the nodes of appropriate degree diagonal-
norm SBP operators.

  2. 2.

Use the prolonged DGD entropy variables in place of the SBP entropy variables
in an entropy-conservative/stable SBP discretization.

  3. 3.

Apply the transposed prolongation operator to the SBP residual vector to map
the SBP equations to the DGD equations.

We verified the accuracy of DGD discretization by solving the steady vortex
problem. Furthermore, we demonstrated that the DGD and SBP discretizations
have comparable efficiency when the $L^{2}$ error is measured versus number of
degrees of freedom. We also showed that the DGD discretizations produced
superconvergent estimates of a drag functional.

The entropy-conservation and -stability properties were confirmed by solving
the unsteady vortex problem. The entropy was conserved up to the accuracy of
the iterative methods for the entropy-conservative scheme, and the
mathematical entropy decreased monotonically for the entropy-stable schemes.

Our investigation of the spectra of the linearized entropy-stable DGD method
revealed both positive and negative results. First, the positive result: the
spectral radius of the DGD discretization is relatively insensitive to order
of accuracy, making it well suited to conditionally stable time-marching
methods. On the other hand, some eigenvalues had positive real parts for the
high-order discretizations ($p\geq 2$), which indicates that the present
entropy-stable DGD method is not immune to the linear instability observed in
other entropy-stable SBP methods Gassner2020stability ; Ranocha2020preventing
.

Future work will include an investigation of conservation and shock capturing.
Our preliminary numerical experiments suggest that the DGD discretization is
conservative and predicts shock speeds correctly; however, a formal proof that
the scheme converges to the relevant weak solution has not been established.
In addition, shock capturing methods are needed to prevent non-physical
states.

## Acknowledgments

G. Yan was supported by the National Science Foundation under Grant No.
1554253, and S. Kaur was supported by the National Science Foundation under
Grant No. 1825991. The authors gratefully acknowledge this support. We also
thank RPI’s Scientific Computation Research Center for the use of computer
facilities.

## Declaration

### Funding

G. Yan and J. Hicken were supported by the National Science Foundation grant
number No.1554253. S. Kaur was supported by the National Science Foundation
grant number No.1825991

### Conflicts of interest/Competing interests

The authors declare that there have no conflict of interest.

### Availability of data and material

All data for the current study if available from the corresponding author upon
request.

### Code availability

The software used to generate all data for the current study is available from
the corresponding author upon request.

##  Appendix A Proof of Lemma 1

We begin with the right-hand side of the mass-matrix identity, and we recall
that $(\mathsf{P}_{k})_{qi}=\phi_{i}(\bm{x}_{q})$:

$\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{H}_{k}\mathsf{P}_{k}\right]_{ij}=\sum_{k=1}^{K}\sum_{q=1}^{n_{q}}\phi_{i}(\bm{x}_{q})(\mathsf{H}_{k})_{qq}\phi_{j}(\bm{x}_{q})=\sum_{k=1}^{K}\int_{\Omega_{k}}\phi_{i}\phi_{j}\,d\Omega=\tilde{\mathsf{H}}_{ij},$

where we used the assumption that the diagonal entries in $\mathsf{H}_{k}$ and
the nodes $X_{k}$ define a $2p$ exact quadrature for $\Omega_{k}$, and the
fact that $\phi_{i},\phi_{j}\in\mathbb{P}_{p}(\Omega_{k})$ when restricted to
element k.

For the $\tilde{\mathsf{Q}}_{x}$ identity, we will use the SBP property (15).
Specifically, for diagonal-norm SBP operators, (15) implies

$\sum_{q=1}^{n_{q}}(\mathsf{Q}_{x,k})_{rq}\phi_{i}(\bm{x}_{q})=(\mathsf{H}_{k})_{rr}\frac{\partial\phi_{i}}{\partial
x}(\bm{x}_{r}),\qquad\forall\,\bm{x}_{r}\in X_{k}.$

Note that the SBP operator differentiates $\phi_{i}$ exactly — when restricted
to element $k$ — since this function is constructed from the degree $p$
polynomial basis $\\{\mathcal{V}_{j}\\}_{j=1}^{n_{p}}$. Therefore, beginning
with the right-hand side of the $\tilde{\mathsf{Q}}_{x}$ identity, we find

$\displaystyle\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{Q}_{x,k}\mathsf{P}_{k}\right]_{ij}$
$\displaystyle=\sum_{k=1}^{K}\sum_{r=1}^{n_{q}}\sum_{q=1}^{n_{q}}\phi_{i}(\bm{x}_{r})(\mathsf{Q}_{x,k})_{rq}\phi_{j}(\bm{x}_{q})$
$\displaystyle=\sum_{k=1}^{K}\sum_{r=1}^{n_{q}}\phi_{i}(\bm{x}_{r})(\mathsf{H}_{k})_{rr}\frac{\partial\phi_{j}}{\partial
x}(\bm{x}_{r})$
$\displaystyle=\sum_{k=1}^{K}\int_{\Omega_{k}}\phi_{i}\frac{\partial\phi_{j}}{\partial
x}\,d\Omega=(\tilde{\mathsf{Q}}_{x})_{ij}.$

In the penultimate step above, we applied the diagonal-norm quadrature to the
degree $2p-1$ polynomial $\phi_{i}\frac{\partial\phi_{j}}{\partial x}$.

Finally, the $\tilde{\mathsf{E}}_{x}$ identity can be derived using the SBP
operator property (16):

$\displaystyle\left[\sum_{k=1}^{K}\mathsf{P}_{k}^{T}\mathsf{E}_{x,k}\mathsf{P}_{k}\right]_{ij}$
$\displaystyle=\sum_{k=1}^{K}\sum_{r=1}^{n_{q}}\sum_{q=1}^{n_{q}}\phi_{i}(\bm{x}_{r})(\mathsf{E}_{x,k})_{rq}\phi_{j}(\bm{x}_{q})$
$\displaystyle=\sum_{k=1}^{K}\int_{\Omega_{k}}\phi_{i}\phi_{j}n_{x}\,d\Gamma=(\tilde{\mathsf{E}}_{x})_{ij},$

where we again used the fact that the DGD basis functions are degree $p$
polynomials when restricted to an individual element. This concludes the
proof.

## References

  * (1) MFEM: Modular finite element methods library.  http://mfem.org (2019).  DOI 10.11578/dc.20171025.1248
  * (2) Banks, J., Aslam, T., Rider, W.: On sub-linear convergence for linearly degenerate waves in capturing schemes.  Journal of Computational Physics 227(14), 6985–7002 (2008).  DOI 10.1016/j.jcp.2008.04.002
  * (3) Banks, J., Hagstrom, T.: On galerkin difference methods.  Journal of Computational Physics 313 (2016).  DOI 10.1016/j.jcp.2016.02.042
  * (4) Carpenter, M.H., Fisher, T.C., Nielsen, E.J., Frankel, S.H.: Entropy Stable Spectral Collocation Schemes for the Navier–Stokes Equations: Discontinuous Interfaces.  SIAM Journal on Scientific Computing 36(5), B835–B867 (2014).  DOI 10.1137/130932193.  URL http://dx.doi.org/10.1137/130932193
  * (5) Chan, J.: On discretely entropy conservative and entropy stable discontinuous galerkin methods.  Journal of Computational Physics 362, 346–374 (2018).  DOI 10.1016/j.jcp.2018.02.033
  * (6) Crean, J., Hicken, J.E., Del Rey Fernández, D.C., Zingg, D.W., Carpenter, M.H.: Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements.  Journal of Computational Physics (2017).  (in revision)
  * (7) Dafermos, C.M.: Hyperbolic Conservation Laws in Continuum Physics, vol. 325.  Springer Berlin Heidelberg, Berlin, Heidelberg (2010).  DOI 10.1007/978-3-642-04048-1.  URL http://dx.doi.org/10.1007/978-3-642-04048-1
  * (8) Fisher, T.C.: High-order l2 stable multi-domain finite difference method for compressible flows.  Ph.D. thesis, Purdue University (2012).  URL http://libproxy.rpi.edu/login?url=https://www.proquest.com/docview/1221239207?accountid=28525.  Copyright - Database copyright ProQuest LLC; ProQuest does not claim copyright in the individual underlying works; Last updated - 2020-05-14
  * (9) Fisher, T.C., Carpenter, M.H.: High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains.  Journal of Computational Physics 252, 518–557 (2013).  DOI 10.1016/j.jcp.2013.06.014.  URL http://dx.doi.org/10.1016/j.jcp.2013.06.014
  * (10) Fisher, T.C., Carpenter, M.H., Nordström, J., Yamaleev, N.K., Swanson, C.: Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions.  Journal of Computational Physics 234, 353–375 (2013).  DOI 10.1016/j.jcp.2012.09.026.  URL http://dx.doi.org/10.1016/j.jcp.2012.09.026
  * (11) Friedrich, L., Schnücke, G., Winters, A.R., Fernández, D.C.D.R., Gassner, G.J., Carpenter, M.H.: Entropy stable space–time discontinuous Galerkin schemes with summation-by-parts property for hyperbolic conservation laws.  Journal of Scientific Computing 80(1), 175–222 (2019)
  * (12) Friedrich, L., Winters, A.R., Fernández, D.C.D.R., Gassner, G.J., Parsani, M., Carpenter, M.H.: An entropy stable h/p non-conforming discontinuous Galerkin method with the summation-by-parts property.  Journal of Scientific Computing 77(2), 689–725 (2018)
  * (13) Gassner, G.J., Svärd, M., Hindenlang, F.J.: Stability issues of entropy-stable and/or split-form high-order schemes (2020)
  * (14) Gassner, G.J., Winters, A.R., Kopriva, D.A.: A well balanced and entropy conservative discontinuous Galerkin spectral element method for the shallow water equations.  Applied Mathematics and Computation 272, 291–308 (2016)
  * (15) Hagstrom, T., Banks, J.W., Buckner, B.B., Juhnke, K.: Discontinuous galerkin difference methods for symmetric hyperbolic systems.  Journal of Scientific Computing 81(3), 1509–1526 (2019)
  * (16) Hicken, J.E.: Entropy-stable, high-order summation-by-parts discretizations without interface penalties.  Journal of Scientific Computing 82(2), 50 (2020).  DOI 10.1007/s10915-020-01154-8.  URL https://doi.org/10.1007/s10915-020-01154-8
  * (17) Hicken, J.E., Del Rey Fernández, D.C., Zingg, D.W.: Multidimensional Summation-by-Parts Operators: General Theory and Application to Simplex Elements.  SIAM Journal on Scientific Computing 38(4), A1935–A1958 (2016).  DOI 10.1137/15m1038360.  URL http://dx.doi.org/10.1137/15m1038360
  * (18) Hicken, J.E., Zingg, D.W.: Summation-by-parts operators and high-order quadrature.  Journal of Computational and Applied Mathematics 237(1), 111–125 (2013).  DOI 10.1016/j.cam.2012.07.015.  URL http://dx.doi.org/10.1016/j.cam.2012.07.015
  * (19) Hicken, J.E., Zingg, D.W.: Dual consistency and functional accuracy: a finite-difference perspective.  Journal of Computational Physics 256, 161–182 (2014).  DOI 10.1016/j.jcp.2013.08.014.  URL http://dx.doi.org/10.1016/j.jcp.2013.08.014
  * (20) Ismail, F., Roe, P.L.: Affordable, entropy-consistent Euler flux functions II: Entropy production at shocks.  Journal of Computational Physics 228(15), 5410–5436 (2009).  DOI 10.1016/j.jcp.2009.04.021.  URL http://dx.doi.org/10.1016/j.jcp.2009.04.021
  * (21) Jacangelo, J., Banks, J.W., Hagstrom, T.: Galerkin differences for high-order partial differential equations.  SIAM Journal on Scientific Computing 42(2), B447–B471 (2020).  DOI 10.1137/19M1259456.  URL https://doi.org/10.1137/19M1259456
  * (22) Li, R., Ming, P., Sun, Z., Yang, Z.: An arbitrary-order discontinuous galerkin method with one unknown per element.  Journal of Scientific Computing 80(1), 268–288 (2019).  DOI 10.1007/s10915-019-00937-y.  URL https://doi.org/10.1007/s10915-019-00937-y
  * (23) Parsani, M., Carpenter, M.H., Fisher, T.C., Nielsen, E.J.: Entropy Stable Staggered Grid Discontinuous Spectral Collocation Methods of any Order for the Compressible Navier–Stokes Equations.  SIAM Journal on Scientific Computing 38(5), A3129–A3162 (2016).  DOI 10.1137/15m1043510.  URL http://dx.doi.org/10.1137/15m1043510
  * (24) Ranocha, H., Gassner, G.J.: Preventing pressure oscillations does not fix local linear stability issues of entropy-based split-form high-order schemes (2020)
  * (25) Ranocha, H., Glaubitz, J., Öffner, P., Sonar, T.: Stability of artificial dissipation and modal filtering for flux reconstruction schemes using summation-by-parts operators.  Applied Numerical Mathematics 128, 1–23 (2018).  DOI 10.1016/j.apnum.2018.01.019.  URL http://dx.doi.org/10.1016/j.apnum.2018.01.019
  * (26) Ranocha, H., Sayyari, M., Dalcin, L., Parsani, M., Ketcheson, D.I.: Relaxation runge–kutta methods: Fully discrete explicit entropy-stable schemes for the compressible euler and navier–stokes equations.  SIAM Journal on Scientific Computing 42(2), A612–A638 (2020)
  * (27) Rojas, D., Boukharfane, R., Dalcin, L., Del Rey Fernández, D.C., Ranocha, H., Keyes, D.E., Parsani, M.: On the robustness and performance of entropy stable collocated discontinuous Galerkin methods.  Journal of Computational Physics 426, 109891 (2021).  DOI 10.1016/j.jcp.2020.109891.  URL http://www.sciencedirect.com/science/article/pii/S0021999120306653
  * (28) Shadpey, S., Zingg, D.W.: Entropy-stable multidimensional summation-by-parts discretizations on hp-adaptive curvilinear grids for hyperbolic conservation laws.  Journal of Scientific Computing 82(3), 1–46 (2020)
  * (29) Tadmor, E.: Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems.  Acta Numerica 12, 451–512 (2003).  DOI 10.1017/s0962492902000156.  URL http://dx.doi.org/10.1017/s0962492902000156