remarkRemark hypothesisHypothesis claimClaim Invariant domain preserving
approximations of hyperbolic systemsV. Carlier and F.Renac

# Invariant domain preserving high-order spectral discontinuous approximations
of hyperbolic systems

Valentin Carlier École normale supérieure, F-75230 Paris, France
<EMAIL_ADDRESS>Florent Renac DAAA, ONERA, Université Paris Saclay
F-92322 Châtillon, France<EMAIL_ADDRESS>
###### Abstract

We propose a limiting procedure to preserve invariant domains with time
explicit discrete high-order spectral discontinuous approximate solutions to
hyperbolic systems of conservation laws. Provided the scheme is discretely
conservative and satisfy geometric conservation laws at the discrete level, we
derive a condition on the time step to guaranty that the cell-averaged
approximate solution is a convex combination of states in the invariant
domain. These states are then used to define local bounds which are then
imposed to the full high-order approximate solution within the cell via an a
posteriori scaling limiter. Numerical experiments are then presented with
modal and nodal discontinuous Galerkin schemes confirm the robustness and
stability enhancement of the present approach.

###### keywords:

Hyperbolic systems, Convex invariant set, Limiting, Spectral discontinuous
method, Discontinuous Galerkin method

65M12, 65M70, 35L65

##  1 Introduction

Let $D\subset\mathbb{R}^{d}$ be an open domain with $d$ the space dimension.
We are interested here in high-order numerical solutions to hyperbolic systems
of conservations law

