## 1 Introduction

The objective of this paper is to continue the work started in Clayton et al. [6] where we introduced an explicit and invariant-domain preserving method to approximate the Euler equations equipped with an equation of state that is either tabulated or is given by an expression that is so involved that elementary Riemann problems cannot be solved exactly or efficiently. The method proposed in [6]

is invariant-domain preserving in the sense that it guarantees that the density and the internal energy of the approximation are positive. A key ingredient of the method is a local approximation of the equation of state using a covolume ansatz (Noble-Abel equation of state) from which upper bounds on the maximum wave speed are derived for every elementary Riemann problem. The downside is that the method is first-order accurate in space. In the present paper we propose a technique that increases the accuracy in space to second-order and preserves the invariant-domain properties. We also introduce a functional that can be used as an entropy surrogate. The functional in question is shown to be increasing across shocks in the Riemann problem involved in the construction of the local maximum wave speed. We show how to use this functional to limit the internal energy from below. This feature is useful for equations of state that are tabulated or interpolated from experimental data since in this case no natural notion of entropy is available.

It is sometimes possible, though possibly expensive, to solve Riemann problems when the equation of state is analytic. For instance this is done in Colella and Glaz [7, §1], Ivings et al. [23], Quartapelle et al. [32]. This cannot be done with tabulated equations of state because the information on the pressure is incomplete. Several attempts to develop methods working with an arbitrary or tabulated equation of state have been reported in the literature. One way to do so consists of using approximate Riemann solvers like those found in Dukowicz [8], [7, §2], Roe and Pike [34], Pike [31], and Lee et al. [25]. One can also simplify the Riemann problem by using flux splitting techniques as those reported in Toro et al. [39]. We also refer to Saurel et al. [35], Banks [1], Dumbser and Casulli [9], Dumbser et al. [10] where approximation techniques are developed using approximate Riemann solvers for various equations of state. Some of these techniques guarantee positivity of the density, but little else is guaranteed in general. The method introduced in Clayton et al. [6] is based on a graph viscosity technique using upper bounds on the maximum wave speed in the Riemann problem. Instead of using the two-shock approximation of the Riemann solution, as done in most methods based on approximate solvers, the method proposed in [6]

approximates the pressure in the Riemann fan by the covolume equation of state and subsequently estimates guaranteed upper bounds on the maximum wave speed. This in turn ensures that the internal energy in the proposed algorithm is positive in addition to the density being positive. If it happens that the equation of state is of covolume type, then the method also preserves the minimum principle on the specific entropy. A technique based on similar principles is reported in

Wang and Li [40] where the authors use a stiffened gas equation of state to approximate the pressure in the Riemann fan. To the best of our knowledge, the method proposed in the present paper is among the very first ones that are provably invariant-domain preserving for complex or tabulated equations of state and high-order accurate in space.The paper is organized as follows. In Section 2, we introduce the mathematical model of interest and the corresponding notation. We also briefly discuss the assumptions we make on the equation of state and recall results from [6] that are needed for this work (we recall the first-order approximation of the Euler equations in §2.3). In Section 3

, we construct a provisional update that is higher-order accurate in space. This update is based on a high-order graph viscosity using an entropy commutator and an activation function. Then, in Section

4, we apply a novel convex limiting technique that corrects the invariant-domain violations of the provisional higher-order method. The final result is an approximation technique that is robust, formally second-order accurate in space, provably invariant-domain preserving, and works for every equation of state (tabulated or analytic) that satisfies the mild assumptions stated in §2.2. Finally in Section 5, the method is verified with analytical solutions and published benchmarks and is validated with experiments. A short conclusion is given in §6.## 2 Preliminaries

We formulate the problem and introduce notation in this section. We also recall essential results from [6] for completeness.

### 2.1 The Euler equations

Let be a bounded polyhedron in . Given some initial data and initial time , we look for solving the compressible Euler equations in the weak sense: equationparentequation

(1a) | ||||

(1b) | ||||

(1c) |

The components of the dependent variable

(considered to be a column vector) are the density,

, the momentum, , and the total mechanical energy, . We also introduce the velocity and the specific internal energy, . To simplify the notation later on we introduce the flux , where is the identity matrix.### 2.2 Equation of state

In (1) the function denotes the pressure which we assume
to be given by an *oracle*. This means we assume to have no a priori
knowledge of the function itself apart from some mild structural
assumptions stated below. We only assume that for a given state we
are able to retrieve a pressure in a suitable way, for example by
evaluating an arbitrary analytic expression, or by deriving a value from
tabulated experimental data. More precisely, in the numerical illustrations
reported in §5, we consider oracles given by
analytic functions for the ideal gas, van der Waals, Jones–Wilkins–Lee
and Mie-Gruneisen equations of state. We also use the SESAME
database developed at Los Alamos National Laboratory [27]
to test the method with experimental tabulated data.

Throughout the paper we assume that the domain of definition of the oracle is the set given by

(2) |

