Analysis of the SBP-SAT Stabilization for Finite Element Methods Part II: Entropy Stability

by   Rémi Abgrall, et al.
Universität Zürich

In the research community, there exists the strong belief that a continuous Galerkin scheme is notoriously unstable and additional stabilization terms have to be added to guarantee stability. In the first part of the series [6], the application of simultaneous approximation terms for linear problems is investigated where the boundary conditions are imposed weakly. By applying this technique, the authors demonstrate that a pure continuous Galerkin scheme is indeed linear stable if the boundary conditions are done in the correct way. In this work, we extend this investigation to the non-linear case and focusing on entropy conservation. Switching to entropy variables, we will provide an estimation on the boundary operators also for non-linear problems to guarantee conservation. In numerical simulations, we verify our theoretical analysis.



There are no comments yet.


page 1

page 2

page 3

page 4


Analysis of the SBP-SAT Stabilization for Finite Element Methods Part I: Linear problems

A pure Galerkin scheme is notoriously unstable. To remedy this issue, st...

Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions for the compressible Navier-Stokes equations

Entropy stable schemes ensure that physically meaningful numerical solut...

A Dissipation Theory for Three-Dimensional FDTD with Application to Stability Analysis and Subgridding

The finite-difference time-domain (FDTD) algorithm is a popular numerica...

An unfitted finite element method using level set functions for extrapolation into deformable diffuse interfaces

We explore a new way to handle flux boundary conditions imposed on level...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In the first paper [6], the authors demonstrated that a pure Galerkin scheme is indeed linear stable if the boundary procedure is done in the correct way. Here, the usage of simultaneous approximation terms (SATs) and imposing the boundary conditions weakly are essential for the investigation.
In this work, we extend this study to the non-linear case and focusing on entropy stability. In [5]

, the author presented a way to build entropy conservative schemes by adding correction term to it. The term works at every degree of freedoms (dofs) and compensates the entropy production/destruction at the dofs but has no influence on the accuracy and the performance of the scheme.

To build now entropy stable Galerkin schemes, we combine both ideas together. We add the correction term to the scheme in the non-linear case and further develop a new SAT boundary procedure to guarantee, in total, entropy stability.
The proposed boundary operators coincide with the classical ones used the linear case [6, 14, 35, 18, 32]) but further extend the SAT approach to non-linear problems. Our boundary approach is not restricted to the continuous Galerkin method but can also be applied using other schemes (e.g. discontinuous Galerkin[18], flux reconstruction [33] or finite differences [14] ). Therefore, the paper is organized as follows:
We explain the residual distribution scheme and the connection to the continuous Galerkin framework. Then, we repeat the concept about entropy where we follow Harten’s approach [23]. Next, we explain how to construct entropy conservative/stable schemes using entropy correction terms developed in [3] and further extended and applied in [5, 7]. These terms will be essential in our further studies. In the next section 3, we shortly summarize the main results about linear stability of the pure Galerkin scheme from [6] and extend in the following section 4 our investigation to the non-linear case. We give a recipe how to construct the boundary operators to guarantee entropy stability for the Galerkin scheme. In numerical simulations, we validate our analysis. Furthermore, we demonstrate also that the correction term can even rescue schemes which are not linear stable by construction. At the end, we summarize our results and give a short outlook.

2 Residual Distribution Schemes

In this section, we shortly introduce/repeat the residual distribution (RD) schemes as they are also well-known in literature, see [12, 1, 4, 3] and references there in. RD provides a unifying framework including some of the up-to-date high order schemes like continuous and discontinuous Galerkin methods and flux reconstruction schemes [3, 5]. The selection of approximation /solution space and the definition of the residuals specifies the scheme completely and thus the properties of the considered methods. In this paper, our focus lies only on the continuous Galerkin scheme. However, we start by repeating the general approach here also to include/present the entropy correction terms as they are originally developed in [3], extended and applied in [5, 7].

2.1 Residual Distribution - Notation and Basic Formulation

We are interested in the numerical approximation of a hyperbolic problem


with suited initial and boundary conditions. Later, we will focus on the boundary condition more precisely, but for the explanation of RD this is not important.

RD will be used for the discretisation in space. For sake of simplicity, we explain the RD approach for the steady state problem of (1) and in the case of a globally continuous approximation. It is

with a Dirichlet condition on the inflow part of the boundary, to fix ideas:

The domain is split into subdomains (e.g triangles/quads in two dimensions, tetrahedrons/hex in 3D). We denote by the generic element of the mesh and by the characteristic mesh size. For the boundary elements, we denote them . Then, the degrees of freedom (DoFs) are defined in each : we have a set of linear forms acting on the set of polynomials of degree such that the linear mapping is one-to-one. The set denote the set of degrees of freedom in all elements.