(1) $\left\\{\begin{array}[]{ll}\partial_{t}{\bf u}+\nabla\cdot{\bf f}({\bf
u})=0,&\text{ in }D\times(0,\infty),\\\ {\bf u}({\bf x},0)={\bf u}_{0}({\bf
x}),&\text{ in }D,\end{array}\right.$

where ${\bf u}({\bf x},t)$ represents the vector of conserved variables with
values in the set of states $\Omega^{a}\subset\mathbb{R}^{m}$ which is assumed
to be convex. The flux tensor ${\bf f}=({\bf f}_{1},\dots,{\bf
f}_{d}):\Omega^{a}\ni{\bf u}\mapsto{\bf f}({\bf u})\in\mathbb{R}^{m\times d}$
is assumed to be smooth.

Solutions to 1 may develop discontinuities in finite time even if the initial
data is smooth, therefore the equations are to be understood in the sense of
distributions. Nevertheless, in this setting we lose uniqueness of the
solution and 1 must be supplemented with further admissibility criteria. We
here focus on entropy inequalities on some twice differentiable strictly
convex function $\eta:\Omega^{a}\rightarrow\mathbb{R}$ associated with a
smooth entropy flux ${\bf q}:\Omega^{a}\rightarrow\mathbb{R}^{d}$ satisfying

(2) $\boldsymbol{\eta}^{\prime}({\bf u})^{\top}{\bf f}_{i}^{\prime}({\bf
u})={\bf q}_{i}^{\prime}({\bf u})^{\top}\quad\forall{\bf u}\in\Omega^{a},\quad
1\leq i\leq d.$

A weak solution to 1 is called an entropy weak solution if for every entropy
pair of 1 we have

(3) $\frac{\partial\eta({\bf u})}{\partial t}+\nabla\cdot{\bf q}({\bf u})\leq
0,$

in the sense of distributions. Classical (smooth) solutions respect this
condition with an equality, since one can apply the chain rule with 2. The
inequality for discontinuous solutions comes from a vanishing viscosity
argument by adding a parabolic perturbation to 1, the regularizing effect
allows to have a smooth unique solution in this case. Then, under some
structural assumptions on 1, vanishing viscosity approximations converge to an
entropy measure valued solution to 1 and 3 [9]. This result is in particular
based on the existence of convex invariant domains
$\mathcal{B}\subset\Omega^{a}$ for 1: if ${\bf u}$ is in $\mathcal{B}$, then
it remains in $\mathcal{B}$ almost everywhere in $D\times(0,\infty)$ [28, 13,
24, 40]. This property generalizes the notion of maximum principle for scalar
equations. Numerical methods keeping this property at the discrete level are
called invariant domain preserving (IDP).

We are here interested in the approximation of 1 using high-order
discontinuous spectral methods (see, e.g., [7, 14, 11, 6, 8] and references
therein) where the solution to 1 is sought under the form of discontinuous
piecewise truncated series of analytic functions over a partition of the
domain $D$. Such methods have been applied to a wide range of applications
[41, 45], and have the potential to achieve high-order accuracy efficiently on
modern parallel architectures [25, 12]. Unfortunately these approximations
suffer from spurious oscillations around discontinuities of the exact solution
due to Gibbs phenomenon [10, 17] that may cause the approximate solution to
become locally nonphysical, leading to robustness issues. A large body of
research has been proposed to address such issues with, e.g., solution and
flux limiters [31, 49, 21], entropy conservative subcell flux differencing
[11, 6, 8] artificial viscosity [20, 2], shock-capturing terms [26, 23].

We here focus on a posteriori limiters from [49, 50] scaling the cellwise
approximate solution around its cell average, thus allowing to preserve
positivity of the solution (i.e., with $\mathcal{B}=\Omega^{a}$) and maximum
principles in scalar problems, while preserving conservativity of the method.
Under some strong assumptions on the mesh, it is indeed possible to derive a
condition on the time step of the scheme so that the cell-averaged solution
remains in the invariant domains on Cartesian and simplicial grids [49, 50,
27], or on unstructured quadrangular straight-sided grids [38]. We here extend
this limiting technique to an IDP limiter for a broad class of spectral
discontinuous methods with explicit time stepping on general unstructured
meshes with possibly curved elements. Provided, the discretization method is
conservative and satisfies geometric conservation laws [43, 29] at the
discrete level, we propose a condition on the time step to guaranty that the
cell-averaged approximate solution is a convex combination of states lying in
the required invariant domains. These states are then used to define local
bounds which are then imposed to the full high-order approximate solution
within the cell via the scaling limiter. This strategy is closely related to
convex limiting [21] based on first-order IDP approximations defining local
bounds and then forcing the high-order approximation to satisfy these bounds
through flux limiting [4, 48]. This approach has been applied to finite
element approximations in [19] and discontinuous Galerkin spectral element
method in [34] among others.

The objective of this paper is hence to derive a CFL condition on the time
step and to propose an iterative algorithm for its evaluation. These are based
first on the existence of a state in the invariant region which satisfies a
trivial flux balance over each mesh cell boundaries. We call this state the
pseudo-equilibrium state and then use tricks from [36] to expand the the cell-
averaged approximate solution is a convex combination of states lying in the
invariant domains. The CFL condition is also based on the existence of a
quadrature rule to evaluate the cell-averaged solution that includes the
traces of the numerical solution used to evaluate numerical fluxes in the
scheme. This latter result generalized the work in [51] on triangular grids to
general curved polyhedral elements.

The paper is organized as follows. In Section 2 we introduce the notion of
invariant domain and invariant domain preserving Riemann solver. In Section 3
we will state and prove our theorem on the existence of a CFL for high order
schemes, present and discuss the limiting strategy. Section 4 will present
various schemes that satisfy the hypothesis of our work and state the CFL
precisely for those. Numerical experiments will be presented in Section 5, and
the conclusions follow in Section 6.

##  2 Approximate Riemann solvers

In this section we present some basic notions [5, 22] on Riemann problem,
approximate Riemann solver (ARS), and convex invariant domain that will be
used in the remainder of this work. Throughout this section, ${\bf n}$ in
$\mathbb{R}^{d}$ is a given unit vector.

###  2.1 Riemann problem and invariant domains

Let two states ${\bf u}_{L}$ and ${\bf u}_{R}$ in $\Omega^{a}$, it is
convenient for the present analysis to consider the Riemann problem in the
direction ${\bf n}$:

(4) $\partial_{t}{\bf u}+\partial_{x}{\bf f}({\bf u})\cdot{\bf n}=0,\text{ in
}\mathbb{R}\times(0,\infty),\quad{\bf u}(x,0)=\left\\{\begin{array}[]{ll}{\bf
u}_{L},&\text{if }x<0,\\\ {\bf u}_{R},&\text{if }x>0.\end{array}\right.$

We will integrate 4 over the space time slab
$[-\frac{h}{2},\frac{h}{2}]\times[0,\Delta t]$ with $h>0$ and $\Delta t>0$ the
space and time steps. We suppose here that all the Riemann problems we
consider have a self-similar entropy weak solution
$\mathcal{W}(\frac{x}{t};{\bf u}_{L},{\bf u}_{R},{\bf n})$. Let introduce the
self-similar variable $\xi=\frac{x}{\Delta t}$ and assume that there exist
$\sigma_{L}$, $\sigma_{R}$ such that: $\mathcal{W}(\xi;{\bf u}_{L},{\bf
u}_{R},{\bf n})={\bf u}_{L}$ for $\xi<\sigma_{L}$ and $\mathcal{W}(\xi;{\bf
u}_{L},{\bf u}_{R},{\bf n})={\bf u}_{R}$ for $\xi>\sigma_{R}$. We then define
the maximum wave speed in 4 by

(5) $|\lambda|({\bf u}_{L},{\bf u}_{R},{\bf
n})=\max{}(|\sigma_{L}|,|\sigma_{R}|),$

and for $\frac{\Delta t}{h}|\lambda|({\bf u}_{L},{\bf u}_{R},{\bf
n})\leq\frac{1}{2}$, we define the average over the Riemann fan [21, 22]

(6) ${\bar{\bf u}}({\bf u}_{L},{\bf u}_{R},{\bf n},\Delta
t):=\frac{1}{h}\int_{-\frac{h}{2}}^{\frac{h}{2}}\mathcal{W}\Big{(}\frac{x}{\Delta
t};{\bf u}_{L},{\bf u}_{R},{\bf n}\Big{)}dx=\frac{{\bf u}_{L}+{\bf
u}_{R}}{2}-\Delta t\big{(}{\bf f}({\bf u}_{R})-{\bf f}({\bf
u}_{L})\big{)}\cdot{\bf n}.$

Finally, we will use the definition of invariant domain from [21]: a convex
set $\mathcal{B}\subset\Omega^{a}$ is an invariant domain 1 if for all ${\bf
u}_{L}$ and ${\bf u}_{R}$ in $\mathcal{B}$, we have

${\bar{\bf u}}\big{(}{\bf u}_{L},{\bf u}_{R},{\bf n},\tfrac{\Delta
t}{h}\big{)}\in\mathcal{B}\quad\forall\;\frac{\Delta t}{h}|\lambda|({\bf
u}_{L},{\bf u}_{R},{\bf n})\leq\frac{1}{2}.$

###  2.2 Two-point numerical fluxes and approximate Riemann solvers

The discretization of 1 will rely on two-point numerical fluxes [32, 22] for
the approximation of ${\bf f}\cdot{\bf n}$ and we assume them to be consistent
and conservative:

(7) ${\bf h}({\bf u},{\bf u},{\bf n})={\bf f}({\bf u})\cdot{\bf n},\quad{\bf
h}({\bf u}_{L},{\bf u}_{R},{\bf n})=-{\bf h}({\bf u}_{R},{\bf u}_{L},-{\bf
n})\quad\forall{\bf u},{\bf u}_{L},{\bf u}_{R}\in\Omega^{a},$

and Lipschitz continuous. We also define the notion of IDP two-point flux in
the following definition.

######  Definition 2.1.

A two-point flux is said to be invariant domain preserving (IDP) for
$\mathcal{B}$ an invariant domain if we have

${\bf u}-\frac{\Delta t}{h}\big{(}{\bf h}({\bf u},{\bf u}_{R},{\bf n})-{\bf
h}({\bf u}_{L},{\bf u},{\bf n})\big{)}\in\mathcal{B}\quad\forall\;{\bf
u}_{L},{\bf u},{\bf u}_{R}\in\mathcal{B},$

under the half CFL condition

$\frac{\Delta t}{h}\max(|\lambda|({\bf u}_{L},{\bf u},{\bf n}),|\lambda|({\bf
u},{\bf u}_{R},{\bf n}))\leq\frac{1}{2}.$

We now introduce the notion of ARS and IDP ARS, which will be used to derive
IDP two-point fluxes.

######  Definition 2.2 (Approximate Riemann solver).

An ARS is a self-similar function $\mathcal{W}^{a}(\frac{x}{t};{\bf
u}_{L},{\bf u}_{R},{\bf n})$, used to approximate the solution
$\mathcal{W}(\frac{x}{t};{\bf u}_{L},{\bf u}_{R},{\bf n})$ of the Riemann
problem 4, that is consistent with the integral form of 1 [22, 16]: for any
$\frac{\Delta t}{h}|\lambda|({\bf u}_{L},{\bf u}_{R},{\bf n})\leq\frac{1}{2}$,
we have

(8)
$\frac{1}{h}\int_{-\frac{h}{2}}^{\frac{h}{2}}\mathcal{W}^{a}\Big{(}\frac{x}{\Delta
t};{\bf u}_{L},{\bf u}_{R},{\bf n}\Big{)}dx=\frac{{\bf u}_{L}+{\bf
u}_{R}}{2}-\Delta t\big{(}{\bf f}({\bf u}_{R})-{\bf f}({\bf
u}_{L})\big{)}\cdot{\bf n}.$

Using consistency in 7, and setting $\xi=\frac{x}{\Delta t}$, one defines a
two-point flux from an ARS as

(9a) $\displaystyle{\bf h}_{\mathcal{W}^{a}}({\bf u}_{L},{\bf u}_{R},{\bf n})$
$\displaystyle={\bf f}({\bf u}_{L})\cdot{\bf
n}-\int_{-\lambda}^{0}\big{(}\mathcal{W}^{a}(\xi;{\bf u}_{L},{\bf u}_{R},{\bf
n})-{\bf u}_{L}\big{)}d\xi$ (9b) $\displaystyle={\bf f}({\bf u}_{R})\cdot{\bf
n}+\int_{0}^{\lambda}\big{(}\mathcal{W}^{a}(\xi;{\bf u}_{L},{\bf u}_{R},{\bf
n})-{\bf u}_{R}\big{)}d\xi,$

where $\lambda=|\lambda|({\bf u}_{L},{\bf u}_{R},{\bf n})$. Both definitions
are equivalent due to 8. We can now define the notion of IDP ARS.

######  Definition 2.3.

An ARS $\mathcal{W}^{a}(\xi;{\bf u}_{L},{\bf u}_{R},{\bf n})$ is IDP for
$\mathcal{B}$ an invariant domain if we have

$\frac{1}{\lambda}\int_{-\lambda}^{0}\mathcal{W}^{a}(\xi;{\bf u}_{L},{\bf
u}_{R},{\bf
n})d\xi\in\mathcal{B},\quad\frac{1}{\lambda}\int_{0}^{\lambda}\mathcal{W}^{a}(\xi;{\bf
u}_{L},{\bf u}_{R},{\bf n})d\xi\in\mathcal{B}\quad\forall{\bf u}_{L},{\bf
u}_{R}\in\mathcal{B}.$

As a consequence
$\frac{1}{2\lambda}\int_{-\lambda}^{\lambda}\mathcal{W}^{a}(\xi,{\bf
u}_{L},{\bf u}_{R},{\bf n})d\xi$ is also in $\mathcal{B}$. We have the
following results linking IDP ARS and two-point numerical flux.

######  Lemma 2.4 (Interface invariant domain preservation [5]).

The ARS $\mathcal{W}^{a}$ is IDP for $\mathcal{B}$ iff. for all ${\bf
u}_{L},{\bf u}_{R}$ in $\mathcal{B}$, and $\frac{\Delta t}{h}|\lambda|({\bf
u}_{L},{\bf u}_{R},{\bf n})\leq 1$, we have

(10a) $\displaystyle{\bf u}_{L}$ $\displaystyle-\frac{\Delta t}{h}\big{(}{\bf
h}_{\mathcal{W}^{a}}({\bf u}_{L},{\bf u}_{R},{\bf n})-{\bf f}({\bf
u}_{L})\cdot{\bf n}\big{)}\in\mathcal{B},$ (10b) $\displaystyle{\bf u}_{R}$
$\displaystyle-\frac{\Delta t}{h}\big{(}{\bf f}({\bf u}_{R})\cdot{\bf n}-{\bf
h}_{\mathcal{W}^{a}}({\bf u}_{L},{\bf u}_{R},{\bf n})\big{)}\in\mathcal{B}.$

######  Proof 2.5.

Let consider 10a, a similar argument holds for 10b. From 9a, we have

${\bf u}_{L}-\frac{\Delta t}{h}\big{(}{\bf h}_{\mathcal{W}^{a}}({\bf
u}_{L},{\bf u}_{R},{\bf n})-{\bf f}({\bf u}_{L})\cdot{\bf
n}\big{)}=\big{(}1-\frac{\Delta t}{h}\lambda\big{)}{\bf u}_{L}+\frac{\Delta
t}{h}\lambda\frac{1}{\lambda}\int_{-\lambda}^{0}\big{(}\mathcal{W}^{a}(\xi,{\bf
u}_{L},{\bf u}_{R},{\bf n})d\xi,$

and since $0<\frac{\Delta t}{h}\lambda\leq 1$ this is a convex combination of
states in $\mathcal{B}$ and therefore is in $\mathcal{B}$. Conversely, taking
$\frac{\Delta t}{h}\lambda=1$ the above integrals are $\mathcal{B}$ iff. 10a
and 10b hold and we conclude from Definition 2.3.

This allows us to state the following result.

######  Lemma 2.6.

Let ${\bf h}_{\mathcal{W}^{a}}$ and ${\bf h}_{{\mathcal{W}}^{b}}$ be two-point
fluxes from two different ARS that are IDP for $\mathcal{B}$. Then, we have

(11) ${\bf u}-\frac{\Delta t}{h}\big{(}{\bf h_{\mathcal{W}^{a}}}({\bf u},{\bf
u}_{R},{\bf n})-{\bf h_{{\mathcal{W}}^{b}}}({\bf u}_{L},{\bf u},{\bf
n})\big{)}\in\mathcal{B}\quad\forall{\bf u}_{L},{\bf u},{\bf
u}_{R}\in\mathcal{B},$

under the half CFL condition

$\frac{\Delta t}{h}\max\big{(}|\lambda|({\bf u}_{L},{\bf u},{\bf
n}),|\lambda|({\bf u},{\bf u}_{R},{\bf n})\big{)}\leq\frac{1}{2}.$

######  Proof 2.7.

We rewrite 11 as

$\displaystyle{\bf u}-\frac{\Delta t}{h}\big{(}{\bf h_{\mathcal{W}^{a}}}({\bf
u},{\bf u}_{R},{\bf n})-{\bf h_{{\mathcal{W}}^{b}}}({\bf u}_{L},{\bf u},{\bf
n})\big{)}$ $\displaystyle=\frac{1}{2}\Big{(}{\bf u}-2\frac{\Delta
t}{h}\big{(}{\bf h_{\mathcal{W}^{a}}}({\bf u},{\bf u}_{R},{\bf n})-{\bf
f}({\bf u})\cdot{\bf n}\big{)}\Big{)}$ $\displaystyle+\frac{1}{2}\Big{(}{\bf
u}-2\frac{\Delta t}{h}\big{(}{\bf f}({\bf u})\cdot{\bf n}-{\bf
h_{{\mathcal{W}}^{b}}}({\bf u},{\bf u}_{R},{\bf n})\big{)}\Big{)},$

and apply Lemma 2.4 to both terms of the right-hand side with $2\frac{\Delta
t}{h}\lambda\leq 1$.

The previous lemma with ${\mathcal{W}}^{b}=\mathcal{W}^{a}$ also proves that
the three-point scheme built from an IDP ARS is also IDP [5].

##  3 Invariant domain preserving limiter

We here state and prove our main results on the existence of an explicit
condition on the time step to ensure that the cell-averaged solution from a
high-order spectral discontinuous scheme is IDP. In Section 3.1 we clarify the
schemes we are considering in this work. Our results are based on the
existence of a pseudo-equilibrium state allowing a balance of the numerical
fluxes at faces of each element which is introduced in Section 3.2 where we
prove its existence. The main result givin the condition on the time step is
given in Section 3.3. A limiting strategy based on convex bounds is described
in Section 3.4, while Section 3.5 introduces a fast algorithm to evaluate the
time step.

###  3.1 Cell-averaged fully discrete scheme

We now describe the main properties of the numerical methods we are
considering in this work. We consider here discretely conservative high-order
approximations of 1. Without loss of generality, we use an explicit forward
Euler discretization in time. High-order time integration is then performed
using strong-stability preserving Runge-Kutta methods [18] that are convex
combinations of explicit first-order schemes in time and thus keep their
stability properties. For the spatial discretization, the approximate solution
${\bf u}_{h}({\bf x},t)$ is defined locally over each element $\kappa$ of the
partition $\mathcal{T}_{h}$ of the domain $D$ in a local function cell space
${\cal V}_{h}^{p}(\kappa)$. By ${\bf u}_{h}^{(n+1)}(\cdot)={\bf
u}_{h}(\cdot,t^{(n+1)})$ we denote the solution at time
$t^{(n+1)}=t^{(n)}+\Delta t^{(n)}$ with $t^{(0)}=0$ and $\Delta t^{(n)}>0$ the
time step. The approximate solution is assumed to satisfy the following
relation for the cell-averaged solution $\langle{\bf u}_{h}\rangle_{\kappa}$:

(12) $\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa}=\langle{\bf
u}_{h}^{(n)}\rangle_{\kappa}-\Delta
t^{(n)}\sum_{k=1}^{N_{f}}s_{k}^{\kappa}{\bf h}\big{(}{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)}\quad\forall\kappa\in\mathcal{T}_{h},$

where the ${\bf x}_{k}^{\kappa}$ are some points on the faces $f$ in
$\partial\kappa$ and ${\bf u}_{h}^{\pm}({\bf
x}_{k}^{\kappa},t^{(n)})=\lim_{\epsilon\rightarrow 0^{+}}{\bf u}_{h}({\bf
x}_{k}^{\kappa}\pm\epsilon{\bf n}_{k}^{\kappa},t^{(n)})$ denote evaluations of
the traces of the solutions at ${\bf x}_{k}^{\kappa}$ (see Fig. 1). The
$s_{k}^{\kappa}>0$ are local contributions to $\frac{|f|}{|\kappa|}$ with
$|f|$ and $|\kappa|$ approximations of the face surface and element volume,
and we introduce

(13) ${\cal S}^{\kappa}:=\sum\limits_{k=1}^{N_{f}}s_{k}^{\kappa}.$

${\bf u}_{h}^{-}$${\bf u}_{h}^{+}$${\bf
n}_{k}^{\kappa}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color<EMAIL_ADDRESS>x}_{k}^{\kappa}}$${\bf
y}_{i}^{\kappa}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color@gray@fill{0.625}\bullet}$${\color[rgb]{0.625,0.625,0.625}\definecolor[named]{pgfstrokecolor}{rgb}{0.625,0.625,0.625}\pgfsys@color@gray@stroke{0.625}\pgfsys@color<EMAIL_ADDRESS>Figure 1: Notations for $d=2$ on a quadrangle: definitions of the unit outward
normal vector ${\bf n}_{k}^{\kappa}$, element quadrature nodes ${\bf
y}_{i}^{\kappa}$ (black bullets), surface quadrature node ${\bf
x}_{k}^{\kappa}$ (gray bullets), and inner and outer traces ${\bf
u}_{h}^{\pm}$ at ${\bf x}_{k}^{\kappa}$.

The geometrical quantities depend on the numerical method under consideration
and examples will be given in Section 4. By ${\bf h}$ we denote a consistent,
conservative 7, and IDP (see Definition 2.1) two-point flux. The cell-averaged
solution $\langle{\bf u}_{h}^{(n)}\rangle_{\kappa}$ is supposed to be
evaluated through a suitable quadrature rule, that includes $N_{v}$ volume
quadrature points ${\bf y}_{i}^{\kappa}$ in $\kappa$ together with the $N_{f}$
surface points ${\bf x}_{i}^{\kappa}$ on $\partial\kappa$ (see Fig. 1)
introduced in 12:

(14) $\langle{\bf
u}_{h}^{(n)}\rangle_{\kappa}=\sum\limits_{i=1}^{N_{v}}\nu_{i}^{\kappa}{\bf
u}_{h}({\bf
y}_{i}^{\kappa},t^{(n)})+\sum\limits_{i=1}^{N_{f}}\beta_{i}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{i}^{\kappa},t^{(n)}),$

where the weights $\nu_{1\leq i\leq N_{v}}^{\kappa}\geq 0$ and $\beta_{1\leq
i\leq N_{f}}^{\kappa}>0$ are assumed to satisfy

(15)
$\sum\limits_{i=1}^{N_{v}}\nu_{i}^{\kappa}+\sum\limits_{i=1}^{N_{f}}\beta_{i}^{\kappa}=1.$

Lemma 3.1 and Corollary 3.4 give a theoretical basis on existence of such a
quadrature on polygonal or polyhedral mesh elements. An explicit example is
also given in [51] in the case of triangles and in Section 4.

######  Lemma 3.1.