We henceforth refer to as the admissible set. One of the
objectives of the paper
is to guarantee that
the approximation is high-order
accurate and leaves invariant.
The inequality found in the definition of is the so called
*maximum compressibility* condition. The constant can be set to
zero if the user has no a priori knowledge about the maximum compressibility
of the fluid under consideration. We recall, however, that a large class of
analytic equations of state, such as the Noble-Abel, the Mie-Gruneisen
(with the Hugoniot locus as the reference curve), and the
Noble-Abel-Stiffened-Gas equations of state all involve a maximum
compressibility constant. Finally, we assume that the oracle
returns a non-negative pressure,

(3) |

The assumption can be weakened, but for the sake of simplicity, we refrain from doing so in this paper. We leave this extension for future works.

### 2.3 First-order time and space approximation

The time and space discretization proposed in Clayton et al. [6] is based on [16]. This method is in some sense a discretization-agnostic generalization of an algorithm introduced by Lax [24, p. 163]. We denote by the current time, , and we denote by the current time step size; that is . Without going into the details of the space approximation, we assume that the current approximation of is a collection of states , where the index set

is used to enumerate all the degrees of freedom of the approximation, and

is in for all . The update at is obtained as follows:(4) |

The quantity is called the lumped mass and we assume that for all . The vector encodes the space discretization. The index set is called the local stencil. This set collects only the degrees of freedom in that interact with (). We assume that the -valued coefficients are such that approximates at some grid point , and the consistency error in space scales optimally with respect to the mesh size for the considered approximation setting. We require that the method be conservative; more precisely, we assume that

(5) |

Concrete expressions for and are given in [19, §4] for continuous and discontinuous finite elements as well as for finite volumes. The computations reported at the end of this paper are done with piecewise linear continuous finite elements. But to stay general, we continue with the abstract discretization-agnostic notation introduced above.

For completeness we now recall how the graph viscosity is defined in [6]. Given and , we set . For , we set , , , and . The non-negativity assumption on the pressure (3) and the assumptions on the density () implies that

(6) |

Then we consider the following Riemann problem:

(7) |

with left data and right data . This problem is well-posed because . Its complete solution is given in Clayton et al. [6, §4]. We denote by the maximum wave speed in this problem. Let be a nontrivial convex subset of . We say that is an invariant set for (7) if for every pair of Riemann data in the solution of (7) takes values in . We then have:

###### Theorem 1 ([6, Thm. 4.6])

###### Remark 1 (Invariant domains)

Theorem 1 asserts that every convex invariant set in
is invariant by the update procedure (4)
with the artificial viscosity defined in (8).
This means that if for all , then
for all . We say that the method is
invariant-domain preserving *for *. Notice in particular that
the method is invariant-domain preserving for the admissible set
. The admissible set may not be the smallest
invariant domain. For instance, if the oracle admits a mathematical
entropy , then the approximation defined above is also
invariant-domain preserving for the set ; see
[16, Cor. 4.2]. By slight abuse of
terminology and provided the context is unambiguous we will simply call a
method invariant-domain preserving without quantifying the precise convex
set.

A source code to compute is available online [5]. A notable drawback of the graph viscosity (8) is that it reduces the space accuracy of the method to first order. The remainder of the paper is concerned with constructing a higher-order approximation and a respective novel convex limiting technique that is invariant-domain preserving.

###### Remark 2 (Pressure approximation)

Notice that the oracle is only invoked to compute the left and right values of the pressure in (7). The pressure in the Riemann fan is approximated by two covolume equations of state. There is one covolume equation of state for each side of the contact discontinuity since on the left of the contact wave and on the right.

## 3 Provisional higher-order method

In this section, we introduce a provisional higher-order method that extends the update (4) to second order when using linear finite elements for the discretization. A convex limiting technique for this provisional update is introduced in §4. The provisional method is based on the work reported in [18, 28].

### 3.1 The method

Higher-order accuracy in space requires using the consistent mass matrix instead of the lumped mass matrix for the discretization of the time derivative. By reducing dispersive errors, the consistent mass matrix is known to yield superconvergence at the grid points; see Christon et al. [4], [14], [28, Sec. 3.4], Thompson [38]. Let be the mass matrix with entries , then can be approximated by where , is the Kronecker symbol, and . This approximation has been shown in [4, 14] to be superconvergent for piecewise-linear continuous finite elements. Let with , then using this approximation we have . Notice that , because

. The skew-symmetry of the summand in

is used in §4 for limiting purposes.Let denote the high-order update (here the superindex reminds us that the update is higher order accurate). Then, for every , the provisional high-order update is given by: equationparentequation

(9a) | |||

(9b) |

The high-order graph viscosity coefficient , defined in Section 3.2, shares the same properties as its low-order counterpart which are necessary for conservation:

(10) |

### 3.2 Entropy commutator