The solution will be approximated by some element from the space defined by


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 comes finally into play. Before we focus now on the RD scheme and the definition, we want to make two small remarks.

Remark 2.1.

  1. Two cases for the solution space (3) are normally considered and . For the continuity requirement, we need an additional condition on the splitting of the domain e.g. a conformal triangulation while in the second case the conformity can be dropped, see [3, 5]. In this paper, we will make use of .

  2. Furthermore, as basis functions we are working either with Lagrange interpolation where the degrees of freedom are associated to points in

    or Bézier polynomials. Both are basis of but the the Bézier polynomials are a viable option for the unsteady case, see [2]. We note that for any the following condition will be fulfilled

(a) Step 1: Compute the total residual
(b) Step 2: Split the total residual
(c) Step 3: Combine the residuals
Figure 1: Illustration of the three steps(total residual, nodal residuals, gathering residuals) of the RD approach for linear triangular elements.

We specify now the RD scheme by the following three steps:

  1. Define for any the total residual111This term is also referred as fluctuation term in literature [8]. of (2)


    where denotes a quadrature formula to calculate the integrals, is a consistent numerical flux in the outward direction , and is the generic state on the opposite faces of . In the case of a globally continuous approximation, .

  2. 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.

  3. For a face on , we also define boundary residuals which satisfies a conservation relation similar to the previous one, namely

  4. The resulting scheme if finally obtained by summing all sub-residuals of one degree of freedom from different elements . It is


    (6) allows us to calculate the coefficients in our numerical approximation (4).

    In case of , (4) will be split for any (DoF) and we can finally write the discretization of (2). For any , it is


The RD scheme is described by (7) where the second sum is empty when , in which case we get (6). In the above mentioned literature it is shown that some of the nowadays most favourable schemes either based on finite volume or finite element approaches (e.g. SUPG, DG, FR, FV-WENO, etc.) can be recast with this framework (7).

With (7) we obtain the RD scheme, what is left over is the exact definition of the residuals which will finally specifies the method and the framework we are working in, the continuous Galerkin scheme. Therefore, we chose for the solution space and the residuals are defined by


where we used Gauss theorem in the second line. The whole spatial discretisation is now determined. The extension to the unsteady case can be done in a similar way as it is described inter alia in [34, 4]. As it is well known that any finite element technique used to semidiscretize (1) will yield to a formulation



denotes the vector of degrees of freedom,

is the approximation of (here: in RD formulation) and is a mass matrix222In the finite difference community is called norm matrix and is classically abbreviated with , c.f. [35, 27].. In case of continuous elements, this matrix is sparse but is not block diagonal, contrarily to the discontinuous Galerkin methods. It is well-known that the continuous Galerkin scheme suffers from his stability issues. Therefore, it is common to add stabilization terms to the scheme as the are introduced for example in [9] and which are applied in the RD framework already in [1, 4, 3]. However, we follow a different approach in this paper and will renounce these classical stabilisation techniques. In order to do this, we still need some results known from the literature, which we will briefly repeat here.

2.2 Entropy framework

Since, especially nonlinear conservation laws may have infinity weak solution, one has to find additional criterion to select the physically relevant solution among all the weak solutions. This criterion is based on the concept of entropy, see [21, 25]. In the maravous work [23], Harten described the symmetrizability of systems of conservation laws which possses an entropy function and repeated also some well-known conditions and properties of entropy functions.
A scalar function is an entropy function for (1) if

  1. The function satisfies


    where is some scalar function called entropy flux in the direction.

  2. The function is a convex function of .

is the entropy pair with entropy flux and a weak solution of (1) is called entropy solution if satisfies additionally


for all entropies of (1).
In the following, we assume to know a entropy function of (1) and is called entropy variable. It is given by . Because of the convexity of , the mapping is one-to-one and instead of working with in (1) we can work directly with by setting . Then, by introducing the new variable in place of a symmetrization of (1) will be accomplished, see [23, 26] for details. In abuse of notation we won’t introduce a new notation for this change of the variable, but we can also rewrite (1) and (10) in terms of the entropy variable. Also by following [23, 36] we know the following relations are fulfilled where we write everything in terms of the entropy variable


where is the flux potential. For an better understanding, we shortly repeat the main steps of the relations (10)-(13) from [23]. (10) will be derived by using (12) and (13)