Let $\kappa$ be a compact subset of $\mathbb{R}^{d}$ and
$\mathcal{P}_{\kappa}$ a finite dimensional subspace of ${\cal
C}^{0}(\kappa,\mathbb{R})$, that contains the constant function $f\equiv
1_{\kappa}$. Suppose that there exists a quadrature $(\varpi_{i}^{\kappa},{\bf
y}_{i}^{\kappa})_{1\leq i\leq N_{v}}$ with positives weights
$\varpi_{i}^{\kappa}>0$ and points ${\bf y}_{i}^{\kappa}$ in $\kappa$ that
integrates exactly products of functions in $\mathcal{P}_{\kappa}$:

(16) $\int_{\kappa}f({\bf y})g({\bf
y})dV=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})g({\bf y}_{i}^{\kappa})\quad\forall
f,g\in\mathcal{P}_{\kappa}.$

Then, for given points $({\bf x}_{i}^{\kappa})_{1\leq i\leq N_{f}}$ in
$\kappa$ there exist nonnegative $(\nu_{i}^{\kappa})_{1\leq i\leq N_{v}}$ and
positive $(\beta_{i}^{\kappa})_{1\leq i\leq N_{f}}$ coefficients such that

(17) $\langle
f\rangle_{\kappa}=\sum\limits_{i=1}^{N_{v}}\nu_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})+\sum\limits_{i=1}^{N_{f}}\beta_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})\quad\forall f\in\mathcal{P}_{\kappa}.$

######  Proof 3.2.

It is known that $(f,g)\mapsto\int_{\kappa}f({\bf y})g({\bf y})dV$ defines a
scalar product on ${\cal C}^{0}(\kappa,\mathbb{R})$ and therefore on
$\mathcal{P}_{\kappa}$. Then using the Riesz representation theorem, every
linear form $\varphi:\mathcal{P}_{\kappa}\rightarrow\mathbb{R}$ can be
represented using this scalar product: there exists
$f_{\varphi}\in\mathcal{P}_{\kappa}$ such that for every
$g\in\mathcal{P}_{\kappa}$, $\varphi(g)=\int_{\kappa}f_{\varphi}({\bf
y})g({\bf y})dV$, then let
$\alpha^{\varphi}_{i}=\varpi_{i}^{\kappa}f_{\varphi}({\bf y}_{i}^{\kappa})$,
we obtain for every $g\in\mathcal{P}_{\kappa}$,
$\varphi(g)=\sum_{i=1}^{N_{v}}\alpha_{i}^{\varphi}g({\bf y}_{i}^{\kappa})$.
Now, since $f\mapsto\sum_{i=1}^{N_{f}}s_{k}^{\kappa}f({\bf x}_{i}^{\kappa})$
defines a linear form on $\mathcal{P}_{\kappa}$, for any $s_{k}^{\kappa}>0$,
it can be represented this way: there exist $(\alpha_{i}^{\kappa})_{1\leq
i\leq N_{v}}$ such that:

(18) $\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})=\sum\limits_{i=1}^{N_{v}}\alpha_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})\quad\forall f\in\mathcal{P}_{\kappa}.$

As the constant function is in $\mathcal{P}_{\kappa}$, we have for $f$ in
$\mathcal{P}_{\kappa}$, $\langle
f\rangle_{\kappa}=\sum_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})$, so for $\varepsilon_{\kappa}>0$

$\displaystyle\langle
f\rangle_{\kappa}=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})$
$\displaystyle=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})-\varepsilon_{\kappa}\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})+\varepsilon_{\kappa}\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})$
$\displaystyle=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})-\varepsilon_{\kappa}\sum\limits_{i=1}^{N_{v}}\alpha_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})+\varepsilon_{\kappa}\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})$
$\displaystyle=\sum\limits_{i=1}^{N_{v}}(\varpi_{i}^{\kappa}-\varepsilon_{\kappa}\alpha_{i}^{\kappa})f({\bf
y}_{i}^{\kappa})+\varepsilon_{\kappa}\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa}).$

Since the $\varpi_{i}^{\kappa}$ are positive, for
$\varepsilon_{\kappa}=\min_{\\{i:\alpha_{i}^{\kappa}>0\\}}(\frac{\varpi_{i}^{\kappa}}{\alpha_{i}^{\kappa}})>0$,
the $(\varpi_{i}^{\kappa}-\varepsilon_{\kappa}\alpha_{i}^{\kappa})$ are
nonnegative, then 17 holds with

(19)
$\nu_{i}^{\kappa}=\varpi_{i}^{\kappa}-\varepsilon_{\kappa}\alpha_{i}^{\kappa}\geq
0\quad\forall 1\leq i\leq
N_{v},\quad\beta_{i}^{\kappa}=\varepsilon_{\kappa}s_{i}^{\kappa}>0\quad\forall
1\leq i\leq N_{f}.$

######  Remark 3.3.

Note that $\kappa$ is usually a polyhedron and $\mathcal{P}_{\kappa}$ a
polynomial space, we can subdivide $\kappa$ into simplices and since there are
quadrature rules integrating exactly arbitrary order polynomials on simplices,
the previous lemma can be applied. Likewise, the existence of the quadrature
in Lemma 3.1 is required only for modal methods or when the DOFs are not
defined at the ${\bf x}_{k}^{\kappa}$ in 12, so the present framework also
holds for non polynomial nodal approximations with DOFs at the faces as in
[8].

The following corollary allows to explicitely define the quadrature 14 for
modal polynomial based methods.

######  Corollary 3.4.

Given a basis $(\phi_{j})_{1\leq j\leq N_{p}}$ of $\mathcal{P}_{\kappa}$ which
is orthonormal with respect to the inner product (i.e.,
$\int_{\kappa}\phi_{i}({\bf x})\phi_{j}({\bf x})dV=\delta_{ij}$ where
$\delta_{ij}$ is the Kronecker symbol), then the $\alpha_{i}^{\kappa}$ in 19
read

(20)
$\alpha_{i}^{\kappa}=\varpi_{i}^{\kappa}\sum\limits_{j=1}^{N_{p}}\sum\limits_{k=1}^{N_{f}}s_{k}^{\kappa}\phi_{j}({\bf
x}_{k}^{\kappa})\phi_{j}({\bf y}_{i}^{\kappa}).$

######  Proof 3.5.

We know that there exists $g\in\mathcal{P}_{\kappa}$ such that (see proof of
Lemma 3.1)

(21) $\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})g({\bf y}_{i}^{\kappa})\quad\forall f\in\mathcal{P}_{\kappa}.$

Expanding $g\equiv\sum\limits_{k=1}^{N_{p}}g_{k}\phi_{k}$ in the orthonormal
basis and using 21 with $f\equiv\phi_{j}$ we get

$\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}\phi_{j}({\bf
x}_{i}^{\kappa})=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}\phi_{j}({\bf
y}_{i}^{\kappa})\sum\limits_{k=1}^{N_{p}}g_{k}\phi_{k}({\bf
y}_{i}^{\kappa})=g_{j},$

for all $1\leq j\leq N_{p}$, by orthonormality of the basis. Substituting $g$
in 21 by its expansion in the basis gives

$\sum\limits_{i=1}^{N_{f}}s_{i}^{\kappa}f({\bf
x}_{i}^{\kappa})=\sum\limits_{i=1}^{N_{v}}\varpi_{i}^{\kappa}f({\bf
y}_{i}^{\kappa})\sum\limits_{j=1}^{N_{p}}\sum\limits_{k=1}^{N_{f}}s_{k}^{\kappa}\phi_{j}({\bf
x}_{k}^{\kappa})\phi_{j}({\bf y}_{i}^{\kappa})\quad\forall
f\in\mathcal{P}_{\kappa},$

and we conclude by comparing this result with 18.

Finally, the scheme is assumed to preserve uniform states in the following
sense:

(22) $\sum\limits_{k=1}^{N_{f}}s_{k}^{\kappa}{\bf n}_{k}^{\kappa}=0,$

which is a discrete version of $\tfrac{1}{|\kappa|}\oint_{\partial\kappa}{\bf
n}dS=0$ for a closed contour. Relation 22 is closely related to the discrete
geometric conservation laws [43] and is required for the numerical scheme to
preserve free-stream states [29]. In Section 5 we will present schemes that
satisfies assumptions 12, 14 and 22.

###  3.2 The pseudo-equilibrium state

We first introduce the Rusanov flux [39] :

(23) ${\bf h}_{\lambda}({\bf u}_{L},{\bf u}_{R},{\bf n})=\frac{{\bf f}({\bf
u}_{L})\cdot{\bf n}+{\bf f}({\bf u}_{R})\cdot{\bf
n}}{2}-\frac{\lambda}{2}({\bf u}_{R}-{\bf u}_{L}),$

for $\lambda\geq|\lambda|({\bf u}_{L},{\bf u}_{R},{\bf n})$ defined in 5. Note
that the Rusanov flux is derived from the following ARS:

${\cal W}_{\lambda}(\xi,{\bf u}_{L},{\bf u}_{R},{\bf
n})=\left\\{\begin{array}[]{ll}{\bf u}_{L},&\xi<-\lambda,\\\ \frac{{\bf
u}_{L}+{\bf u}_{R}}{2}-\frac{1}{2\lambda}({\bf f}({\bf u}_{R})\cdot{\bf
n}-{\bf f}({\bf u}_{L})\cdot{\bf n}),&-\lambda<\xi<\lambda,\\\ {\bf
u}_{R},&\lambda<\xi,\\\ \end{array}\right.$

so from 6 we know that it is IDP, see also [13]. We now state a result that
will allow us to rewrite 12 with updates of three-point schemes.

######  Lemma 3.6 (pseudo-equilibrium state).

Suppose that the numerical scheme satisfies 12, 14 and 22. Let $\mathcal{B}$
be a invariant domain and suppose that the internal traces ${\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})_{1\leq k\leq N_{f}}$ are in
$\mathcal{B}$, then there exists ${\bf u}^{\star}_{\kappa}={\bf
u}^{\star}_{\kappa}(t^{(n)})$ in $\mathcal{B}$ and
${\lambda}_{\kappa}^{\star}={\lambda}_{\kappa}^{\star}(t^{(n)})>0$ finite such
that

(24) $\sum_{k=1}^{N_{f}}s_{k}^{\kappa}{\bf
h}_{\lambda_{\kappa}^{\star}}\big{(}{\bf u}^{\star}_{\kappa},{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)}=0,$

with ${\bf h}_{\lambda_{\kappa}^{\star}}$ defined in 23 with
$\lambda=\lambda_{\kappa}^{\star}$ where

(25) $\lambda^{\star}_{\kappa}\geq\max_{1\leq k\leq
N_{f}}\big{(}|\lambda|({\bf u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa})\big{)},$

and the pseudo-equilibrium state is defined by

(26) ${\bf
u}^{\star}_{\kappa}=\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\bigg{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})-\frac{{\bf f}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}^{\kappa}}{\lambda^{\star}_{\kappa}}\bigg{)},\quad\tilde{\gamma}_{k}^{\kappa}:=\frac{s_{k}^{\kappa}}{{\cal
S}^{\kappa}}.$

######  Proof 3.7.

For the sake of clarity, we remove the time dependance of ${\bf u}_{h}$ since
all evaluations are done at $t^{(n)}$ and write ${\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa})$ for ${\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})$. We
first remark that from 13 and 26, we have

(27) $\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}=1.$

We introduce the following two sequences:

(28) $\left\\{\begin{array}[]{l}{\bf
u}^{\star}_{0}=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa}),\quad\lambda_{0}=\max\limits_{1\leq k\leq
N_{f}}\big{(}|\lambda|({\bf u}^{\star}_{0},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}),{\bf n}_{k}^{\kappa})\big{)}\Big{)},\\\
\lambda_{p+1}=\max\Big{(}\lambda_{p},\frac{1}{d({\bf u}^{\star}_{p},\partial
B)}+\max\limits_{1\leq k\leq N_{f}}\big{(}|\lambda|({\bf u}^{\star}_{p},{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf n}_{k}^{\kappa})\big{)}\Big{)},\quad
p\geq 0,\\\ {\bf
u}^{\star}_{p+1}=\displaystyle\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\bigg{(}\frac{{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa})+{\bf u}^{\star}_{p}}{2}-\frac{{\bf
f}\big{(}{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa})\big{)}\cdot{\bf
n}_{k}^{\kappa}-{\bf f}({\bf u}^{\star}_{p})\cdot{\bf
n}_{k}^{\kappa}}{2\lambda_{p+1}}\bigg{)},\quad p\geq 0,\end{array}\right.$