Using the entropy-viscosity methodology introduced in Guermond et al. [17], we now discuss the construction of the high-order graph viscosity coefficient that is used in the high-order update (9b). A complicating factor in this construction is that the pressure is given by an oracle. We thus do not have a priori knowledge on the equation of state and entropies of the system.

Recall that is an entropy pair for the Euler equations if

(11) |

where is the flux of the system (1). Since one may not have access to entropies for tabulated equations of state, we are going to use at every and every time one entropy pair associated with the following flux

(12) |

Here, with (see (6)). We use the following shifted Harten entropy pair in the numerical tests reported in §5:

(13) |

where . Then, we estimate “entropy production” by inserting the approximate solution into a discrete counterpart of (11) which we write as follows:

(14) |

where are the shape functions of the finite element approximation. Notice that the above identity holds true for every smooth function and for every . Substituting into (14), we estimate the local entropy residual as follows:

(15) |

The residual (15) can be thought of as a measure of how well the discrete solution satisfies (14) in each local stencil . We then define the normalized entropy residual as follows: equationparentequation

(16a) | ||||

(16b) |

By definition, the normalized residual has values in and, for piecewise linear finite elements, behaves like where is the typical meshsize. Finally we set

(17) |

where the activation function is defined as follows:

(18) |

and is the rectified linear activation function. Notice that , and for all for a chosen . As a result, one recovers if the entropy residual is larger than . Note that for (the identity also holds for ). When we have . The numerical tests reported in §5 are done with (). An activation function with the same purpose is used in Persson and Peraire [30, Eq. (8)]. Up to a translation, the activation function therein behaves like in the interval .

###### Remark 3 (Entropy shift)

## 4 Convex Limiting

The high-order update (9) is not guaranteed to be invariant-domain preserving; in particular, it is not guaranteed to stay in the admissible set . In order to correct this defect we now discuss a new convex limiting technique that re-establishes invariant-domain preservation for the final high-order update, .

### 4.1 Key observation

A key observation is that one can rewrite (4) as follows: equationparentequation

(19a) | |||

(19b) |

That is, the low-order update is a convex combination (under the appropriate CFL condition) of the local state and the auxiliary states ,

(20) |

The main result established in [6, Thm. 4.6] (and summarized in Theorem 1) is that under the CFL condition stated in Theorem 1 and the definition (8) for , the states are in provided this is already the case of the states . This is done by proving that the states are space averages of the solution to a Riemann problem. We refer the reader to [16, Sec. 3.3], [18, Sec. 3.2] and [19, Sec. 3.2] where this is discussed in detail. We are going to use the states to define local bounds in space and time to perform the limiting of the high-order states .

### 4.2 Entropy surrogate

Our goal is to use the methodology introduced
in [17] to perform the convex limiting of the update
. In this context the use of an oracle
with little a priori knowledge on the equation of state poses a
significant challenge as it makes it impossible to properly define an
entropy. We resolve the impasse by introducing an artificial
*surrogate* entropy that has the right mathematical properties for the
convex limiting methodology to be applied; see
Theorem 2.

For any admissible state and , we define

(21) |

where and are the density and specific internal energy of the state . Furthermore, for every index , we set

(22) |

where , and . The following result is the key motivation for the definition of an entropy surrogate.

###### Lemma 1

###### Proof

We omit the superscript in the proof to simplify the notation. The solution to the extended Riemann problem (7), , is given in Clayton et al. [6]. In particular, we have that if and if . Here, , and is the speed of the contact wave. Let , with the convention that the index is for the left state and is for the right state. Assume that the -wave is a shock wave. Let be the state before the shock and let be an arbitrary state connected to through a shock curve. With the notation for the specific volume (not to be confused with the time step), and since along the wave curve, the Rankine-Hugoniot condition implies that (see Godlewski and Raviart [12, Eq. (4.8), p. 144]), , where by slight abuse of notation we renamed the pressure in (7) by setting . We then infer that only depends on along the shock curve; more precisely, we have

Notice that this function is well defined only on the interval with

(24) |

and it is nonegative on the interval with . Notice that since we assumed that . We now show that the function is nonnegative and increasing on the shock curve for all and . This will then prove the assertion for the choices and , because and

. Setting , we have ; hence, the pressure, , is a monotone increasing function of along shock curves. By definition of shock curves, starting from the state the pressure increases, we conclude that also increases along shock curves, or . Actually, the pressure is finite only in the range ; hence, showing that is nonnegative increasing on shock curves, is equivalent to showing that

is nonnegative decreasing on shock curves. Recall that the specific entropy for the covolume equation of state is given by . A fundamental property of the specific entropy is that it is an increasing function along shocks. That is, is a decreasing function over the interval . This also means that is a decreasing function. We now have

Let be the defined by . Then we have

Computing the derivative of , we see,

Note that

Then as for , , and , and as is a decreasing function. Thus, by the choice of the constants and , we have that and . Hence, is nonnegative increasing on shock curves. This completes the proof.

We now define the surrogate entropy. For all , we set equationparentequation

(25a) |