Finally, the concept of entropy has been introduced the consideration of numerical schemes [15] and one tries to find/ construct schemes which fulfill a discrete counterpart of (11) which can be written in term of the RD framework


where (15) represents a consistent numerical entropy flux [3]. One speaks about entropy conservative schemes if the equality in (14) holds. In case of an inequality, the scheme is called entropy stable.

2.3 Entropy Correction Terms

In [3], the author present an ansatz to build entropy conservative / stable schemes in a general framework. A correction term is added to the scheme at every degree of freedom and will ensure that the schemes fulfill the discrete entropy condition, but simultaneously has no effect on the conservation relation. In [5, 7], a re-interpretation and applications of these terms can be found. However, a pure continuous Galerkin scheme is not entropy stable at all. Therefore, the correction term is added to the scheme and we shorty repeat the main idea of his procedure from [3].
Therefore, represents the approximation of the entropy variable in . As we already mentioned the mapping is one-to-one and we work directly with . We are adding the correction term to our residual at every degree of freedom. It is

where the discrete entropy condition


holds. In [3], the following correction term is introduced


with is the number of degrees of freedom in ,

Adding (16) to our residual (here: continuous Galerkin) will provide that our constructed scheme fulfills the entropy condition (15) and, by construction of , the conservation relation is guaranteed, since

is satisfied. The error behavior of can be controlled through the used numerical quadrature, see [3] for details and later.

Remark 2.2.

The correction term can be seen also as the approximation of a second order derivation and a further stabilisation term. However, the term works directly on the degrees of freedom and guarantee that the semidiscrete entropy inequality is fulfilled. Through the conservation relation and the application of high order quadrature formulas the accuracy of the scheme is not affected. Actually, these terms are only compensating the entropy production/destruction of the original scheme in the semidiscrete setting.

3 Linear Stability for a pure Galerkin Scheme

In [6], linear stability was investigated for a pure Galerkin scheme using simultaneous approximation terms (SATs). The origin of the SAT approach lies in the finite difference community, together with summation by parts (SBP) operators they are used to prove linear stability [28, 35]. The main idea is to impose the boundary conditions weakly where also boundary operators (penalty terms) are involved. The operators are determined to guarantee linear stability in the semi-discrete setting.
The SBP-SAT technique has already been applied in the discontinuous Galerkin (DG) and Flux Reconstruction (FR) framework to prove stability results, see inter alia [18, 19, 11, 10, 33] and references therein. However, to the best of our knowledge, we were the first ones who applied this technique together with a pure Galerkin scheme in [6]. We were able to prove linear stability and to derive conditions on the boundary operators. Here, we will repeat shortly the main ideas and results from [6] to further extend them to the non-linear case in the next section.

3.1 SATs in the Galerkin Framework

For simplicity reasons, we are considering the following linear advection equation


where and are initial and boundary data in and .
Instead of having an extra equation for the boundary condition as in (17), the condition is enforced weakly by some term which is called Simultaneous Approximation Term. One considers

In a Galerkin FE based discretization the solution is approximated by where are basis functions and are the coefficients. Let us further assume that are Lagrange polynomials where the degrees of freedoms are associated to points in the interval. Consider the variational formulation of (17) with test function and inserting the approximation yields




We consider now


If we set on both element boundaries degrees of freedom, then we obtain

The same holds if a quadrature rule is applied which is sufficiently exact. In this case, it can be shown that by imposing the boundary conditions weakly together with a suitable boundary operators the pure Galerkin scheme is energy stable, see [6, Proposition 3.4].
The result about stability can further generalized to universal FE semi-discretisations (9). For a general linear problem (scalar or systems) the formulation (9) can be written with penalty terms as


where is the boundary operators and the Theorem is formulated as follows:

Theorem 3.1 ([6], Theorem 3.5).

Applying a general FE semidiscretization (21) together with the SAT approach to a linear equation and let the mass matrix of (21) be symmetric. If the boundary operator together with the discretization can be chosen such that


has only non-positive eigenvalues

, then the scheme is (energy) stable.

Therefore, the determination of the operators is essential to guarantee stability. One is inspired by the continuous energy analysis as it can be seen in [6] and also in the following.

3.2 Estimation of the Boundary Operator

Here, we present now away to estimate the boundary operators using the continuous energy analysis. We consider as in [6] the symmetric linear hyperbolic system


where are the Jacobian matrices of the system, the matrix and the vector are known. is the number of boundary conditions to satisfy. are constant and the system (23

) is symmetrizable, i.e. it exists a symmetric and invertible matrix

such that for any vector the matrix

is symmetric with . Using the matrix , one can introduce new variables . The original variable can be expressed as and the original system (23) will become