where $d({\bf u},\partial\mathcal{B})=\inf_{{\bf
v}\in\partial\mathcal{B}}\|{\bf u}-{\bf v}\|$ is the distance from ${\bf u}$
to the boundary of $\mathcal{B}$. We will show that both sequences converge
and the limits satisfy 24 and 25. We first remark that $(\lambda_{p})$ is non
decreasing and therefore converges to some $\lambda^{\star}_{\kappa}$ in
$\mathbb{R}\cup\\{+\infty\\}$ and that for all $p\geq 0$, ${\bf
u}^{\star}_{p}$ is in $\mathcal{B}$ since $\lambda_{p+1}\geq|\lambda|({\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf u}^{\star}_{p},{\bf n}_{k}^{\kappa})$ so
we can use 6 with $\Delta t=\frac{1}{2\lambda_{p+1}}$. Using successively 22
and the definition of $\tilde{\gamma}_{k}^{\kappa}$ in 26, then 27, we have

$\displaystyle{\bf u}^{\star}_{p+1}$
$\displaystyle=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\bigg{(}\frac{{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa})+{\bf u}^{\star}_{p}}{2}-\frac{{\bf f}({\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa}))\cdot{\bf
n}_{k}^{\kappa}}{2\lambda_{p+1}}\bigg{)}$ $\displaystyle=\frac{{\bf
u}^{\star}_{p}}{2}+\frac{1}{2}\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf
x}_{k}^{\kappa})-\frac{1}{2\lambda_{p+1}}\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
f}({\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa}.$

Then the first step of the sequence 28 for ${\bf u}^{\star}_{\kappa}$ reads

${\bf u}^{\star}_{1}=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf
x}_{k}^{\kappa})-\frac{1}{2\lambda_{1}}\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
f}({\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa},$

so applying the recurrence relation $p$ more times, we obtain

(29) ${\bf
u}^{\star}_{p+1}=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf
x}_{k}^{\kappa})-\sum\limits_{i=1}^{p+1}\frac{1}{2^{i}\lambda_{(p+2-i)}}\times\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
f}({\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa}.$

Let show that $\sum_{i=1}^{p}\frac{1}{2^{i}\lambda_{(p+1-i)}}$ above converges
to $\frac{1}{\lambda^{\star}_{\kappa}}$. Assume $\lambda^{\star}_{\kappa}$ is
finite and let $\epsilon>0$, since $\lambda_{p}$ converges to
$\lambda^{\star}_{\kappa}$, there exists $p_{0}$ such that
$|\frac{1}{\lambda_{p}}-\frac{1}{\lambda^{\star}_{\kappa}}|<\epsilon$ for all
$p>p_{0}$. Then for $p>p_{0}$, we set

$\displaystyle\sum\limits_{i=1}^{p}\frac{1}{2^{i}\lambda_{(p+1-i)}}-\frac{1}{\lambda^{\star}_{\kappa}}$
$\displaystyle=\sum\limits_{k=1}^{p}\frac{1}{2^{(p+1-k)}\lambda_{k}}-\frac{1}{\lambda^{\star}_{\kappa}}$
$\displaystyle=\sum\limits_{k=1}^{p_{0}}\frac{1}{2^{(p+1-k)}\lambda_{k}}+\sum\limits_{k=p_{0}+1}^{p}\frac{1}{2^{(p+1-k)}\lambda_{k}}-\frac{1}{\lambda^{\star}_{\kappa}}$
$\displaystyle=\frac{1}{2^{p+1}}\sum\limits_{k=1}^{p_{0}}\frac{2^{k}}{\lambda_{k}}+\\!\\!\\!\sum\limits_{k=p_{0}+1}^{p}\frac{2^{k}}{2^{p+1}}\left(\frac{1}{\lambda_{k}}\\!-\\!\frac{1}{\lambda^{\star}_{\kappa}}\right)-\frac{1}{\lambda^{\star}_{\kappa}}\bigg{(}1\\!-\\!\\!\\!\sum\limits_{k=p_{0}+1}^{p}\frac{2^{k}}{2^{p+1}}\bigg{)}$

and using a triangle inequality we obtain

$\displaystyle\left|\sum\limits_{i=1}^{p}\frac{1}{2^{i}\lambda_{(p+1-i)}}-\frac{1}{\lambda^{\star}_{\kappa}}\right|\leq$
$\displaystyle\frac{1}{2^{(p+1)}}\sum\limits_{k=1}^{p_{0}}\frac{2^{k}}{\lambda_{k}}+\\!\\!\sum\limits_{k=p_{0}+1}^{p}\frac{2^{k}}{2^{p+1}}\left|\frac{1}{\lambda_{k}}-\frac{1}{\lambda^{\star}_{\kappa}}\right|$
$\displaystyle+\frac{1}{\lambda^{\star}_{\kappa}}\bigg{(}1-\\!\\!\sum\limits_{k=p_{0}+1}^{p}\frac{2^{k}}{2^{p+1}}\bigg{)}$
$\displaystyle\leq$
$\displaystyle\frac{1}{2^{(p+1)}}\sum\limits_{k=1}^{p_{0}}\frac{2^{k}}{\lambda_{k}}+\sum\limits_{i=1}^{p-p_{0}}\frac{1}{2^{i}}\epsilon+\frac{1}{\lambda^{\star}_{\kappa}}\bigg{(}1-\sum\limits_{i=1}^{p-p_{0}}\frac{1}{2^{i}}\bigg{)}$
$\displaystyle=$
$\displaystyle\frac{1}{2^{p+1}}\sum\limits_{k=1}^{p_{0}}\frac{2^{k}}{\lambda_{k}}+\left(1-\frac{1}{2^{p-p_{0}}}\right)\epsilon+\frac{1}{2^{p-p_{0}}}\frac{1}{\lambda^{\star}_{\kappa}}$

which clearly converges to $\epsilon$ when $p\rightarrow\infty$, proving our
statement and 26.

Let now prove that $\lambda^{\star}_{\kappa}$ is finite by contradiction.
Suppose that $\lambda_{p}\rightarrow+\infty$, then ${\bf
u}^{\star}_{p}\rightarrow\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa})$, but the application

${\bf u}\mapsto\max_{1\leq k\leq N_{f}}(|\lambda|({\bf u},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}),{\bf n}_{k}^{\kappa}))+\frac{1}{d({\bf u},\partial B)}$

is continuous in the interior of $\mathcal{B}$ and is hence locally bounded
around $\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa})$, implying that $\lambda_{p}$ is bounded and
$\lambda^{\star}_{\kappa}$ is finite which is a contradiction. By 28, we
obviously have $\lambda_{p+1}\geq\max_{1\leq k\leq N_{f}}\big{(}|\lambda|({\bf
u}^{\star}_{p},{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf
n}_{k}^{\kappa})\big{)}$ and passing the inequality to the limit we obtain 25.

Let now prove that ${\bf u}^{\star}_{\kappa}$ is in $\mathcal{B}$, since for
all $p\geq 0$, ${\bf u}^{\star}_{p}$ is in $\mathcal{B}$, we already know that
${\bf u}^{\star}_{\kappa}$ is in the closure ${\bar{\mathcal{B}}}$. Now by
contradiction, assuming that ${\bf u}^{\star}_{\kappa}$ is not in
$\mathcal{B}$, we necessarily have ${\bf u}^{\star}_{\kappa}$ in
$\partial\mathcal{B}$, so $d({\bf
u}^{\star}_{p},\partial\mathcal{B})\rightarrow 0$ inducing
$\lambda_{p}\rightarrow+\infty$ by 28 which a contradiction. It remains to
prove 24. Using 22, we add $\frac{1}{\lambda^{\star}_{\kappa}}{\bf f}({\bf
u}^{\star}_{\kappa})\cdot\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
n}_{k}^{\kappa}=0$ to 26 and get

${\bf
u}^{\star}_{\kappa}=\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\bigg{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa})-\frac{{\bf f}({\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa}+{\bf f}({\bf
u}^{\star}_{\kappa})\cdot{\bf
n}_{k}^{\kappa}}{\lambda^{\star}_{\kappa}}\bigg{)}.$

Moving ${\bf u}^{\star}_{\kappa}$ to the right-hand side and using 27 we
obtain

$\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\bigg{(}{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa})-{\bf u}^{\star}_{\kappa}-\frac{{\bf f}({\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa}+{\bf f}({\bf
u}^{\star}_{\kappa})\cdot{\bf
n}_{k}^{\kappa}}{\lambda^{\star}_{\kappa}}\bigg{)}=0,$

and multiplying the above quantity by
$-\frac{\lambda^{\star}_{\kappa}}{2}\sum_{i=1}^{N_{f}}s_{i}^{\kappa}=-\frac{\lambda^{\star}_{\kappa}}{2}{\cal
S}^{\kappa}$ we finally get

$\sum_{k=1}^{N_{f}}s_{k}^{\kappa}\bigg{(}\frac{{\bf f}({\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}))\cdot{\bf n}_{k}^{\kappa}+{\bf f}({\bf
u}^{\star}_{\kappa})\cdot{\bf
n}_{k}^{\kappa}}{2}-\lambda^{\star}_{\kappa}\frac{({\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa})-{\bf u}^{\star}_{\kappa})}{2}\bigg{)}=0,$

which is exactly 24 from the definition of the Rusanov flux 23.

###  3.3 Invariant domain preserving schemes

Using Lemma 3.6 we now state and prove the main result of this work in the
theorem below.

######  Theorem 3.8 (Time step condition).

Assume that the numerical scheme satisfies 12, 14 and 22 and assume that
$u_{h}^{\pm}({\bf x}_{k}^{\kappa},t^{(n)})_{1\leq k\leq N_{f}}$ and
$u_{h}({\bf y}_{i}^{\kappa},t^{(n)})_{1\leq i\leq N_{v}}$ are in
$\mathcal{B}$, then under the following condition on the time step

(30) $\Delta t^{(n)}\max\limits_{\kappa\in\mathcal{T}_{h}}\max_{1\leq k\leq
N_{f}}\frac{s_{k}^{\kappa}}{\beta_{k}^{\kappa}}\max\Big{(}\lambda^{\star}_{\kappa},|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)}\Big{)}\leq\frac{1}{2},$

where $\lambda^{\star}_{\kappa}$ is defined by 25, $\langle
u_{h}^{(n+1)}\rangle_{\kappa}$ is also in $\mathcal{B}$.

######  Proof 3.9.

Once again we remove the time dependance of ${\bf u}_{h}$ for the sake of
clarity, except when explicitly needed in the evaluation of $\langle{\bf
u}_{h}^{(n+1)}\rangle_{\kappa}$. Using Lemma 3.6 we add the trivial quantity
$\Delta t^{(n)}\times$ 24 to 12 and use 14 to get

$\displaystyle\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa}=$
$\displaystyle\langle{\bf u}_{h}^{(n)}\rangle_{\kappa}-\Delta
t^{(n)}\sum_{k=1}^{N_{f}}s_{k}^{\kappa}\big{(}{\bf h}({\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}),{\bf u}_{h}^{+}({\bf x}_{k}^{\kappa}),{\bf n}_{k})-{\bf
h}_{\lambda^{\star}_{\kappa}}({\bf u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}),{\bf n}_{k}^{\kappa})\big{)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{v}}\nu_{i}{\bf u}_{h}({\bf
y}_{i}^{\kappa})+\sum_{k=1}^{N_{f}}\Big{(}\beta_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa})$ $\displaystyle-\Delta
t^{(n)}s_{k}^{\kappa}\big{(}{\bf h}({\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf
u}_{h}^{+}({\bf x}_{k}^{\kappa}),{\bf n}_{k})-{\bf
h}_{\lambda^{\star}_{\kappa}}({\bf u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa}),{\bf n}_{k}^{\kappa})\big{)}\Big{)}$ (31) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{v}}\nu_{i}{\bf u}_{h}({\bf
y}_{i}^{\kappa})+\sum_{k=1}^{N_{f}}\beta_{k}^{\kappa}{\cal U}_{k}^{\kappa,n},$

where, from Lemma 2.6 and the condition 30, the updates

(32) ${\cal U}_{k}^{\kappa,n}:={\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa})-\frac{\Delta
t^{(n)}s_{k}^{\kappa}}{\beta_{k}^{\kappa}}\Big{(}{\bf h}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf u}_{h}^{+}({\bf x}_{k}^{\kappa}),{\bf
n}_{k}^{\kappa}\big{)}-{\bf h}_{\lambda^{\star}_{\kappa}}\big{(}{\bf
u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa}),{\bf
n}_{k}^{\kappa}\big{)}\Big{)},$

are in $\mathcal{B}$. Then by 15 $\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa}$
is a convex combination of quantities in $\mathcal{B}$, which concludes the
proof.

######  Remark 3.10.

Since the set of states $\Omega^{a}$ is in general a convex invariant domain,
Theorem 3.8 can then be applied to $\mathcal{B}=\Omega^{a}$ to ensure
robustness of the scheme.

######  Remark 3.11.

In the case where we are using the quadrature on the volume defined by Lemma
3.1, from 20 and 19 and the definition
$\varepsilon_{\kappa}=\min_{\\{i:\alpha_{i}^{\kappa}>0\\}}(\frac{\varpi_{i}^{\kappa}}{\alpha_{i}^{\kappa}})$,
we have

