Consider a hyperbolic conservation law
in space dimensions such as the compressible Euler equations of gas dynamics, where are the conserved variables, the fluxes, and , the time and space coordinates, respectively. Of course, the conservation law has to be equipped with appropriate initial and boundary conditions.
Given a convex entropy and entropy fluxes fulfilling , smooth solutions of (1) satisfy and the entropy inequality
is used as additional admissibility criterion for weak solutions, cf. [dafermos2010hyperbolic]. Because of the convexity of , the entropy variables and the conservative variables can be used interchangeably.
Ever since the seminal work of Tadmor [tadmor1987numerical], there has been some interest in techniques to mimic (2) for semidiscretisations of hyperbolic conservation laws. Some recent contributions are, e.g. [carpenter2016towards, ranocha2017shallow, wintermeyer2017entropy, bohm2018entropy, friedrich2018entropy, sjogreen2018high, chan2018discretely]. Recently, relaxation Runge–Kutta methods have been proposed to transfer such semidiscrete entropy conservation/dissipation results to fully discrete schemes [ketcheson2019relaxation, ranocha2019relaxation].
In this article, the correction terms enforcing entropy conservation of numerical methods in the general class of residual distribution (RD) schemes proposed by Abgrall [abgrall2018general] and modifications suggested in [ranocha2018strong] are studied deeper. They are characterised as solutions of certain optimisation problems and different weightings are introduced. Additionally, new applications and generalisations are developed and compared in numerical experiments. Finally, the basic idea of such an optimisation approach is applied to grid refinement and coarsening, resulting in entropy stable/dissipative transfer operations.
This article is structured as follows. Firstly, the numerical schemes and entropy correction terms are introduced in Section 2, starting with residual distribution schemes in Section 2.1. Thereafter, discontinuous element based schemes such as discontinuous Galerkin (DG) methods are described in Section 2.2 and the characterisations of entropy correction terms as solutions of certain optimisation problems are developed. Generalisations to entropy inequalities, multiple linear constraints, and kinetic energy preservation for the compressible Euler equations are developed in Section 3. Numerical examples using all these schemes are presented in Section 4. Next, the application of such optimisation ideas to grid refinement and coarsening is developed in Section 5, including numerical examples. Finally, the results of this study are summed up in Section 6, conclusions are drawn and some directions of further research are presented.
2 Entropy Corrections for Numerical Schemes
In this section, existing formulations of entropy correction terms for residual distribution and discontinuous element based schemes are described and an interpretation in terms of a quadratic minimisation problem is presented.
2.1 Nodal Formulation: Residual Distribution Schemes
The first introductions of residual distribution (RD) schemes, also denoted by fluctuation splitting schemes, can be found in Roe’s seminal work [roe1981approximate] and in the paper by Ni [ni1981multiple]. Since then, further developments have been done for generalisation and to reach high order in the discretisation, cf. [abgrall2018general, abgrall2011construction]
and references therein. The main advantage of the RD approach is the abstract formulation of the schemes. One works only with the degrees of freedom (DoFs). The selection of approximation/solution space and the definition of the residuals specifies the scheme completely and thus the properties of the considered methods. Today, the RD ansatz provides a unifying framework including some — if not most — of the up-to-date used high order methods like continuous and discontinuous Galerkin methods and flux reconstruction schemes[abgrall2018general, abgrall2018connection]. Also the schemes which will be considered in Section 2.2 can be recast into this framework. However, to follow the approach of Abgrall [abgrall2018general], where the entropy correction term is introduced for the first time, the same notation and the general RD approach is applied which will be shortly repeated in the following subsection.
Residual Distribution Schemes
For simplicity, we will explain the RD approach only for the steady state problem
of a hyperbolic conservation law (1) with suitable boundary conditions. RD is applied to discretise (3) in space. A discretisation of (3) will be considered and the correction term is working on this. A possible temporal discretisation for time-dependent problems will be not considered in the following part. Several time-integration methods such as deferred correction and Runge–Kutta schemes to discretise (1) fully together with an RD approach can be found in [abgrall2017high2, ricchiuto2010explicit, abgrall2017high] and references therein, but this is not topic of this paper.
The domain is split into subdomains (e.g triangles in two dimensions). The term denotes the generic element of the mesh and is used for the characteristic mesh size. For the boundary elements, is applied. Then, the degrees of freedom (DoFs) are defined in respect to the splitting and the weights in each . For each , the set of DoFs of linear forms acting on the set of polynomials of degree such that the linear mapping is one-to-one and denotes the set of degrees of freedom in all elements.
The solution will be approximated by an element from the space defined as
A linear combination of basis functions will be used to describe the numerical solution
where the coefficients must be found by a numerical method. Therefore, the residuals come finally into play. Now, the RD scheme can be formulated by the following three steps to calculate the coefficients , which are also illustrated in Figure 1:
Define for any the total residual111This term is also referred as fluctuation term in the literature [abgrall2003high]. of (3) as an approximation of
In the following, will be used to denote the discrete evaluation of integrals by some quadrature rule. For continuous Galerkin schemes, one can choose
Other examples are given in [abgrall2018general] and in Example 2.3 below.
Split the total residual into sub-residual for each degree of freedom , so that the sum of all the contributions over an element is the fluctuation term itself, i.e.
where denotes the boundary residual. With (10) the RD scheme is described. To specify completely the method (FV, DG, etc.), the solution space (4) (and therefore also the basis) has to be chosen and the exact definition of the residuals has to be given. Varying selections yield different methods. However, the conservation relation has to be always fulfilled. For any element and any ,
where is the restriction of in the element , is the restriction of on the other side of the local edge/face of , and is the
-th component of the outer unit normal vectorat . In addition, is a consistent numerical flux, i.e. , and summation over repeated indices is implied. Similar equations (11) hold for the boundary residuals, see [abgrall2018general] for details.
Further properties of the scheme can be written in terms of the residuals. Here, the focus lies on the reinterpretation and extension of the entropy correction term introduced in [abgrall2018general].
Entropy Correction Term
In [abgrall2018general], the author presented an approach to construct entropy conservative/stable schemes in a general framework. Therefore, a correction term is added to the scheme at every degree of freedom to ensure that the scheme fulfils discretely the entropy condition (2). In terms of RD, an entropy conservative scheme222An entropy stable semidiscretisation has an inequality in (2). Here, the steady state case (3) is considered. fulfils
where is a numerical entropy flux and is the entropy variable at the Dof .
In addition, these correction terms have to be chosen such that they do not violate the conservation relation. The correction term is added to the residual at every degree of freedom, such that the corrected residual
fulfils the discrete entropy condition (12). In [abgrall2018general], the following correction term is introduced
The relation defines a linear system of equation with always at least two unknowns. It is enough to show that (14) with (15) is a valid solution. Results concerning possible contradictions of the constraints (16) are given below in Remark 2.6.
The conservation relation for the new scheme is guaranteed because of
The entropy condition is satisfied, since
It is obvious that fulfils the entropy condition (12). ∎
The error behaviour of can be controlled through the used numerical quadrature, see [abgrall2018general] for details.
To finish this section, an example describing a specific scheme will be presented.
To express a discontinuous Galerkin scheme in the RD framework, one has to specify the solution space and the residual in detail. The space is given by (4) and the residual is given by
The boundary residuals are
2.2 Operator Formulation: Discontinuous Element Based Schemes
Here, the focus will lie on element based discretisations using an operator formulation as in [ranocha2018strong]. Therefore, the domain is partitioned into non-overlapping elements and the following discrete operators are used on each element (dropping the elemental index for convenience).
A symmetric and positive definite mass matrix , approximating the scalar product via .
Derivative matrices , approximating the partial derivative .
A restriction/interpolation operator, performing interpolation to the boundary nodes at via .
A symmetric and positive definite boundary mass matrix , approximating the scalar product on .
Multiplication operators , , representing the multiplication of functions on the boundary by the -th component of the outer unit normal at .
Together, the restriction and boundary operators approximate the boundary integral with respect to the outer unit normal as in the divergence theorem, i.e.
If the SBP property
is fulfilled, integration by parts (the divergence theorem) is mimicked on a discrete level, cf. [fernandez2014generalized]. For example, some discontinuous Galerkin schemes can be formulated in this way [gassner2013skew], allowing the transfer of stability results established at the continuous level to the discretisation, cf. [svard2014review, fernandez2014review] and references cited therein.
The general semidiscretisations considered here can be written as
where are volume terms discretising and are surface terms implementing interface/boundary conditions weakly. For the -th conserved variable , , the corresponding semidiscretisation is
A central nodal DG scheme using the numerical (surface) fluxes can be obtained by choosing
A flux differencing or split form discretisation using symmetric numerical volume fluxes and numerical surface fluces can be obtained using [fisher2013high, gassner2016split]
where the upper indices indicate the grid node. Examples using such a notation can be found in [ranocha2018thesis, ranocha2017shallow, ranocha2018numerical] and references cited therein.
The discretisation should be (locally) conservative, i.e. it should satisfy
where is the numerical flux for the -th conserved variable in coordinate direction . An entropy conservative semidiscretisation results in
where a summation over the repeated index is implied and are numerical entropy fluxes corresponding to , cf. [tadmor1987numerical, tadmor2003entropy].
The basic idea of Abgrall [abgrall2018general] is to enforce (28) for any semidiscretisation via the addition of a correction term that is consistent with zero and does not violate the conservation relation (27), resulting in
Using the mass matrix for discrete integration as proposed in [ranocha2018strong], the correction term is
If the denominator of in (30) is zero, the numerical solution is constant in the element because of the Cauchy Schwarz inequality (since and are linearly dependent for each in that case). Then, the discontinuous element scheme reduces to a finite volume scheme using the numerical fluxes and (28) has to hold for entropy conservative numerical fluxes , as described in the following Remark 2.6. Additionally, the term multiplied by vanishes and becomes zero.
Such a correction is only possible, if the constraints (27) and (28) do not contradict each other. In particular, (28) has to hold whenever all are constant (proportional to ). This is satisfied for the schemes investigated in this section, if the numerical surface fluxes are entropy conservative and the corresponding numerical entropy fluxes in the sense of Tadmor [tadmor1987numerical, tadmor2003entropy], i.e. if
where is the flux potential . Indeed, if the numerical solution is constant in an element, DG type schemes such as in Example 2.6 result in a finite volume scheme
which is conservative because of (since the divergence theorem has to hold for constants because of consistency). Additionally,
Hence, it suffices to consider one boundary node. There,
Equation (36) can be reformulated as
Since is symmetric and positive definite, there is a unique solution , which is given by [nocedal1999numerical, Section 16.1]
for some . Hence, must satisfy the constraints and must be in the span of such that the coefficients of are the same for all . This is obviously true for as defined in (30). ∎
Using the same argument used in the proof of Theorem 2.7, one obtains
There are some differences between the role of the correction terms described in this section and the terms used in the previous section. Firstly, since the correction (30) is added to the other side of the hyperbolic equation, the sign differs from the correction term (14) of [abgrall2018general]. Secondly, (30) is a correction for the pointwise time derivative of while (14) is a correction for an integrated version. Loosely speaking, they are related via
Additionally, the role of the indices differs: is a correction term for the -th variable at all grid nodes while is a correction term at the grid node for all variables (if a DG type setting is used for the RD scheme).
Using the notation of RD schemes, the entropy correction term (30) corresponds to
in accordance with (40). Hence, (30) uses an integral weighting (by the quadrature rule) instead of a summation without weights. While (30) is the optimal correction with respect to the mass matrix , its corresponding integral correction (42) is optimal with respect to the norm induced by .
2.3 Finite Difference and Global Spectral Collocation Schemes
Classical single block finite difference and spectral collocation schemes can be interpreted as RD or DG schemes described in Sections 2.1 and 2.2 with one element. In that case, the entropy corrections described above yield globally conservative and globally entropy conservative/stable schemes.
Sadly, global conservation and a global entropy inequality do not imply any sort of convergence towards an entropy solution of scalar conservation laws, even if the scheme converges.
Consider Burgers’ equation
with periodic boundary conditions and the initial condition
The unique entropy solution contains a stationary shock at and a rarefaction wave starting at .
Central periodic finite difference and Fourier collocation schemes can be represented by a skew-symmetric derivative operatorand mass matrix . Hence, there is no difference between the approaches (14) and (30).
Since is constant, a classical central scheme yields a stationary numerical solution. The entropy correction vanishes, too, since is zero (because of periodic boundary conditions). Hence, the same stationary numerical solution is obtained if the entropy correction is added.
If an element based scheme is used instead, the numerical solution with entropy correction term cannot be stationary if the grid is fine enough (obtained by increasing the number of elements), since the difference of numerical entropy boundary fluxes does not vanish if the element contains exactly one initial discontinuity.
2.4 Related Schemes
By Theorem 2.7, the entropy correction terms can be interpreted as solution of a quadratic minimisation problem with equality constraints. The idea of solving such a problem can also be exploited for many different applications.
A similar approach has been used in [hicken2018family] to construct entropy conservative numerical fluxes. In their work, the authors focus on the quadratic optimization problem
where 1.2 1.2 ⋅ ψ¯ff^num
The approach using a quadratic optimisation problem can be extended to more general constraints that are linear for a given state , e.g. some special form of the kinetic energy for the Euler equations, cf. Section 3.3. A combination of split forms and correction terms similar to the ones described hitherto has been presented in [singh2019kinetic].
In this section, some generalisations of the entropy correction terms based on the interpretations as a quadratic optimisation problem will be developed. Here, the notation of Section 2.2 for discontinuous element based schemes will be used.
3.1 Inequality Constraints
In many applications to hyperbolic conservation laws, the main interest lies in an entropy inequality
instead of entropy conservation, resulting in some kind of stability estimates. For example, even if the baseline scheme is not necessarily entropy dissipative in general, it can be dissipative in some cases. In that case, it could be beneficial to preserve this dissipation introduced by the baseline scheme. Moreover, it could be possible to obtain better approximations with smaller corrections if some entropy dissipation is allowed. Instead of (36) in Theorem 2.7, such an optimisation problem is
holds instead of the entropy equality (28).
Since is symmetric and positive definite, there is a unique solution. Apply the active set method described in [nocedal1999numerical, Section 16.4] with feasible initial value of given in (30) to
If the Lagrange multiplier is , the unique solution has been found.
Otherwise, drop the constraint and compute the solution of
which is obviously .
Finally, note that the Lagrange multiplier associated with is , where is given in (30) and satisfies . ∎
3.2 Multiple Linear Constraints
A generalisation of the approach described in Theorem 2.7 to multiple linear constraints is straightforward. This is demonstrated for two constraints
in addition to (27). Here, is the difference of desired and current values for a linear constraint similar to in (30). For example, could be caused by a correction for the kinetic energy for the Euler equations, cf. in (66) below.
It can be desirable to satisfy entropy (in-) equalities for multiple entropies. Based on Remark 2.6, the numerical surface fluxes should be entropy conservative for both entropies. However, this is in general not possible. Indeed, for a scalar conservation law and a fixed entropy, the entropy conservative numerical flux is uniquely determined as
3.3 Kinetic Energy for the Euler Equations
Consider the compressible Euler equations in two space dimensions (the extension to three space dimensions is straightforward)
where is the density of the gas, its speed, the momentum, the specific total energy, and the pressure. The total energy can be decomposed into the internal energy and the kinetic energy , i.e. . For a perfect gas,
where is the ratio of specific heats. For air, will be used, unless stated otherwise.
Following [ranocha2018thesis, Section 7.4], the kinetic energy satisfies
A numerical flux is kinetic energy preserving, if
The corresponding numerical flux for the kinetic energy (approximating the conservative part of the kinetic energy equation) is [ranocha2018thesis, Eq. (7.66)]
Using , a kinetic energy preserving semidiscretisation mimicking (61) has to satisfy (cf. [ranocha2018thesis, Section 7.4])