Multiplying (24) form the left by we obtain the system


which is symmetric since . Focusing on the boundary treatment of the problem and using the weak formulation, we get for the system (23):


where is the boundary indicator function and is our boundary projection operator. The estimation of the boundary operator is driven by the energy balance for the weak formulation (26) in the continuous setting. We define the global energy of the solution of (23) by


where we take the symmetrizable of the system (23) under consideration. By multiplying (26) and integrating over , we obtain


We reformulate the left side of (28).


Combining (28) and (29), we get the following energy balance:


For stability, the energy does not increase, i.e. which is guaranteed if the integral term of the left-hand side of (30) is non-negative. Therefore,


has to be fulfilled. Inequality (31) imposes the restrictions on the choice of the projection operator . Switching to the discrete Finite Element setting, the information of (31) is used to determine and lead to stable pure Galerkin schemes as it is considered in [6].

4 Investigation of the Boundary Operator - Entropy Stability

As we have seen in section 3 the key to guarantee energy stability using the SAT approach is a proper definition of the boundary operators. Up-to-now exact integration is always considered in this context. Therefore, the continuous setting is applied to estimate the boundary operators for linear problems using the energy method and the results transform directly from the continuous to the discrete framework.
Here, we focus on the more general approach and describe how to estimate the boundary operators for non-linear hyperbolic systems.
As we have seen in 3, one key is the symetrizability of the system as it is already described in [17]. Therefore, the energy approach can be adapted as it is described in [35]. As it was already named before, the symetrizability of a system is already guaranteed if a entropy function exists [23] but then, we can further extend the estimations considering in general non-linear conservation laws (1).
Here, we give now -up to our knowledge- a first estimation on the boundary operators to extend stability results.
Extension of the linear case which is explained before for systems. Here, we consider the nonlinear case. Therefore, we have again . For the sake of simplicity, we are considering homogeneous boundary condition and we do further assume that the boundary matrix

is the identity matrix. Therefore, the boundary conditions reads


which means nothing else that the incoming waves are set to zero.
Later we restrict ourself to two dimensions but for now on we are still in . We have the following conservation law


where . If we impose the boundary conditions in the weak form by modifying the right side of the system (1) gets


where is the boundary indicator function and is the boundary projection operator which will be specified.
In section 2.3 we already introduced the concept of entropy. Assuming that we have an entropy pair where is a convex entropy function and is an entropy flux for (1). A symmetrization of (1) will be accomplished by introducing the new variable in place of . is called entropy variable and is given by . The mapping mapping is one-to-one and instead of working with in (33) we are working with by setting . In abuse of notation we won’t introduce a new notation for this change of the variable, but we can rewrite (33) in terms of the entropy variable and also using this to determine the boundary projection operator . Therefore, we get a dependence in of on the right side of (33). Also by following [23, 36], we know the following relations are fulfilled where we write everything in terms of the entropy variable . Now, let us consider the weak formulation (33) in terms of the entropy variable and multiply directly by and integrate in space. We obtain

Using (10), we get


Using Gauß-Green, we can re-write the second term and get


For stability (35), we have to determine in a way that the left side is not positive or

Therefore, we have a closer look on . This is a scalar function and by assuming smoothness we can apply the mean value theorem333In abuse of notation we use for the derivative of and do not apply everywhere the normal vector . Similar also for the relation (10).. It is


Using (10) by

yields in (36)


To estimate we apply (37) and get


is a constant factor whereas we have to calculate the rest depending on the . The term can be calculated either with a numerical quadrature formula or exact. We will demonstrate both calculations in section 5. Nevertheless, two questions automatically rise:

  • What is the connection between this estimation and the classical SAT approach for linear problems?

  • How does this estimation affect the stability properties of our continuous Galerkin method?

The first question can be easily answered. This approach directly simplifies to the linear case as it is described before using for the entropy variable .
The second question is more difficult. As we mentioned in section 3 the exactness of the quadrature rule is essential to transform the estimation from the continuous to the discrete setting. For a non-linear flux functions this is impossible and we get because of this directly stability issues even for scalar equations. Here, the correction term introduced in section 2.3 comes into play. It balance on every inner degree of freedom the entropy error. It allows to have exact calculation of this part like exact integrations are performed and hence the only remaining term is the boundary term. However, it further can rescue also in the linear case the scheme if for example the quadrature is not exact enough. We will present one example also for this case in the following section.
However, the correction term guarantees that in the discrete setting the entropy condition is fulfilled and we can prove the following:

Theorem 4.1.