(33)
$\frac{s_{k}^{\kappa}}{\beta_{k}^{\kappa}}=\frac{1}{\varepsilon_{\kappa}}=\max_{\\{i:\alpha_{i}^{\kappa}>0\\}}(\frac{\alpha_{i}^{\kappa}}{\varpi_{i}^{\kappa}})=\max_{1\leq
i\leq N_{v}}\sum_{j=1}^{N_{p}}\sum_{l=1}^{N_{f}}s_{l}^{\kappa}\phi_{j}({\bf
x}_{l}^{\kappa})\phi_{j}({\bf y}_{i}^{\kappa})\quad\forall 1\leq k\leq N_{f},$

and the CFL condition 30 now reads

(34) $\Delta
t^{(n)}\max\limits_{\kappa\in\mathcal{T}_{h}}\frac{1}{\varepsilon_{\kappa}}\max\Big{(}\lambda^{\star}_{\kappa},|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)}\Big{)}\leq\frac{1}{2}.$

###  3.4 Limiting strategy

Convex limiting enforces the numerical solution to preserve invariant domains
[21] through quasiconcave constraints [1]. Let recall that a function
$\psi:\mathcal{B}\rightarrow\mathbb{R}$ is quasiconcave iff. for every family
of convex coefficients $(\lambda_{i})\geq 0$, with $\sum\lambda_{i}=1$, we
have $\psi(\sum_{i}\lambda_{i}{\bf u}_{i})\geq\min_{i}\psi({\bf u}_{i})$ for
all ${\bf u}_{i}$ in $\mathcal{B}$. From Theorem 3.8 we see that for any
quasiconcave function $\psi$ we have

(35) $\psi(\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa})\geq
m_{\kappa}^{\psi}:=\min\Big{(}\psi\big{(}{\bf u}_{h}({\bf
y}_{i}^{\kappa},t^{(n)})\big{)}_{1\leq i\leq N_{v}},\psi\big{(}{\cal
U}_{k}^{\kappa,n}\big{)}_{1\leq k\leq N_{f}}\Big{)},$

where the updates ${\cal U}_{k}^{\kappa,n}$ are defined in 32. We now limit
the solution around its cell-average $\langle{\bf
u}_{h}^{(n+1)}\rangle_{\kappa}$ so that it satisfies the same bounds and we
rely on scaling limiters introduced in [50, 49] to enforce the bounds from
quasiconcave functions to points where ${\bf u}_{h}$ needs to be evaluated.
The limited solution is thus defined as

(36) $\tilde{{\bf u}}_{h}^{(n+1)}\equiv(1-\theta_{\kappa}){\bf
u}_{h}^{(n+1)}+\theta_{\kappa}\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa},$

where

$\theta_{\kappa}=\\!\\!\min\limits_{{\bf z}\in({\bf y}_{1\leq i\leq
N_{v}}^{\kappa})\cup({\bf x}_{1\leq i\leq
N_{f}}^{\kappa})}\\!\\!\max\big{\\{}0\leq t\leq 1:\psi\big{(}(1-t){\bf
u}_{h}({\bf z},t^{(n+1)})+t\langle{\bf
u}_{h}^{(n+1)}\rangle_{\kappa}\big{)}\geq m_{\kappa}^{\psi}\big{\\}}.$

This strategy may be applied to a finite family $(\psi_{i})_{1\leq i\leq
n_{c}}$ of $n_{c}$ quasiconcave functions by using the minimum value
$\theta_{\kappa}=\min\\{\theta_{\kappa}(\psi_{i}):1\leq i\leq n_{c}\\}$. The
limiter 36 is then applied locally to each cell $\kappa$ and preserves high-
order accuracy of the scheme in smooth domains [50]. The cell-average is not
modified, $\langle\tilde{\bf u}_{h}^{(n+1)}\rangle_{\kappa}=\langle{\bf
u}_{h}^{(n+1)}\rangle_{\kappa}$, and cellwise discrete conservation 12 still
holds.

###  3.5 Practical evaluation of ${\bf u}^{\star}_{\kappa}$

The evaluation of $m_{\kappa}^{\psi}$ in 35 requires to first evaluate both
$\lambda^{\star}_{\kappa}$ and ${\bf u}^{\star}_{\kappa}$ as limits of the
sequences in 28 which might be cumbersome. Here we propose another strategy
where we only need to check the convergence of $(\lambda_{p})$ to ensure that
$({\bf u}^{\star}_{p})$ has also converged to ${\bf u}^{\star}_{\kappa}$ which
is in $\mathcal{B}$. Let introduce the sequences

(37)
$\left\\{\\!\\!\\!\begin{array}[]{l}\lambda_{0}=\frac{1}{\vartheta}\max\limits_{1\leq
k\leq N_{f}}\\!|\lambda|\big{(}{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)},{\bf
u}^{\star}_{0}=\\!\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),\\\
\lambda_{p+1}=\max\Big{(}\lambda_{p},\frac{1}{\vartheta}\max\limits_{1\leq
k\leq N_{f}}\big{(}|\lambda|({\bf u}^{\star}_{p},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa})\big{)}\Big{)},\quad p\geq 0,\\\
{\bf
u}^{\star}_{p+1}=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}\Big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})-\frac{{\bf f}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}^{\kappa}}{\lambda_{p+1}}\Big{)},\quad p\geq 0,\end{array}\right.$

where

$\displaystyle\vartheta$ $\displaystyle=\max\Big{\\{}0\leq t\leq 1:{\bf
u}^{\star}_{0}-\frac{t}{\tilde{\lambda}_{1}}\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
f}\big{(}{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}^{\kappa}\in\mathcal{B}\Big{\\}},$ $\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=\max\limits_{1\leq k\leq N_{f}}\Big{(}|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)},|\lambda|\big{(}{\bf
u}^{\star}_{0},{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)}\Big{)}.$

Since ${\bf u}^{\star}_{0}\in\mathcal{B}$, we have $0<\vartheta\leq 1$ and
then ${\bf u}^{\star}_{1}\in\mathcal{B}$ by definition of $\vartheta$. Now
$\lambda_{p}$ is increasing and ${\bf u}^{\star}_{p+1}\in[{\bf
u}^{\star}_{0},{\bf u}^{\star}_{p}]\subset[{\bf u}^{\star}_{0},{\bf
u}^{\star}_{1}]$ for all $p\geq 1$ which is enough to prove convergence of
both sequences and ensure that ${\bf u}^{\star}_{p}\in\mathcal{B}$ and since
$[{\bf u}^{\star}_{0},{\bf u}^{\star}_{1}]$ is closed the limit is also in
$\mathcal{B}$ (no need to add the distance term as in 28). Obviously, the
limits satisfy 24, 25 and 26. But this time, if $\lambda_{p+1}=\lambda_{p}$
for $p\geq 1$, then ${\bf u}^{\star}_{p+1}={\bf u}^{\star}_{p}$ and from this
point both sequences are stationary. Now, the evaluation of ${\bf
u}^{\star}_{p}$ is really fast and we only need to evaluate $\max_{1\leq k\leq
N_{f}}(|\lambda|({\bf u}^{\star}_{p},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}))$.

It is possible to use a local $\lambda^{\star k}_{\kappa}$ at each vertex
${\bf x}_{\kappa}^{k}$ in $\partial\kappa$ in the above algorithm to lower the
artificial dissipation of the Rusanov flux 23 and thus avoid the induced
restriction on the time step as well as to reduce over-diffusion of the
updates ${\cal U}_{k}^{\kappa,n}$ in 32. We thus look for $\lambda^{\star
k}_{\kappa}$, $1\leq k\leq N_{f}$, and ${\bf u}^{\star}_{\kappa}$ satisfying

(38) $\sum_{k=1}^{N_{f}}s_{k}^{\kappa}{\bf h}_{\lambda^{\star
k}_{\kappa}}\big{(}{\bf u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)}=0,$

and

(39) $\lambda^{\star k}_{\kappa}\geq|\lambda|\big{(}{\bf
u}^{\star}_{\kappa},{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)}.$

We therefore introduce new sequences as

(40)
$\left\\{\begin{array}[]{ll}\lambda_{0}^{k}=\frac{1}{\vartheta}|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)},\;{\bf
u}^{\star}_{0}=\sum\limits_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),\\\
\lambda_{p+1}^{k}=\max\Big{(}\lambda_{p}^{k},\frac{1}{\vartheta}|\lambda|\big{(}{\bf
u}^{\star}_{p},{\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)}\Big{)},\;p\geq 0,\\\ {\bf
u}^{\star}_{p+1}=\sum\limits_{k=1}^{N_{f}}\frac{\tilde{\gamma}_{k}^{\kappa}\lambda_{p+1}^{k}}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda_{p+1}^{i}}\Big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})-\frac{{\bf f}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}^{\kappa}}{\lambda_{p+1}^{k}}\Big{)},\quad p\geq 0,\end{array}\right.$

where the index $k$ ranges from $1$ to $N_{f}$ and

$\displaystyle\vartheta$ $\displaystyle=\max\Big{\\{}0\leq t\leq
1:\sum\limits_{k=1}^{N_{f}}\frac{\tilde{\gamma}_{k}^{\kappa}\lambda_{1}^{k}}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda_{1}^{i}}\Big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})-t\frac{{\bf f}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}^{\kappa}}{\tilde{\lambda}_{1}^{k}}\Big{)}\in\mathcal{B}\Big{\\}},$
$\displaystyle\tilde{\lambda}_{1}^{k}$
$\displaystyle=\max\Big{(}|\lambda|\big{(}{\bf u}_{h}^{-}({\bf
x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf
n}_{k}^{\kappa}\big{)},|\lambda|\big{(}{\bf u}^{\star}_{0},{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf n}_{k}^{\kappa}\big{)}\Big{)}.$

Now all the sequences $(\lambda^{k}_{p})_{p}$ are non-decreasing and will
converge. Also for $p\geq 1$, we compute

$\Big{(}\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p+1}\Big{)}{\bf
u}^{\star}_{p+1}-\Big{(}\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p}\Big{)}{\bf
u}^{\star}_{p}=\sum_{k=1}^{N_{f}}\tilde{\gamma}_{k}^{\kappa}(\lambda_{p+1}^{k}-\lambda_{p}^{k}){\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),$

so ${\bf u}^{\star}_{p+1}$ can be recast as a convex combination

${\bf
u}^{\star}_{p+1}=\frac{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p}}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p+1}}{\bf
u}^{\star}_{p}+\sum_{k=1}^{N_{f}}\frac{\tilde{\gamma}_{k}^{\kappa}(\lambda^{k}_{p+1}-\lambda^{k}_{p})}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p+1}}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}).$

with $N_{f}+1$ positive weights such that

$\frac{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p}}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p+1}}+\sum_{k=1}^{N_{f}}\frac{\tilde{\gamma}_{k}^{\kappa}(\lambda^{k}_{p+1}-\lambda^{k}_{p})}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{i}_{p+1}}=1,$

and by recurrence ${\bf u}^{\star}_{p+1}$ is also a convex combination of
${\bf u}^{\star}_{1}$ and the ${\bf u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})$.
Therefore, if ${\bf u}^{\star}_{1}$ is in $B$, then ${\bf u}^{\star}_{p}$ is
also in $\mathcal{B}$ for all $p\geq 1$ and stays in a compact subset of
$\mathcal{B}$ which ensures that the $\lambda^{k}_{p}$ are bounded. So they
converge to some finite $\lambda^{\star k}_{\kappa}$ and ${\bf u}^{\star}_{p}$
converges to

(41) ${\bf
u}^{\star}_{\kappa}=\sum_{k=1}^{N_{f}}\frac{\tilde{\gamma}_{k}^{\kappa}\lambda^{\star
k}_{\kappa}}{\sum_{i=1}^{N_{f}}\tilde{\gamma}_{i}\lambda^{*i}_{\kappa}}\Big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})-\frac{{\bf f}\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)})\big{)}\cdot{\bf
n}_{k}}{\lambda^{\star k}_{\kappa}}\Big{)}\in\mathcal{B}.$

By definition of the $\lambda_{p}^{k}$, 39 holds and with the definition of
${\bf u}^{\star}_{\kappa}$, 38 is also satisfied, while Theorem 3.8 still
holds with $\lambda_{\kappa}^{\star k}$ instead of $\lambda_{\kappa}^{\star}$
in 30. We summarize these results in the following lemma.

######  Lemma 3.12.

Suppose that the numerical scheme satisfies 12, 14 and 22, let $\mathcal{B}$
be a invariant domain and suppose that for all $1\leq i\leq N_{f}$, ${\bf
u}_{h}^{\pm}(x_{i}^{\kappa},t^{(n)})$ are in $\mathcal{B}$, then ${\bf
u}^{\star}_{\kappa}$ defined by 41 is in $\mathcal{B}$, satisfies 38, and
there exists a family of finite and positive estimates
$({\lambda}_{\kappa}^{\star k})_{1\leq k\leq N_{f}}$ satisfying 39.

##  4 Examples of high-order spectral discontinuous methods

We here review some high-order spectral discontinuous approximations of 1
which satisfy the assumptions of discrete conservation 12, existence of
quadrature rule 14 and preservation of uniform states 22. As a consequence
there exists pseudo-equilibrium states ${\bf u}_{\kappa}^{\star}$ such that
Lemma 3.6, Lemma 3.12, and Theorem 3.8 hold, and the limiter 36 can be
applied. In the following, we consider a partition $\mathcal{T}_{h}$ of
$D\subset\mathbb{R}^{d}$, composed of non-overlapping and non-empty elements
$\kappa$, and by ${\cal F}_{h}$ we denote the set of faces in the partition.
The approximate solution is sought under the form

(42) ${\bf u}_{h}({\bf x},t)=\sum_{k=1}^{N_{p}}\phi^{\kappa}_{k}({\bf x}){\bf
U}^{\kappa}_{k}(t)\quad\forall{\bf
x}\in\kappa,\,\kappa\in\mathcal{T}_{h},\,\forall t\geq 0,$

where the basis functions $\phi^{\kappa}_{k}$ span the function space ${\cal
V}_{h}^{p}(\kappa)$ restricted onto $\kappa$ and the discrete scheme may be
written as

(43) $M_{k}^{\kappa}\frac{{\bf U}_{k,n+1}^{\kappa}-{\bf
U}_{k,n}^{\kappa}}{\Delta t^{(n)}}+{\bf R}_{k}^{\kappa}({\bf
u}_{h}^{(n)})=0\quad\forall\kappa\in\mathcal{T}_{h},\;1\leq k\leq
N_{p},\;n\geq 0,$

where ${\bf U}_{k,n}^{\kappa}={\bf U}_{k}^{\kappa}(t^{(n)})$ and the
$M_{k}^{\kappa}$ are the entries of the mass matrix.

###  4.1 Discontinuous Galerkin Method

The first numerical scheme we describe here is the discontinuous Galerkin
method with modal basis [37, 3, 7]. We look for approximate solutions in the
function space of discontinuous polynomials

${\cal V}_{h}^{p}=\bigoplus\limits_{\kappa\in\mathcal{T}_{h}}{\cal
V}_{h}^{p}(\kappa)=\\{\phi\in L^{2}(D):\;\phi|_{\kappa}\circ{\bf
x}_{\kappa}\in\mathcal{P}^{p}(\hat{K})\;\forall\kappa\in\mathcal{T}_{h}\\},$

where $\mathcal{P}^{p}(\hat{K})$ is a polynomial space over a reference
element $\hat{K}$, Each physical element $\kappa$ is the image of $\hat{K}$
through the mapping ${\bf x}={\bf x}_{\kappa}({\boldsymbol{\xi}})$ with
${\boldsymbol{\xi}}=(\xi_{1},\dots,\xi_{d})$. Likewise, each face $f$ in
${\cal F}_{h}$ is the image of a reference face $\hat{F}$ through the mapping
${\bf x}={\bf x}_{f}(\xi_{1},\dots,\xi_{d-1})$. We suppose that we have a
quadrature $({\boldsymbol{\xi}}_{i}^{V},\omega_{i}^{V})_{1\leq i\leq N_{v}}$
on $\hat{K}$ and denote ${\bf y}_{i}^{\kappa}={\bf
x}_{\kappa}({\boldsymbol{\xi}}_{i})$ (see Fig. 1). Similarly, we consider a
quadrature $({\boldsymbol{\xi}}_{k}^{f},\omega_{k}^{f})_{1\leq k\leq n_{f}}$
on $\hat{F}$ and denote the faces of $\kappa$ $(f_{\kappa}^{j})_{1\leq j\leq
N_{\kappa}}$ where $N_{\kappa}$ is the number of faces of $\kappa$ and the
$f_{\kappa}^{j}\in{\cal F}$ are distinct. Then define the ${\bf
x}_{i}^{\kappa}$ in 12 by ${\bf x}_{(j-1)n_{f}+k}^{\kappa}={\bf
x}_{f^{j}_{\kappa}}({\boldsymbol{\xi}}_{k}^{f^{j}_{\kappa}})$ for $1\leq j\leq
N_{\kappa}$ and $1\leq k\leq n_{f}$ and $N_{f}=N_{\kappa}\times n_{f}$. We
further define Jacobians of the transformations by $J_{\kappa}({\bf x})=|{\bf
x}_{\kappa}^{\prime}({\boldsymbol{\xi}})|$ and $J_{f}({\bf x})=|{\bf
x}_{f}^{\prime}(\xi_{1},\dots,\xi_{d-1})|$.

The DG method consists in defining a discrete weak formulation of 1 by
multiplying it with test functions $\phi_{k}^{\kappa}$ spanning ${\cal
V}_{h}^{p}$ and integrating over $\kappa$, using integration by parts and
approximating ${\bf f}({\bf u}_{h})\cdot{\bf n}$ by two-point numerical fluxes
and the integrals by the quadrature rules. The space discretization in 43
reads

(44) ${\bf R}_{k}^{\kappa}({\bf
u}_{h})=-\sum_{i=1}^{N_{v}}\omega_{i}^{V}J_{\kappa}({\bf y}_{i}^{\kappa}){\bf
f}\big{(}{\bf u}_{h}({\bf
y}_{i}^{\kappa},t^{(n)})\big{)}\cdot\nabla\phi_{k}^{\kappa}({\bf
y}_{i}^{\kappa})\\\ +\sum_{i=1}^{N_{f}}\omega_{i}^{f}J_{f}({\bf
x}_{i}^{\kappa}){\bf h}\big{(}{\bf u}_{h}^{-}({\bf
x}_{i}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf x}_{i}^{\kappa},t^{(n)}),{\bf
n}_{i}^{\kappa}\big{)}\phi_{k}^{\kappa}({\bf x}_{i}^{\kappa}).$

We use an orthonormal basis such that
$M_{\kappa}^{k}=|\kappa|:=\sum_{i=1}^{N_{v}}\omega_{i}^{V}J_{\kappa}({\bf
y}_{i}^{\kappa})$ and further seting $\phi_{k}^{\kappa}=1_{\kappa}$ the
indicator function of $\kappa$, the first sum vanishes and we obtain 12 with
$s_{k}^{\kappa}=\frac{\omega_{k}^{f}J_{f}({\bf x}_{k}^{\kappa})}{|\kappa|}$
and $\beta_{k}^{\kappa}=\varepsilon_{\kappa}s_{k}^{\kappa}$ in 14 where
$\varepsilon_{\kappa}$ is evaluated from 33. Then by 34, the scheme 43 and 44
is IDP under the condition

(45) $\Delta
t^{(n)}\max_{\kappa\in\mathcal{T}_{h}}\frac{1}{\varepsilon_{\kappa}}\max_{1\leq
k\leq N_{f}}\max\Big{(}\lambda^{\star}_{\kappa},|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)})\big{)}\Big{)}\leq\frac{1}{2}.$

###  4.2 Discontinuous Galerkin Spectral Element Method

In the DGSEM, the reference element is an hypercube :
$\hat{K}=I^{d}:=\\{{\boldsymbol{\xi}}=(\xi_{1},\dots,\xi_{d}):\;-1\leq\xi_{j}\leq
1\\}$ and the polynomial space $\mathcal{P}^{p}(I^{d})$ is formed by tensor
products of polynomials of degree at most $p$ in each direction. The
approximate solution is sought under the form 42 where $({\bf
U}^{\kappa}_{k})_{1\leq k\leq N_{p}}$ are the $N_{p}=(p+1)^{d}$ DOFs in the
element $\kappa$ with indexing

$k=\overline{k}(i_{1},\dots,i_{d}):=1+\sum_{j=1}^{d}i_{j}(p+1)^{j-1}\quad
0\leq i_{1},\dots,i_{d}\leq p.$

We define a basis $(\phi^{\kappa}_{k})_{1\leq k\leq N_{p}}$ of ${\cal
V}_{h}^{p}(\kappa)$ by using tensor products: $\phi^{\kappa}_{k}({\bf
x})=\phi^{\kappa}_{k}({\bf
x}_{\kappa}({\boldsymbol{\xi}}))=\Pi_{j=1}^{d}\ell_{i_{j}}(\xi_{j})$, where
$\ell_{0\leq i\leq p}$ denote the $i$th Lagrange interpolation polynomial
associated to $\zeta_{i}$ the $i$th Gauss-Lobatto quadrature node with
$\zeta_{0}=-1<\zeta_{1}<\dots<\zeta_{p}=1$ (see Fig. 2) and by $\omega_{i}$ we
denote the associated weight. We therefore have the following cardinality
relation at quadrature points
${\boldsymbol{\xi}}_{k^{\prime}}=(\xi_{i^{\prime}_{1}},\dots,\xi_{i^{\prime}_{d}})$
in $\hat{K}$: $\phi^{\kappa}_{k}({\bf
x}_{k^{\prime}}^{\kappa})=\phi^{\kappa}_{k}({\bf
x}_{\kappa}({\boldsymbol{\xi}}_{k^{\prime}}))=\delta_{i_{1},i^{\prime}_{1}}\dots\delta_{i_{d},i^{\prime}_{d}}$
with $\delta_{i,i^{\prime}}$ the Kronecker symbol, so the DOFs correspond to
the point values of the solution: ${\bf U}_{k}^{\kappa}(t)={\bf u}_{h}({\bf
y}_{k}^{\kappa},t)$ and interpolation and quadrature points are collocated,
hence $N_{v}=N_{p}$.

$\kappa=\kappa^{-}$$\kappa^{+}$$f$${\bf u}_{h}^{-}$${\bf u}_{h}^{+}$${\bf
n}_{k}^{\kappa}$${\bf x}_{k}^{\kappa}$${\bf
y}_{i}^{\kappa}$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$
Figure 2: Notations for the DGSEM and $d=2$: inner and outer elements,
$\kappa^{-}$ and $\kappa^{+}$, for $d=2$; definitions of traces ${\bf
u}_{h}^{\pm}$ on the interface $f$ and of the unit outward normal vector ${\bf
n}_{k}^{\kappa}$. Element quadrature nodes ${\bf y}_{i}^{\kappa}$ and surface
quadrature node ${\bf x}_{k}^{\kappa}$ that are also included in the ${\bf
y}_{i}^{\kappa}$.

Let introduce the discrete derivative matrix with entries
$D_{ij}=\ell_{j}^{\prime}(\zeta_{i})$ with $0\leq i,j\leq p$. The DGSEM
discretization takes the form 43 with
$M_{k}^{\kappa}=\omega_{k}^{V}J_{\kappa}({\bf k}_{k}^{\kappa})$,
$\omega_{k}^{V}=\Pi_{j=1}^{d}\omega_{i_{j}}$ for
$k=\overline{k}(i_{1},\dots,i_{d})$ and

$\displaystyle{\bf R}_{k}^{\kappa}({\bf u}_{h})$
$\displaystyle=2\omega_{k}^{V}\sum_{j=1}^{d}\sum_{l=0}^{p}D_{i_{j}l}{\bf
h}_{sym}\big{(}{\bf U}_{k}^{\kappa},{\bf
U}_{k^{\prime}_{j}}^{\kappa},\\{J_{\kappa}\nabla\xi_{j}\\}_{(k,k^{\prime}_{j})}\big{)}$
(46) $\displaystyle+\sum_{i=1}^{N_{f}}\phi_{k}^{\kappa}({\bf
x}_{i}^{\kappa})\omega_{i}^{f}J_{f}({\bf x}_{i}^{\kappa})\Big{(}{\bf
h}\big{(}{\bf U}_{k,n}^{\kappa},{\bf u}_{h}^{+}({\bf x}_{i}^{\kappa},t),{\bf
n}_{i}^{\kappa}\big{)}-{\bf f}\big{(}{\bf U}_{k,n}^{\kappa}\big{)}\cdot{\bf
n}_{i}^{\kappa}\Big{)}$

where for $d=3$ $k^{\prime}_{1}=\overline{k}(l,i_{2},i_{3})$,
$k^{\prime}_{2}=\overline{k}(i_{1},l,i_{3})$ and
$k^{\prime}_{3}=\overline{k}(i_{1},i_{2},l)$, hence
$k^{\prime}_{j}=k+(i^{\prime}_{j}-i_{j})(p+1)^{j-1}$, $N_{f}=2d(p+1)^{d-1}$,
$\omega_{i}^{f}=\Pi_{j=1}^{d-1}\omega_{i_{j}}$, and