A continuous Galerkin semi-discretization (9) for a hyperbolic conservation law (1) together with the entropy correction term (16) and the SAT approach with the boundary operator (37) is entropy stable.


The proof is a combination of the above considered calculation together with the properties of the correction term (16). Enforcing the boundary condition weakly using the SAT approach yields a pde (33). The semi-discretization with a continuous Galerkin approach leads to some stability issues due to the entropy production of the space discretization. However, the correction term is constructed in away to delete/ cancel out these. Multiplying with the approximation of the entropy variable and using the above manipulation with the boundary operations gives us finally an inequality (14) in the semi-discrete setting. This shows 4.1. ∎

Let us consider a first example how to construct the boundary operator (37).

Example 4.2.

We are considering the scalar Burgers’ equation in this context. The flux is given by . The entropy pair is and the entropy variable . Therefore, we do not have to change anything and we can directly estimate via . and so,


at the boundary . We will see the result (39) even more often in some other examples later.

Remark 4.3.

Entropy stability is guaranteed in the semidiscrete setting. However, to construct fully discrete entropy stable continuous Galerkin schemes three different ways can be applied but not further investigated here.

  • Using the SBP-SAT technique in time [31, 16] lead to an implicit time integration scheme where the stability results transform directly from the semidiscrete setting to the fully discrete ones.

  • To obtain an explicit method the relaxation technique of Ketcheson [24] can be applied where after each time step a non-linear equation has to be solved.

  • Finally, by using modal filters/artificial viscosity in an adaptive way as presented in [20], fully discrete entropy stable continuous Galerkin can be constructed.

Remark 4.4.

Different from the investigation in this section, there exists also other possibilities to determine the boundary operators. Hence, most if not all them comes with some restriction on the problem or the flux function. If one considers for example a multi-dimensional scalar case, another approach would be the use of a Taylor series expansion in respect to where sufficient smoothness is assumed.
Another possibility is if one has some symmetric behavior of the flux function and can bring all the information to the boundary. This is the key in the investigations of [17, 13, 29] for the linear case and for the nonlinear case [30] where the incompressible Navier-Stokes equation is considered.
Finally, also an recursive relation for the flux function can be found and used to determine the boundary operators. However, these approaches are not part of this investigation.

5 Numerical Examples

We demonstrate that a pure Galerkin scheme is entropy stable if the boundary procedure is done via the SAT approach which is introduced in section 4, and the correction terms are used. We impose the boundary conditions weakly and use an adequate boundary operator. We use the estimation (38) and develop using exact integration and numerical quadrature rules in the experiments (nearly no differences can be seen). If is non-negative, we set the value to zero and for negative we apply this value.
As basis functions either Bernstein or Lagrange polynomials of different orders (second to fourth order) on triangular meshes will be applied. Nearly no difference can be seen by applying either of these two bases.  The time integration is done via a strong stability preserving Runge-Kutta (RK) method with the same order of accuracy as the space discretization, e.g. [22] for details. Our scheme is semi-discrete entropy stable. Nevertheless, if shocks appear oscillations are seen and stability issues can appear especially for higher orders. One can increase the parameters in the free penalty terms or add additional streamline or diffusion terms to the correction terms to increase the stability properties (up-to-now we have only used correction terms to guarantee entropy conservation).

Also approaches to build fully discrete entropy stable schemes as described in remark 4.3 can be used but all of this is not topic of the current investigation.

5.1 Linear Equations

In the first experiments, we will focus on the correction terms and the influence of these terms on the scheme. The term (16) works specifically on every degree of freedom and is constructed to force the entropy condition (15).

Stabilisation Term for Linear Advection

In the first test, we are studying the same problem as described in [6] and consider the linear advection equation


where is the advection speed and the domain. The initial condition is given by

It is a small bump with hight one and located at . We consider homogeneous boundary conditions we do further assume that the boundary Matrix is the identity matrix. The boundary conditions reads for which means nothing else that the incoming waves are set to zero.
The bump is moving to the right with speed and no vertical movement can be seen. By using a pure Galerkin scheme and the SAT boundary procedure it can be shown that the scheme is stable under the restriction that the used quadrature formula for the mass matrix is exact, c.f. [6]. If we lower the accuracy of the applied quadrature rule the mass matrix is not exact anymore (up to machine precision). We touch only the highest order and so the mass matrix in (9) is not exact whereas the discretization describes exact.
Even with a low CFL number we run into stability issues and the scheme crashes as it is shown in in Figure 2. In figure 1(c), the structure of the bump can still be seen, but, simultaneously, the minimum value is and the maximum value is around . By going further in time it is getting worth.