$\begin{array}[]{c}\\{J_{\kappa}\nabla{\bf\xi}\\}_{(k,k^{\prime}_{j})}=\frac{1}{2}\big{(}J_{\kappa}({\bf
y}_{k}^{\kappa})\nabla\xi_{j}({\boldsymbol{\xi}}_{k})+J_{\kappa}({\bf
y}_{k^{\prime}_{j}}^{\kappa})\nabla\xi_{j}({\boldsymbol{\xi}}_{k^{\prime}_{j}})\big{)},\end{array}$

have been introduced to keep conservation of the scheme [46]. By ${\bf
h}_{sym}$ we denote a two-point flux supposed to be symmetric in the sense
that ${\bf h}_{sym}({\bf u},{\bf v},{\bf n})={\bf h}_{sym}({\bf v},{\bf
u},{\bf n})$. Note that $\phi_{k}^{\kappa}({\bf x}_{i}^{\kappa})=1$ if ${\bf
x}_{i}^{\kappa}={\bf y}_{k}^{\kappa}$ and $\phi_{k}^{\kappa}({\bf
x}_{i}^{\kappa})=0$ else.

The choice of the ${\bf y}_{k}^{\kappa}$ in the quadrature 14 is not unique
and we here use the following decomposition

$\displaystyle\langle{\bf u}_{h}^{(n)}\rangle_{\kappa}$
$\displaystyle=\frac{1}{|\kappa|}\sum_{{\bf
y}_{i}^{\kappa}\in\kappa}\omega_{i}^{V}J_{\kappa}({\bf y}_{i}^{\kappa}){\bf
U}_{i,n}^{\kappa}$ $\displaystyle=\frac{1}{|\kappa|}\sum_{{\bf
y}_{i}^{\kappa}\in\text{int}(\kappa)}\omega_{i}^{V}J_{\kappa}({\bf
y}_{i}^{\kappa}){\bf
U}_{i,n}^{\kappa}+\frac{1}{|\kappa|}\sum_{i=1}^{N_{f}}\tilde{\omega}_{i}^{f}J_{\kappa}({\bf
x}_{i}^{\kappa}){\bf u}_{h}({\bf x}_{i}^{\kappa},t^{(n)}),$

where $\text{int}(\kappa)$ denotes the interior of $\kappa$, while
$\tilde{\omega}_{i}^{f}J_{\kappa}({\bf
x}_{i}^{\kappa})=\frac{1}{d}\omega_{j}^{V}J_{\kappa}({\bf y}_{j}^{\kappa})$ if
${\bf y}_{j}^{\kappa}={\bf x}_{i}^{\kappa}$ is a vertex of the $d$-dimensional
hexahedron, $\tilde{\omega}_{i}^{f}J_{\kappa}({\bf
x}_{i}^{\kappa})=\frac{1}{d-1}\omega_{j}^{V}J_{\kappa}({\bf y}_{j}^{\kappa})$
if ${\bf y}_{j}^{\kappa}={\bf x}_{i}^{\kappa}$ is on some edge, and
$\tilde{\omega}_{i}^{f}J_{\kappa}({\bf
x}_{i}^{\kappa})=\omega_{j}^{V}J_{\kappa}({\bf y}_{j}^{\kappa})$ else.

Summing 43 over $1\leq k\leq N_{v}$ gives for the cell-averaged solution

$\displaystyle\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa}$
$\displaystyle=\langle{\bf u}_{h}^{(n)}\rangle_{\kappa}-\frac{\Delta
t^{(n)}}{|\kappa|}\sum_{k=1}^{N_{v}}{\bf R}^{\kappa}_{k}({\bf u}_{h}^{(n)})$
$\displaystyle=\langle{\bf u}_{h}^{(n)}\rangle_{\kappa}-\frac{\Delta
t^{(n)}}{|\kappa|}\sum_{i=1}^{N_{f}}\omega_{i}^{f}J_{f}({\bf
x}_{i}^{\kappa}){\bf h}\big{(}{\bf u}_{h}^{-}({\bf
x}_{i}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf x}_{i}^{\kappa},t^{(n)}),{\bf
n}_{i}^{\kappa}\big{)},$

by conservation of the DGSEM [15, 46] and providing that the so-called metric
identities are satisfied at the discrete level [29]. This relation can be
identified with 12 with $s_{k}^{\kappa}=\frac{\omega_{k}^{f}J_{f}({\bf
x}_{k}^{\kappa})}{|\kappa|}$. Then we can apply Theorem 3.8 with
$\beta^{\kappa}_{k}=\frac{\tilde{\omega}_{k}^{f}J_{\kappa}({\bf
x}_{k}^{\kappa})}{|\kappa|}$ and the DGSEM scheme 43 is IDP under the
condition

(47) $\Delta t^{(n)}\max\limits_{\kappa\in\Omega_{h}}\max\limits_{1\leq k\leq
N_{f}}\frac{\omega_{k}^{f}J_{f}({\bf
x}_{k}^{\kappa})}{\tilde{\omega}_{k}^{f}J_{\kappa}({\bf
x}_{k}^{\kappa})}\max\Big{(}\lambda^{\star}_{\kappa},|\lambda|\big{(}{\bf
u}_{h}^{-}({\bf x}_{k}^{\kappa},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{k}^{\kappa},t^{(n)})\big{)}\Big{)}\leq\frac{1}{2}.$

###  4.3 Other methods

Properties 12, 14, and 22 also hold for other discretely conservative spectral
difference methods on general curved elements provided the discretization
operators satisfy the metric identities at the discrete level, which imposes
some limits on the order of approximation of the mesh elements compared to the
approximation order of the solution [29]. The limiter 36 may hence be applied
to make these methods invariant domain preserving. We list some examples
below.

The skew-symmetric entropy stable modal discontinuous Galerkin methods [6]
uses skew-hybrized summation-by-parts (SBP) operators allowing conservation
and free-stream preservation under the standard accuracy requirements of
volume and surface quadratures. In [8] multidimensional discretization schemes
based on SBP operators on general curved elements generalize the staggered
finite differences from [11]. The discretizations with curved elements remain
accurate, conservative, and entropy stable. Spectral differences [33] and
staggered Chebyshev [30] methods on curved elements belong to a same family of
conservative approximations satisfying the discrete metric identities. The
methods use two sets of interpolation points for the solution and fluxes and
impose the discrete residuals to be satisfied at solution points. Flux points
contain points at faces of the elements where two-point numerical fluxes are
used. Then, the space derivatives at solution points are evaluated by
differencing the polynomials interpolating the fluxes.

##  5 Numerical experiments

Let consider the compressible Euler equations of gas dynamics. The
conservative variables and fluxes in 1 are

(48) ${\bf u}=\begin{pmatrix}\rho\\\ \rho{\bf v}\\\ \rho E\\\
\end{pmatrix},\quad{\bf f}=\begin{pmatrix}\rho{\bf v}^{\top}\\\ \rho{\bf
v}{\bf v}^{\top}+\mathrm{p}{\bf I}_{d}\\\ (\rho E+\mathrm{p}){\bf v}^{\top}\\\
\end{pmatrix},$

where $\rho$, ${\bf v}$, and $E$ denote the density, velocity vector, and
specific total energy, respectively. The system is closed by defining the
equation of state $\mathrm{p}=\mathrm{p}(\frac{1}{\rho},e)$ with
$e=E-\frac{1}{2}{\bf v}\cdot{\bf v}$ the specific internal energy and the
system is hyperbolic over the set of states $\Omega^{a}=\\{{\bf
u}\in\mathbb{R}^{d+2}:\;\rho>0,{\bf v}\in\mathbb{R}^{d},e>0\\}$. We focus here
on the polytropic ideal gas law $\mathrm{p}=(\gamma-1)\rho e$ where
$\gamma=\frac{C_{p}}{C_{v}}=\frac{7}{5}$ is the ratio of specific heats. The
compressible Euler equations 1 and 48 possess the natural entropy – entropy
flux pair

$\eta=-\rho\mathrm{s},\quad{\bf q}=-\rho\mathrm{s}{\bf
v},\quad\mathrm{s}=C_{v}\ln\Big{(}\frac{p}{\rho^{\gamma}}\Big{)},$

and $\mathcal{B}=\\{{\bf u}\in\Omega^{a}\;\mathrm{s}({\bf
u})\geq\mathrm{s}_{0}\\}$, with $\mathrm{s}_{0}$ in $\mathbb{R}$, is an
invariant domain for 1 and 48 [13]. We use our convex limiting strategy with
the quasiconcave functions $\psi_{1}\equiv\rho$ and $\psi_{2}\equiv\rho e$.

We now test the robustness and efficiency of the CFL condition on simulations
of 1 and 48 with discontinuous solutions on one-dimensional and unstructured
two-dimensional grids. We use the modal DG method in Section 4.1 and the DGSEM
in Section 4.2. The CFL conditions 45 and 47 guaranty the cell-averaged
solution to be IDP and we then apply the limiter 36 to further impose the
high-order solution to be IDP. We will compare the following limiting
strategies: a positivity limiter (POS) which imposes the solution to remain in
$\Omega^{a}$ thus extending [50] to unstructured grids; an IDP limiter which
imposes the solution to remain in the convex hull of the states in 31 and
computing the time step in 45 and 47 with either the global wave estimate
$\lambda_{\kappa}^{\star}$ (IDP) from algorithm 37, or the local wave estimate
$\lambda_{\kappa}^{\star}$ from 40 (IDPloc). Imposing the IDP property may
result in over-limiting of the solution and some strategies are usually
applied such as bound relaxation [21], or subcell smoothness indicator [34].
We here follow the second strategy which relies on the smoothness indicator
from [35] (see [34, Sec. 4.4] for details). Finally, we use the Suliciu
pressure relaxation based numerical flux from [5, Sec. 2.4.6] at interfaces,
while for ${\bf h}_{sym}$ in the DGSEM scheme 46 we use the Kennedy and
Grubber splitting from [15].

###  5.1 Riemann problems

We here consider Riemann problems 4 with initial data given in Tab. 1. We
first consider computations with the DGSEM (see Section 4) and the three
limiting strategies. Results are shown in Figs. 3, 4, and 5. The POS limiter
only ensures that the solution remains in the set of states $\Omega^{a}$ and
does not modify non-physical oscillations, with the IDP and IDPloc strategies
succeed in damping spurious oscillations and result in very close
oscillations. All numerical experiments always show that there is no sensible
improvement to evaluate the pseudo-equilibrium state ${\bf
u}_{\kappa}^{\star}$ with local wave estimates $\lambda_{\kappa}^{\star k}$ in
37 instead of a global estimate $\lambda_{\kappa}^{\star}$ in 40, while it
leads to a more expensive algorithm. Besides, our observation show that using
a local estimates require more iterations for algorithm 40 to converge with a
global average between $2.8$ and $3.2$ iterations evazluated over the whole
computations compared to between $1.11$ and $1.17$ when using 37. In the
latter case, for most of the computations the initial guess ${\bf
u}^{\star}_{0}$ given in 37 satisfies the requirement 38 and 39 so no more
step is needed.

Table 1: Initial conditions of Riemann problems 4 where $x_{0}$ indicates the abscissa separating the states. problem | left state $(\rho_{L},u_{L},p_{L})^{\top}$ | right state $(\rho_{R},u_{R},p_{R})^{\top}$ | $x_{0}$
---|---|---|---
Sod [42] | $(1,0,1)^{\top}$ | $(0.125,0,0.1)^{\top}$ | $0$
Lax | $(0.445,0.698,3.528)^{\top}$ | $(0.5,0,0.571)^{\top}$ | $0$
Toro $4$ [44] | $(5.99924,19.5975,460.894)^{\top}$ | $(5.99242,-6.19633,46.0950)^{\top}$ | $-0.1$

(a) POS

(b) IDP

(c) IDPloc

Figure 3: DGSEM computations of the Sod problem at $t=0.2$ with $p=3$ and
$N=100$ elements for the density (top), velocity (middle), and pressure
(bottom).

(a) POS

(b) IDP

(c) IDPloc

Figure 4: DGSEM computations of the Lax problem at $t=0.2$ with $p=3$ and
$N=100$ elements for the density (top), velocity (middle), and pressure
(bottom).

(a) POS

(b) IDP

(c) IDPloc

Figure 5: DGSEM computations of the Toro 4 problem at $t=0.2$ with $p=3$ and
$N=100$ elements for the density (top), velocity (middle), and pressure
(bottom).

We now reproduce the same numerical experiments but using the classical modal
DG (see Section 4.1) with the IDP strategy. Fig. 6 presents the results for
the threes Riemann problems and we observe a good resolution of the waves with
only lower amplitude oscillations compared to the DGSEM results. Imposing the
IDP property at all quadrature points may result in stronger limiting of the
solution with the modal DG scheme as these are more quadrature points (compare
Figs. 1 and 2). Again, algorithm 37 converges very fast with a global average
between of $1.11$ and $1.18$ iterations evaluated over the whole computations.
Finally note that all the computations (with either DGSEM, or modal DG scheme)
require to apply one of the limiting strategy to avoid non-physical solution.