(a) Initial data
(b) 300 steps
(c) 350 steps
Figure 2: 4-th order scheme in space and time

By applying the correction term the inaccuracy of the mass matrix can be compensated. We run the same test with the correction terms. The scheme remain completely stable as it can be seen in Figure 3.

(a) 300 steps
(b) 350 steps
(c) 1700 steps
Figure 3: 4-th order scheme in space and time with correction

In this test, the correction term works as a classical stabilization term to compensate the inaccuracy of the mass matrix. The term works at every degrees of freedom and is only applied where it is needed.
So, in case of a wrong/imprecise implementation of the Galerkin scheme, the term can rescue the stability properties. Nevertheless, for linear problems it should not be needed if the boundary procedure is done via the SAT approach described here and more specific in the first part of this series [6], see Figure 4.

(a) Initial data
(b) 100 steps
(c) 150 steps
Figure 4: 4-th order scheme in space and time, CFL=0.3

In the next part, we analyze the total amount and the influence of the correction terms on an exact pure Galerkin scheme.

Influence of the Correction Terms

The second test is a linear rotation equation in two space dimensions given by


where is the unit disk in . For the boundary, outflow conditions will be considered and the estimations are done using the techniques from section 4.

For time integration, the third order strong stability preserving scheme SSPRK(3,3) will be used with CFL number and the correction in each step will be done. A pure continuous Galerkin scheme of third order is used and Bernstein polynomials are applied as basis functions. In this test, a small bump which is located around is moving around in a circle. The rotation will be finished at . In figure (4(a)) the initial state is presented where figure (4(b)) shows the used mesh.

(a) Initial data
(b) Mesh
Figure 5: 4-th order scheme in space and time

The mesh contains 3582 triangular elements. We calculate the change of the entropy


after one rotation (1 second). The value is given by . Using a finer grid or smaller CFL number will also lower the entropy production in time. Finally, we like to mention that we renounce on an error plot here since nearly no differences to the results presented in [6] can be seen.

Next, non-linear equations will be considered.

5.2 Burgers’ Equation

Here, we focus first on a non-linear problem. We study Burgers’ equation


We use the same mesh and chose also the boundary condition as before. The bump is now located in the third quadrant and moves diagonal up. We are considering a non-linear problem and a shock is formed in time which can be recognized in figures 6-8. The boundary operator is calculated using the estimation (39) where we set for non-negative the boundary operator to to zero and for negative we apply . This procedure guarantees the inequality 39. Note that we can also determinate the boundary operator directly calculating (37) using a quadrature rule. Indeed, for most of the quadrature rules this value will be exact.

(a) Initial data
(b) 140 steps
(c) 200 steps
(d) 336 steps
Figure 6: 2-th order scheme in space and time
(a) 100 steps
(b) 200 steps
(c) 316 steps
Figure 7: 3-th order scheme in space and time
(a) 100 steps
(b) 200 steps
(c) 316 steps
Figure 8: 4-th order scheme in space and time

Focusing now on the experiment, the CFL number is set to and in figures 6-8 the experiment is printed using several orders. In the smooth regime the fourth order scheme demonstrates the best behavior as it can be seen in figure 8 where accuracy is lost after the formation of the shock and oscillations can be seen. Also in the second and third order schemes 6- 7 these oscillations appear but not as strong as in the 4th order test case. To delete/handle these oscillations further techniques have to be applied but are not topic of these current study. Next, we consider a Burgers’ type of equation where the exact evaluation of the boundary operator using a quadrature formula is not possible anymore. We investigate the differences and also the change of the entropy.

5.3 Burgers’ Type of Equation

The problem is given by


where is the unit disk in . Again, outflow boundary conditions will be considered. In this test case, a small bump located in zero will move up to the left and a shock will appear after a finite time. First, we consider the total change of the entropy. As time integration we choose a SSPRK(3,3) with CFL number and the space space discretisation is again done by a pure continuous Galerkin scheme of third order with Bernstein polynomials. In Table 1, the total change of entropy using the correction term and the SAT procedure for the test using a grid with elements. Here, the operator is calculate via a quadrature formula of order five and for non-negative values, it is set to zero again. The calculations ends right before the shock is build. We recognize the total amount is negative and entropy stability is obtained. If we increase the accuracy by using a finer grid, the change of the entropy will be tend to zero quite fast up to machine precision.

Time Change of entropy
Table 1: Total entropy change of numerical solutions using a continuous Galerkin scheme for the Burgers’ type of equation (44).

Now, we compare the difference by applying a quadrature rule or exact integration to calculate for the estimation of the boundary operator.

Exact Integration

In Figure 9, a 4th order Galerkin scheme with correction term is used. The time integration is done via a strong stability preserving RK method of fourth order using 5 stages. To estimate we determine the operator using exact integration and apply for . The results are presented in figure 9.

(a) Initial Condition
(b) 100 steps
(c) 174 steps
Figure 9: 4-th order scheme in space and time

Quadrature Formula

Running the same experiments using a quadrature formula of order to evaluate leads the results plotted in figure 10

(a) Initial Condition
(b) 100 steps
(c) 174 steps
Figure 10: 4-th order scheme in space and time, CFL=0.1

Obviously, no difference can be seen and also the values are up to machine precision identically.

Remark 5.1.

We run this test with several different order for the Galerkin scheme (Lagrange/Bernstein basis) and also for the quadrature rule no major difference have be seen in the results. Therefore, we conclude that applying a quadrature formula to determine have no effect on the accuracy and stability of the method.

All of our experiments verify that we obtain semi-discrete entropy stability for a pure Galerkin scheme if we apply the correction term together with our here derived boundary procedure. In the linear case, the entropy correction term is not needed but can be applied for stabilization reasons. However, our approach to estimate the boundary operator includes the linear cases used in the SBP-SAT framework by setting the entropy variable . Therefore, our approach can be seen as a natural extension to the linear case.

6 Conclusion

In this paper, we extend the investigation from [6] and present a way to build entropy stable Galerkin schemes. By switching to the entropy variables, we generalize the SAT approach to nonlinear problems and develop new estimations for the boundary operators. Our approach is in accordance with the boundary operators from the linear case, derived and used in the SBP-SAT community (e.g. [35, 18, 32] ).
By applying the boundary operator together with the correction term developed in [3] we are able to build semi-discrete entropy stable Galerkin schemes. Numerical experiments supports our theoretical analysis and demonstrate further that nearly no difference can be seen in the results if the boundary operator is developed using a quadrature rule or exact integration.
Further research in this direction is the construction of fully entropy stable Galerkin schemes using SBP-SAT in time or relaxation Runge-Kutta methods. A more detailed analysis of the constructed boundary operator is also desirable. Finally, the influence of the entropy correction term for more advance problems should also be considered. In [7], a first study is already been done. However, this analysis has to be extended.