(a) Sod

(b) Lax

(c) Toro 4

Figure 6: Modal DG (IDP) computations of Riemann problems with $p=3$ and
$N=100$ elements for the density (top), velocity (middle), and pressure
(bottom).

###  5.2 Double Mach Reflection problem

We now consider the two-dimensional problem of a Mach $10$ shock reflection
over a $30^{\circ}$ wedge [47]. Ahead of the shock, the gaz is at rest and has
a density of $1.4$ and pressure of $1$. Inflow and outflow conditions are
applied at the left and bottom boundaries, while a symmetry condition is
applied at the top boundary. Initially, the shock is located at $x=0$
corresponding to the beginning of the wedge. We use an unstructured mesh with
$132800$ quadrangles to solve the horizontally moving shock interacting with
the inclined wall where slip conditions are applied. In Fig. 7 we can observe
qualitatively similar results for the three tests. The limiting strategy
induce some spurious oscillations, but the computations proved to be robust.
Here again, algorithm 37 proves to be cheap in term of iterations to converge
with global averages of $1.18$ for both the DGSEM (IDP) and modal DG (IDP)
computations.

(a) POS (DGSEM)

(b) IDP (DGSEM)

(c) IDPloc (DGSEM)

(d) IDP (modal DG)

Figure 7: Comparison of the results for the double Mach reflection problem: 39
equispaced density contours between 2.1 and 22.1.

##  6 Conclusions

We here investigate robustness and stability properties of discretely
conservative high-order spectral discontinuous methods with explicit time
stepping for the approximation of hyperbolic systems of conservation laws. We
derive a condition on the time step to guaranty that the cell-averaged
approximate solution is a convex combination of DOFs at preceding time step
and updates of invariant domain preserving and entropy stable three-point
schemes. As a consequence, the cell-averaged solution lies in some convex
invariant domain and we apply a posteriori scaling limiting techniques [50]
that impose to all the DOFs to satisfy the same invariant domain properties.

The condition on the time step is evaluated from the traces of the solution at
faces of the mesh and can be easily evaluated on the fly. Provided the scheme
satisfies the discrete metric identities, the condition is fairly general and
holds for general unstructured grid with possibly curved elements. It relies
on the existence of a so-called pseudo-equilibrium state which satisfies a
flux balance over each mesh element, and the existence of a quadrature rule
including the traces to evaluate the cell-averaged solution. We here prove
their existence in the general case and provide an iterative algorithm to
evaluate the pseudo-equilibrium state. We illustrate these results with the
classical modal discontinuous Galerkin and DGSEM schemes. Numerical
experiments in one and two space dimensions are provided to illustrate the
robustness and stability of the present approach. The extension of this
framework to parabolic systems of conservation laws, e.g., the compressible
Navier-Stokes, is a possible direction of future research.

## References

  * [1] M. Avriel, r-convex functions, Math. Program., 2 (1972), pp. 309–323.
  * [2] G. E. Barter and D. L. Darmofal, Shock capturing with pde-based artificial viscosity for dgfem: Part i. formulation, J. Comput. Phys., 229 (2010), pp. 1810–1827, https://doi.org/https://doi.org/10.1016/j.jcp.2009.11.010.
  * [3] F. Bassi and S. Rebay, High-order accurate discontinuous finite element solution of the 2d Euler equations, J. Comput. Phys., 138 (1997), pp. 251–285.
  * [4] J. P. Boris and D. L. Book, Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works, J. Comput. Phys., 11 (1973), pp. 38–69, https://doi.org/https://doi.org/10.1016/0021-9991(73)90147-2.
  * [5] F. Bouchut, Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws and Well-Balanced Schemes for Sources, Frontiers in Mathematics, Birkhäuser Basel, 2004.
  * [6] J. Chan, Skew-symmetric entropy stable modal discontinuous Galerkin formulations, J. Sci. Comput., 81 (2019), pp. 459–485, https://doi.org/https://doi.org/10.1007/s10915-019-01026-w.
  * [7] B. Cockburn and C. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws ii: general framework, Math. Comput., 52 (1989), pp. 411–435, https://doi.org/https://doi.org/10.1090/S0025-5718-1989-0983311-4.
  * [8] J. Crean, J. E. Hicken, D. C. Del Rey Fernández, D. W. Zingg, and M. H. Carpenter, Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements, J. Comput. Phys., 356 (2018), pp. 410–438, https://doi.org/https://doi.org/10.1016/j.jcp.2017.12.015.
  * [9] R. J. Diperna, Convergence of approximate solutions to conservation laws, Arch. Rat. Mech. Anal., 82 (1983), pp. 27–70.
  * [10] W. S. Don, Numerical study of pseudospectral methods in shock wave applications, J. Comput. Phys., 110 (1994), pp. 103–111.
  * [11] T. C. Fisher and M. H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comput. Phys., 252 (2013), pp. 518–557.
  * [12] M. Franco, J.-S. Camier, J. Andrej, and W. Pazner, High-order matrix-free incompressible flow solvers with GPU acceleration and low-order refined preconditioners, Comput. Fluids, 203 (2020), p. 104541.
  * [13] H. Frid, Maps of convex sets and invariant regions for finite-difference systems of conservation laws, Arch. Rational Mech. Anal., 160 (2001), pp. 245––269, https://doi.org/https://doi.org/10.1007/s002050100166.
  * [14] G. J. Gassner, A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods, SIAM J. Sci. Comput., 35 (2013), pp. A1233–A1253, https://doi.org/10.1137/120890144.
  * [15] G. J. Gassner, A. R. Winters, and D. A. Kopriva, Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations, J. Comput. Phys., 327 (2016), pp. 39–66.
  * [16] S. Godunov, A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics, Math. USSR Sbornik, 47 (1959), pp. 271–306.
  * [17] D. Gottlieb and C.-W. Shu, On the Gibbs phenomenon and its resolution, SIAM Rev., 39 (1997), pp. 644–668.
  * [18] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [19] J.-L. Guermond, M. Nazarov, B. Popov, and I. Tomas, Second-order invariant domain preserving approximation of the Euler equations using convex limiting, SIAM J. Sci. Comput., 40 (2018), pp. A3211–A3239, https://doi.org/10.1137/17M1149961.
  * [20] J.-L. Guermond, R. Pasquetti, and B. Popov, Entropy viscosity method for nonlinear conservation laws, J. Comput. Phys., 230 (2011), pp. 4248–4267, https://doi.org/https://doi.org/10.1016/j.jcp.2010.11.043.  Special issue High Order Methods for CFD Problems.
  * [21] J.-L. Guermond, B. Popov, and I. Tomas, Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems, Comput. Methods Appl. Mech. Engrg., 347 (2019), pp. 143–175, https://doi.org/https://doi.org/10.1016/j.cma.2018.11.036.
  * [22] A. Harten, P. D. Lax, and B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev., 25 (1983), pp. 35–61.
  * [23] A. Hiltebrand and S. Mishra, Entropy stable shock capturing space–time discontinuous Galerkin schemes for systems of conservation laws, Numer. Math., 126 (2014), pp. 103–151.
  * [24] D. Hoff, Invariant regions for systems of conservation laws, Trans. Amer. Math. Soc., 289 (1985), pp. 591–610, https://doi.org/https://doi.org/10.2307/2000254.
  * [25] M. Hutchinson, A. Heinecke, H. Pabst, G. Henry, M. Parsani, and D. Keyes, Efficiency of high order spectral element methods on petascale architectures, in International Conference on High Performance Computing, Springer, 2016, pp. 449–466.
  * [26] J. JAFFRE, C. JOHNSON, and A. SZEPESSY, Convergence of the discontinuous Galerkin finite element method for hyperbolic conservation laws, Math. Models Methods Appl. Sci., 05 (1995), pp. 367–386, https://doi.org/10.1142/S021820259500022X.
  * [27] Y. Jiang and H. Liu, Invariant-region-preserving dg methods for multi-dimensional hyperbolic conservation law systems, with an application to compressible Euler equations, J. Comput. Phys., 373 (2018), pp. 385–409, https://doi.org/https://doi.org/10.1016/j.jcp.2018.03.004.
  * [28] J. S. K. Chueh, C. Conley, Positively invariant regions for systems of nonlinear diffusion equations, Indiana Univ. Math. J., 26 (1977), pp. 373–392.
  * [29] D. A. Kopriva, Metric identities and the discontinuous spectral element method on curvilinear meshes., J. Sci. Comput., 26 (2006), pp. 302–327, https://doi.org/https://doi.org/10.1007/s10915-005-9070-8.
  * [30] D. A. Kopriva and J. H. Kolias, A conservative staggered-grid chebyshev multidomain method for compressible flows, J Comput. Phys., 125 (1996), pp. 244–261, https://doi.org/https://doi.org/10.1006/jcph.1996.0091.
  * [31] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, and J. Flaherty, Shock detection and limiting with discontinuous galerkin methods for hyperbolic conservation laws, Applied Numerical Mathematics, 48 (2004), pp. 323–338, https://doi.org/https://doi.org/10.1016/j.apnum.2003.11.002.
  * [32] P. D. Lax, Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, Society for Industrial and Applied Mathematics, 1973, https://doi.org/10.1137/1.9781611970562, https://arxiv.org/abs/https://epubs.siam.org/doi/pdf/10.1137/1.9781611970562.
  * [33] Y. Liu, M. Vinokur, and Z. Wang, Spectral difference method for unstructured grids i: Basic formulation, J Comput. Phys., 216 (2006), pp. 780–801, https://doi.org/https://doi.org/10.1016/j.jcp.2006.01.024.
  * [34] W. Pazner, Sparse invariant domain preserving discontinuous galerkin methods with subcell convex limiting, Comput. Methods Appl. Mech. Engrg., 382 (2021), p. 113876, https://doi.org/https://doi.org/10.1016/j.cma.2021.113876.
  * [35] P.-O. Persson and J. Peraire, Sub-Cell Shock Capturing for Discontinuous Galerkin Methods, 2006, https://doi.org/10.2514/6.2006-112.
  * [36] B. Perthame and C.-W. Shu, On positivity preserving finite volume schemes for Euler equations, Numer. Math., 73 (1996), pp. 119–130, https://doi.org/10.1007/s002110050187.
  * [37] W. H. Reed and T. R. Hill, Triangular mesh methods for the neutron transport equation, Los Almos Report LA-UR-73-479, (1973).
  * [38] F. Renac, Entropy stable, robust and high-order DGSEM for the compressible multicomponent Euler equations, J. Comput. Phys., 445 (2021), p. 110584, https://doi.org/10.1016/j.jcp.2021.110584.
  * [39] V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, J. Comp. Math. Phys. USSR, 1 (1961), pp. 267–279.
  * [40] D. Serre, Domaines invariants pour les systèmes hyperboliques de lois de conservation, J. Differ. Equ., 69 (1987), pp. 46–62, https://doi.org/10.1016/0022-0396(87)90102-1.
  * [41] C.-W. Shu, High orderWENO and DG methods for time-dependent convection-dominated PDEs: A brief survey of several recent developments, J. Comput. Phys., 316 (2016), pp. 598–613.
  * [42] G. A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, Journal of computational physics, 27 (1978), pp. 1–31.
  * [43] P. D. Thomas and C. K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA J., 17 (1979), pp. 1030–1037, https://doi.org/10.2514/3.61273.
  * [44] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Third Edition, Springer-Verlag Berlin Heidelberg, 2009.
  * [45] Z. J. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. T. Huynh, et al., High-order CFD methods: current status and perspective, International Journal for Numerical Methods in Fluids, 72 (2013), pp. 811–845.
  * [46] N. Wintermeyer, A. R. Winters, G. J. Gassner, and D. A. Kopriva, An entropy stable nodal discontinuous Galerkin method for the two dimensional shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry, J. Comput. Phys., 340 (2017), pp. 200–242.
  * [47] P. Woodward and P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys., 54 (1984), pp. 115–173, https://doi.org/https://doi.org/10.1016/0021-9991(84)90142-6.
  * [48] S. T. Zalesak, Fully multidimensional flux-corrected transport algorithms for fluids, J. Comput. Phys., 31 (1979), pp. 335–362.
  * [49] X. Zhang and C.-W. Shu, On maximum-principle-satisfying high order schemes for scalar conservation laws, J. Comput. Phys., 229 (2010), pp. 3091–3120.
  * [50] X. Zhang and C.-W. Shu, On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, J. Comput. Phys., 229 (2010), pp. 8918–8934.
  * [51] X. Zhang, Y. Xia, and C.-W. Shu, Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes, Journal of Scientific Computing, 50 (2012), pp. 29–62.