P.Ö. has been funded by the the SNF grant (Number 200021_175784) and by the UZH Postdoc grant. This research was initiated by a first visit of JN at UZH, and really started during ST postdoc at UZH. This postdoc was funded by an SNF grant 200021_153604. The Los Alamos unlimited release number is LA-UR-19-32411.


  • [1] R. Abgrall. Residual distribution schemes: current status and future trends. Computers & Fluids, 35(7):641–669, 2006.
  • [2] R. Abgrall. High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices. Journal of Scientific Computing, 73(2-3):461–494, 2017.
  • [3] R. Abgrall. A general framework to construct schemes satisfying additional conservation relations. application to entropy conservative and entropy dissipative schemes. Journal of Computational Physics, 372:640–666, 2018.
  • [4] R. Abgrall, P. Bacigaluppi, and S. Tokareva. A high-order nonconservative approach for hyperbolic equations in fluid dynamics. Computers & Fluids, 19(1-3), 2017.
  • [5] R. Abgrall, E. l. Meledo, and P. Oeffner. On the connection between residual distribution schemes and flux reconstruction. arXiv preprint arXiv:1807.01261, 2018.
  • [6] R. Abgrall, J. Nordström, P. Öffner, and S. Tokareva. Analysis of the sbp-sat stabilization for finite element methods. arXiv preprint arXiv:1912.08108, 2019. Submitted.
  • [7] R. Abgrall, P. Öffner, and H. Ranocha. Reinterpretation and extension of entropy correction terms for residual distribution and discontinuous galerkin schemes. arXiv preprint, 2019. Submitted.
  • [8] R. Abgrall and P. L. Roe. High order fluctuation schemes on triangular meshes. Journal of Scientific Computing, 19(1-3):3–36, 2003.
  • [9] E. Burman and P. Hansbo. Edge stabilization for galerkin approximations of convection–diffusion–reaction problems. Computer Methods in Applied Mechanics and Engineering, 193(15-16):1437–1453, 2004.
  • [10] T. Chen and C.-W. Shu. Review of entropy stable discontinuous galerkin methods for systems of conservation laws on unstructured simplex meshes1.
  • [11] T. Chen and C.-W. Shu. Entropy stable high order discontinuous galerkin methods with suitable quadrature rules for hyperbolic conservation laws. Journal of Computational Physics, 345:427–461, 2017.
  • [12] H. Deconinck, K. Sermeus, and R. Abgrall. Status of multidimensional upwind distribution schemes and applications in aeronautics. AIAA paper, 4079:2007, 2000.
  • [13] A. Ern and J.-L. Guermond. Discontinuous galerkin methods for friedrichs’ systems. i. general theory. SIAM Journal on Numerical Analysis, 44(2):753–778, 2006.
  • [14] D. C. D. R. Fernández, J. E. Hicken, and D. W. Zingg.

    Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations.

    Computers & Fluids, 95:171–196, 2014.
  • [15] T. C. Fisher and M. H. Carpenter. High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains. Technical Report NASA/TM-2013-217971, NASA, NASA Langley Research Center, Hampton, VA 23681-2199, United States, February 2013.
  • [16] L. Friedrich, G. Schnücke, A. R. Winters, D. C. D. R. Fernández, G. J. Gassner, and M. H. Carpenter. Entropy stable space–time discontinuous galerkin schemes with summation-by-parts property for hyperbolic conservation laws. Journal of Scientific Computing, 80(1):175–222, 2019.
  • [17] K. O. Friedrichs. Symmetric positive linear differential equations. Communications on Pure and Applied Mathematics, 11(3):333–418, 1958.
  • [18] G. J. Gassner.

    A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods.

    SIAM Journal on Scientific Computing, 35(3):A1233–A1253, 2013.
  • [19] G. J. Gassner, A. R. Winters, and D. A. Kopriva. A well balanced and entropy conservative discontinuous Galerkin spectral element method for the shallow water equations. Applied Mathematics and Computation, 272:291–308, 2016.
  • [20] J. Glaubitz, P. Öffner, H. Ranocha, and T. Sonar. Artificial viscosity for correction procedure via reconstruction using summation-by-parts operators. In XVI International Conference on Hyperbolic Problems: Theory, Numerics, Applications, pages 363–375. Springer, 2016.
  • [21] E. Godlewski and P.-A. Raviart. Hyperbolic systems of conservation laws. Ellipses, 1991.
  • [22] S. Gottlieb, D. I. Ketcheson, and C.-W. Shu. Strong stability preserving Runge-Kutta and multistep time discretizations. World Scientific, 2011.
  • [23] A. Harten. On the symmetric form of systems of conservation laws with entropy. Journal of computational physics, 49:151–164, 1983.
  • [24] D. I. Ketcheson. Relaxation runge-kutta methods: Conservation and stability for inner-product norms. arXiv preprint arXiv:1905.09847, 2019.
  • [25] P. D. Lax. Shock waves and entropy. In E. H. Zarantonello, editor, Contributions to Nonlinear Functional Analysis, pages 603–634. Mathematics Research Center, University of Madison, Academic Press, 1971.
  • [26] M. S. Mock. Systems of conservation laws of mixed type. Journal of Differential equations, 37(1):70–88, 1980.
  • [27] J. Nordström. Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation. Journal of Scientific Computing, 29(3):375–404, 2006.
  • [28] J. Nordström. A roadmap to well posed and stable problems in computational physics. Journal of Scientific Computing, pages 1–21, 2016.
  • [29] J. Nordström. A roadmap to well posed and stable problems in computational physics. Journal of Scientific Computing, 71(1):365–385, 2017.
  • [30] J. Nordström and C. La Cognata. Energy stable boundary conditions for the nonlinear incompressible navier–stokes equations. Mathematics of Computation, 2018.
  • [31] J. Nordström and T. Lundquist. Summation-by-parts in time. Journal of Computational Physics, 251:487–499, 2013.
  • [32] P. Öffner and H. Ranocha. Error boundedness of discontinuous Galerkin methods with variable coefficients. Journal of Scientific Computing, pages 1–36, 2019.
  • [33] H. Ranocha, P. Öffner, and T. Sonar. Summation-by-parts operators for correction procedure via reconstruction. Journal of Computational Physics, 311:299–328, 2016.
  • [34] M. Ricchiuto and R. Abgrall. Explicit runge–kutta residual distribution schemes for time dependent problems: second order case. Journal of Computational Physics, 229(16):5653–5691, 2010.
  • [35] M. Svärd and J. Nordström. Review of summation-by-parts schemes for initial-boundary-value problems. Journal of Computational Physics, 268:17–38, 2014.
  • [36] E. Tadmor. The numerical viscosity of entropy stable schemes for systems of conservation laws. I. Mathematics of Computation, 49(179):91–103, 1987.