 # 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, stabilization terms are usually added and various formulations can be found in the literature. In this paper, we are also dealing with this problem, but present a different approach. We use the boundary conditions in our investigation in the sense that so called simultaneous approximation terms (SATs) are applied which are frequently used in the finite difference community. Here, the main idea is to impose the boundary conditions weakly. Specific boundary operators are constructed which guarantee stability. The SAT approach has already been used in the discontinuous Galerkin framework, but here we apply it – up to our knowledge – for the first time together with a continuous Galerkin scheme. We demonstrate that a pure continuous Galerkin scheme is stable if the boundary conditions are implemented in the correct way. This contradicts the general perception of stability problems for pure Galerkin schemes. In numerical simulations, we verify our theoretical analysis.

## Authors

##### 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 recent years, significant efforts have been made to construct and develop high-order methods for the solution of hyperbolic balance laws, and most of the common methods are either based on finite difference (FD) or finite element (FE) approaches. In the FE framework, one favorable, if not the most favorable scheme, seems to be the discontinuous Galerkin (DG) method introduced by Reed and Hill [reed1973triangular] because of its good stability properties [cockburn2012discontinuous, hesthaven2002nodal, chenreview, gassner2013skew]. In the stability proofs, so called summation-by-parts (SBP) operators are used [carpenter2014entropy, chen2017entropy, chan2018discretely, kopriva2014energy, gassner2013skew, ranocha2016summation]. SBP operators originate in the FD framework [kreiss1974finite] and lead to a way to demonstrate stability similar to the one in the continuous analysis [fernandez2014review, hicken2016multidimensional, svard2014review]. Together with SBP operators, Simultaneous Approximation Terms (SATs) that impose the boundary conditions weak are applied. The SBP-SAT technique is powerful and universally applicable. Certainly, one of the reasons for the popularity of DG is that the numerical solution is allowed to have discontinuities at the element boundaries, and, since non-linear hyperbolic problems are supporting shocks (like Riemann problems), this property is believed to be desirable. Another reason, maybe the most important, is that DG methods leads to block diagonal mass matrices which are easy to invert.

The difference between a DG approach and continuous Galerkin (CG), besides the structure of the mass matrix, is that in CG the approximated solution is forced to be continuous also over the element boundaries and this restriction seems to be quite strong also in terms of stability where a common belief in the research community is that a pure CG scheme is notoriously unstable111We like to mention that also parts of the authors had this belief before starting the project.. Therefore, stabilization terms have been developed and are frequently applied to remedy this issue [abgrall2019high, burman2010explicit, burman2004edge]. Even there exist some preliminary stability results [layton1983stable, layton1983stable1, gunzburger1977stability] including the procedure at the boundary where the main idea is to switch the norm of the trial space. However, these results may seem forgotten in the community.

In this paper, we also focus on the stability property of a pure Galerkin scheme, but follow a different approach. Our preliminary idea/thought is: If one considers the DG method with one element, the method is stable through the investigation done in the literature mentioned above. There is nothing that says that the approximation space must be a broken polynomial space, the only thing that is needed is that the trial and test function must have some kind of regularity within the elements, so that the divergence theorem (or SBP techniques) can be applied. Continuity is enough. Hence, one can see a CG method as a DG one, with only one element (the union of the simplex) with an approximation space made of polynomials with continuity requirement between the simplex. Hence, what is the difference between these two approaches? The answer to this question points to the procedure at the boundary.

In the stability proofs, the usage of SATs is essential. However, up-to-our-knowledge SATs have never been used together with a pure CG scheme and this is the topic of this paper. We divide the paper as follows: In the second section, we shorty introduce the continuous Galerkin scheme which is used and investigated in the following. Next, we introduce and repeat the main idea of the SAT procedure from the FD framework and extend it to the Galerkin approach. We show that the determination of the boundary operators is essential. In section 4

, we focus on the eigenvalue analysis of the spatial operators and derive conditions from the continuous setting to build adequate boundary operators in the discrete framework. We give some recipes which will be used in section

5 to support our analysis in numerical experiments. Finally, we conclude and discuss future work.

## 2 Continuous Galerkin Scheme

In this section, we shortly introduce the pure continuous Galerkin scheme (CG) as it is also known in the literature [Hughes1, burman2010explicit, kreiss1974finite]. We are interested in the numerical approximation of a hyperbolic problem

 ∂U∂t+divf(U)=0 (1)

with suitable initial and boundary conditions. Later, we will focus on the boundary condition, but for the explanation of CG this is not important. 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. 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

 Vh:=⨁K{Uh∈L2(K),Uh|K∈Pk}∩C0(Ω). (2)

A linear combination of basis functions will be used to describe the numerical solution

 Uh(x)=∑K∈Ωh∑σ∈KUhσφσ|K(x),∀x∈Ω (3)

As basis functions we are working either with Lagrange interpolation where the degrees of freedom are associated to points in

or Bézier polynomials.
To start the discretisation, we apply a Galerkin approach and multiply with a test function and integrate over the domain. This gives

 ∫Ω(Vh)T∂U∂tdx+∫Ω(Vh)Tdivf(U)dx=0. (4)

Using the divergence theorem, we get

 ∫Ω(Vh)T∂U∂tdx−∫Ω(∇Vh)Tf(U)dx+∫∂Ω(Vh)Tf(U)⋅ndγ=0. (5)

By choosing for any and inserting (3), we obtain a system of equations:

 ∑K∈Ωh∑σ′∈K(∂Uhσ′(t)∂t∫Kφσ′(x)φσ(x)dx−∫K∇φσ′(x)f(Uh)dx)=0, (6)

that in practice we compute using a quadrature rule:

 ∑K∈Ωh∑σ′∈K(∂Uhσ′(t)∂t∮Kφσ′(x)φσ(x)dx−∮K∇φσ′(x)f(Uh)dx)=0,

where represents the quadrature rules for the volume and surface integrals.

In this paper, we are considering only linear problems, i.e. the flux is linear in , but may depend on the spatial coordinate. In all the numerical experiences, we will make the spatial dependency simple enough (i.e. typically polynomial in ), so that it will always be possible to find a standard quadrature formula that make the formula exact. In other words, in this paper we proceed such that (6) is always exactly reproduced, unless it is specified.

Using a matrix formulation, we obtain the classical FE framework:

 \matM∂∂tU––h+\matF=0 (7)

where

denotes the vector of degrees of freedom,

is the approximation of and is a mass matrix 222In the finite difference community is called norm matrix and is classically abbreviated with , c.f. [svard2014review, nordstrom2006conservative].. In case of continuous elements, this matrix is sparse but not block diagonal, contrary to the discontinuous Galerkin methods. It is well-known that the continuous Galerkin scheme suffers from stability issues. Therefore, it is common to add stabilization terms to the scheme as for example in [burman2004edge]. However, we follow a different approach in this paper and will renounce these classical stabilisation techniques. In order to do this, we need more known results from the literature, which we will briefly repeat here.

## 3 Weak Boundary Conditions

### 3.1 SATs in SBP-FD framework

Weak boundary conditions implemented using simultaneous approximation terms (SATs) was originally developed in the finite difference (FD) framework. Together with summation-by-parts (SBP) operators it provides a powerful tool for proofs of semidiscrete () stability of linear problems by the energy method, see [svard2014review, fernandez2014review, nordstrom2017roadmap] for details.

Here, we present a short introductory example of the SBP-SAT technique as it presented [nordstrom2006conservative, svard2014review]. Consider the linear advection equation

 ∂u∂t+a∂u∂x=0,0≤x≤1,t>0,u(x,0)=uin(x),u(x,t)=b(x,t) for inflow boundary (8)

where is the initial condition and is the boundary data in that is only define on the inflow part of . In other words, if , the is only set for , and if , this will be for only. A semi-discretisation of (8) is given in terms of SBP operators.

 ∂u––∂t+a\matDu––=\matM−1S––,t>0,u––(0)=u––in, (9)

with are the coefficients of and similarly for . The coefficients correspond to the degrees of freedom in the finite element setting and are used to express the numerical solution (3) in the grid points. is a symmetric positive definite mass matrix which approximates the usual scalar product. The term is a difference matrix. This is exemplified below as

 \matDu––≈∂∂xu and ||u––||2\matM:=u––T\matMu––≈∫10u2(x)dx. (10)

Instead of having an extra equation on the boundary like in (8), the boundary condition is enforced weakly by the term which is called the SAT. We will focus on this operator more precisely and demonstrate how it should be selected to guarantee stability for (9).

###### Definition 3.1.

The scheme (9) is called strongly energy stable if

 ||u––(t)||2\matM≤K(t)(||u––in||2\matM+maxt1∈[0,t]||b0(t1)||2) (11)

holds. The term is bounded for any finite and independent from , and the mesh.

###### Remark 3.2.

The definition 3.1 is formulated in terms of the initial value problem (8) where only one boundary term is fixed. Indeed, extensions are straightforward. If an additional forcing function is considered at the right hand side of (8), we have to include the maximum of this function in (11) in the spirit of , for details see [svard2014review].

As established in [carpenter1994time], we can prove the following:

###### Proposition 3.3.

Let be an SBP operator with that fulfills

 \matQ+\matQT=\matB=diag(−1,0,⋯,0,1). (12)

Let and , and . If with , and with then the scheme (9) is energy stable.

###### Proof.

Multiplying (9) with yields

 u––T\matM∂∂tu––+au––T\matM\matDu––=u––TS––. (13)

 ddt||u––||2\matM=u––T\matM∂∂tu––+∂∂t–uT\matMu––=−au––T(\matQ+\matQT)–u+2u––TS––.

Further, we obtain from (12)

 ddt||u––||2\matM=(au20+2a+τu0(u0−b0))−(au2N−2a−τuN(uN−bN)).

Looking at the cases and , if , we find

 ddt||u––||2\matM≤−a+τ2(1+2τ)b20+a−τ2(1+2τ)b2N≤0.

This shows that boundary operator can be chosen in such way that it guarantees stability for the SBP-SAT approximation of (8). Next, we will apply this technique in the Galerkin framework.

### 3.2 SATs in the Galerkin-Framework

Instead of working with SBP-FD framework we consider now for the approximation of (8), a Galerkin approach. In [carpenter2014entropy, gassner2013skew] it is presented that also the discontinuous Galerkin spectral element method satisfies a discrete summation-by-parts (SBP) property and can be interpreted as an SBP-SAT scheme with diagonal mass matrix. As we described already in section 2, the differences between the continuous and discontinuous Galerkin approach are the solution space (2) and structure of the mass matrix (7) which is not block diagonal in CG. However, the approach with SAT terms can still be used to ensure stability also in case of CG but one has to be precise, as we will explain in the following.

Let us step back to the proof of Proposition 3.3 and have a closer look on it. Essential in the proof is condition (12). Let us focus on this condition for a FE based discretisation of (8) as described also in [nordstrom2006conservative]. We approximate (8) now with where are basis functions and are the coefficients. Let us assume that are Lagrange polynomials where the degrees of freedoms are associated to points in the interval. Introducing the scalar product

 ⟨u,v⟩=∫Iu(x)v(x)dx,

let us consider the variational formulation of (8) with test function and inserting the approximation yields

 ⟨∂∂tuh(t,x),φi(x)⟩+⟨a∂∂xuh(t,x),φi(x)⟩=0,∀i=0,⋯,N

i.e.

 ∫IN∑j=0(∂∂tuhj(t))φj(x)φ(x)dx+a∫IN∑j=0uhj(t)(∂∂xφj(x))φi(x)dx=0.

and finally

 N∑j=0Mi,j(∂∂tuhj(t))+aN∑j=0Qi,juhj(t)=0 (14)

where

 Mi,j=∫Iφj(x)φi(x)dx and Qi,j=∫I(∂∂xφj(x))φi(x)dx. (15)

In matrix formulation (14) can be written

 \matM∂∂tu––+a\matQu––=0

as it is also described in [nordstrom2006conservative]. Let us check (12). Therefore, we consider

 Qi,j+QTi,j=∫I(∂∂xφj(x))φi(x)dx+(∂∂xφi(x))φj(x)dx=∫I∂∂x(φj(x)φi(x))dx=φi(x)φj(x)|10=φi(1)φj(1)−φi(0)φj(0)∀i,j=0,⋯,N. (16)

If the boundaries are included in the set of degrees of freedom, then we obtain

 φi(1)φj(1)−φi(0)φj(0)=⎧⎨⎩1%fori=j=N,−1 for i=j=0,0 elsewhere.

Up to this point exact integrals are considered but the same steps are valid if a quadrature rule is applied with sufficient accuracy, especially (16) has to hold exactly. If this is the case, we can apply the proof of proposition 3.3 and demonstrate:

###### Proposition 3.4.

If the Galerkin method described above together with a SAT approach is applied to (8). If with is used then the described scheme is energy stable. The weak formulation of the problem writes:

 ⟨∂∂tuh(t,x),φi(x)⟩+⟨a∂∂xuh(t,x),φi(x)⟩=τa+(uh(0,t)−bh0(t))φi(0)+τa−(u(1,t)−bN)φi(1)=0,∀i,⋯,N
###### Proof.

For simplicity, we consider the case only.

The SAT techniques adds a penalty term into the approximation (8) on the right side. We focus now on the energy. Therefore, we multiply also with instead of and rearrange the terms. We obtain for the semi-discretization (14):

 N∑i,j,=0Mi,j(∂∂tuhj(t))uhi(t)+aN∑i,j=0Qi,juhj(t)uhi(t)=2aτuh0(t)(uh0(t)−b0(t))

where we used the fact that is valid. By following the steps of the proof of proposition 3.3 and using the above considerations we get the final result. ∎

In the derivation above, we restricted ourselves to one-dimensional problems using Lagrange interpolations. Nevertheless, this shows that a continuous Galerkin method is stable if the boundary condition is enforced by a proper penalty term. For the general FE semi-discretization of (7) it is straightforward to the procedure. For a general linear problem (scalar or systems) the formulation (7)333Remember that the flux is linear with Jacobian . can be written with penalty terms as

 \matM∂∂tU––h+\matQ1\matAU––h=αΠ(U––h) (17)

where is the boundary operator which includes the boundary conditions. represents the flux and the spatial operator. The matrix is usually sparse and we can formulate the following.

###### Theorem 3.5.

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

 (αΠ−\matQ1\matA)+(αΠ−\matQ1\matA)T (18)

has only non-positive eigenvalues , then the scheme is energy stable.

###### Proof.

We use the energy approach and multiply our discretization with instead of and add the transposed equation using . We obtain

 ddt||U––h||2\matM=U––h,T((αΠ−\matQ1\matA)+(αΠ−\matQ1\matA)T)U––h=W–––h,T\matΔW–––h≤0.

###### Remark 3.6.

This theorem yields directly conditions for a FE based stable scheme. If (18) is not fulfilled, stabilization terms have to be added. However, no internal stabilization terms are necessary when (16) holds. For this result to hold, a number of requirements are needed. The distribution of the degrees of freedom should be suitable for the problem and the mesh and the quadrature rule must be chosen to guarantee exactness of all the calculations. In the numerical test, we will present an example of what happens if the quadrature rule is not sufficiently exact.

Furthermore, in case of a non-conservative formulation of the hyperbolic problem or in case of variable coefficients a skew-symmetric formulation must be applied in the way as it is described in

[nordstrom2006conservative, ranocha2017extended, offner2019error] . If the implementation of the continuous Galerkin method is done in such way that the condition (18) are fulfilled, then by applying the SAT technique the method is stable only through our boundary conditions. To our opinion, this is contrary to common belief about continuous Galerkin methods. The only stabilizing factor needed is a proper implementation of boundary conditions. For the linear scalar case, the proof is given in (3.4). In the following, we will extend this theory to more general cases.

## 4 Estimation of the SAT-Boundary Operator

As described before, a proper implementation of the boundary condition is essential for stability. Here, we give a recipe how these SAT boundary operators can be chosen to get a stable CG scheme for different types of problems.

First, we are considering a scalar equation in 2D and transfer the eigenvalue analysis for the spatial operator from the continuous to the discrete setting. Then, we extend our investigation to 1D systems. Using again the continuous setting, we develop estimations for

and transfer the results to the finite element framework. Finally, we extend the investigation to the system case in two dimensions which will also be used in the numerical section later.

### 4.1 Eigenvalue Analysis

In the following part, we derive the conditions on the boundary operators and perform an eigenvalue analysis to get well posedness in the continuous setting. Next, the results are transformed to the discrete framework to guarantee stability of the discrete scheme.

#### Continuous Setting

Consider the initial boundary value problem

 ∂∂tu+a∂∂xu+b∂∂yu =0 x∈Ω, t>0, (19) Bu =g x∈∂Ω, t>0, u =f x∈Ω, t=0.

Without loss of generality, it is enough to consider the homogeneous boundary conditions and we consider the spatial operator

 Du :=(a∂∂x+b∂∂y)u, (20)

considered in the subspace of functions for which . This operator will be dissipative if . Using the Gauss-Green theorem, we obtain

 (21)

The operator will be dissipative if Looking at this question amounts to looking at the spectral properties of the operator . The question rises: How do we guarantee this condition? This is the role of the boundary conditions, i.e. when , we need to impose . For outflow, i.e. we get and using this information, we directly obtain

 ⟨u,Du⟩=∫∂Ωoutanu2ds>0. (22)

We do not go further into details about well posedness, but we recommend [nordstrom2017roadmap, nordstrom2019energy, nordstrom2019spatial] for details. Now, we transfer our analysis to the discrete framework and imitate this behavior discretely.

#### Discrete Setting

We have to approximate the spatial operator and the boundary condition (B.C), i.e. which is approximated an operator of the form . The term approximates where is used to describe weakly. denotes the mass matrix444In the finite difference context this matrix is also called norm matrix. and the derivative matrix, in the context of SBP [svard2014review]. Here, the projection operator works at the boundary points.

Looking at the dissipative nature of amounts to study its spectrum. The eigenvalue problem is

 \matM−1(\matQ−αΠ)~u––=λ~u–– (23)

We denote by , the adjoint of , multiply with and obtain

 ~u––∗,T(\matQ−αΠ)~u––=λ~u––∗,T\matM~u––=λ||~u––||2M. (24)

We transpose (24) and add both equations together. This results in

 ~u––∗,T((\matQ+\matQT)−α(Π+ΠT))~u––:=BT=2Re(λ)||~u––||2M. (25)

The boundary terms (BT) correspond to with a properly chosen . The matrix

 (α(Π+ΠT)−(\matQ+\matQT))≥0

is positive semi-definite. However, the remaining terms in the boundary term makes , i.e. the eigenvalues for the spatial operator have a strictly positive real parts only.

Next, we estimates the boundary operators such that theorem 3.5 is fulfilled and a pure CG scheme is stable. We start with the continuous energy analysis and derive the estimate above. Afterwards, we translate the result to the discrete FE framework as done for the scalar one-dimensional case.

#### 4.1.2 Continuous Energy Analysis for 1D Systems

First, we extend our result to the 1D system case. The problem under consideration is given by:

 ∂U∂t+A∂U∂x=0, x∈(0,1), (26) L0(U)=g0(t), x=0, L1(U)=g1(t), x=1, U(x,0)=U0(x), x∈[0,1],

where is a symmetric matrix, and we assume and to be linear. If is the number of incoming characteristics at and the number of incoming characteristics at , we know that the rank of (resp ) is (resp ). The system (26) admits an energy: if we multiply by , we first get

 ∫10UT∂U∂tdx+∫10UTA∂U∂xdx=0.

The energy satisfies

 dEdt+U(1,t)TAnU(1,t)−U(0,t)TAnU(0,t)=0.

To understand the role of boundary conditions, we follow what is usually done for conservation laws, we consider the weak form of (26): let be a regular vector function in space and time. We multiply the equation by , integrate and get:

 ∫T0∫10φT∂U∂tdxdt−∫T0∫10∂φ∂xTAUdxdt+∫T0φ(1,t)TAnU(1,t)dt−∫T0φ(0,t)TAnU(0,t)dt=0

In order to enforce the boundary conditions weakly, we modify this relation by:

 ∫T0∫10φT∂U∂tdxdt−∫T0∫10∂φ∂xTAUdxdt+∫T0φ(1,t)TAnU(1,t)dt−∫T0φ(0,t)TAnU(0,t)dt=∫T0φ(1,t)TΠ1(L1(U)−g1)dt−∫T0φ(0,t)Π0(L0(U)−g0)dt. (27)

The operators and are selected in such a way that

1. For any , the image of the boundary operator is the same as the image of , i.e. there is no loss of boundary information, and the same applies for and .

2. If , then

 dEdt<0.

A solution to this problem is given by the following: First, let where are the eigenvalues of and

is the matrix which rows are the right eigenvectors of

. Secondly, we have and choose:

 (28)

where are the negative eigenvalues only. Here we have introduced the operator and which are just and written using characteristic variables, and not . We will write

 Π(L(U)−g)=(L0(U)−g0)δx=0+(L1(U)−g1)δx=1

in the sequel.

#### Step 1: Strong Implementation

First, we consider again the strong implementation of the problem.

 ddt||U||2=UTAU∣∣∣x=0−UTAU∣∣∣x=1

For simplicity, we look only at what occurs for . Let where are the eigenvalues of and is the matrix which rows are the right eigenvectors of . Then, we obtain

 UTAU=UTXΛXTU=(XTU)TΛ(XTU)=(W+W−)T(Λ+00Λ−)(W+W−) (29)

with are the ingoing waves and they have the size of the positive eigenvalues . Analogously, are the outgoing waves with size of . A general homogeneous boundary condition is , since with that form we get

 UTAU=(W−)T(Λ−+RTΛ+R)W−≤0 (30)

if the matrix in the bracket is negative semidefinite.

#### Step 2: Weak Implementation

Assume now that we have chosen an such that

 Λ−+RTΛ+R≤0. (31)

The energy is given (remember )

 ∫ΩUT∂U∂tdx+∫ΩUTA∂U∂xdx=∫∂ΩUTΠ(W+−RW−)dx (32)

Second, we consider everything at the boundary

 12ddt||U||2=−12∫∂ΩUTAUdx+∫∂ΩUTΠ(W+−RW−)dx+∫∂Ω(W+−RW−)TΠTUdx (33)

Since , we focus only on the part. We define such that and get for the integrands

 W+Λ+W++W−Λ−W−+(W+)T~Π(W+−RW−)+(W+−RW−)T~ΠTW+ (34)

Collecting the terms, we obtain

 (W+W−)T(Λ++~Π+~ΠT−~ΠR−RT~ΠTΛ−)=:WB(W+W−) (35)

We must select such that the matrix is negative definite. Now, let us use the fact that we have the strong condition (31). By adding and subtracting, we obtain

 (W+W−)T(Λ++~Π+~ΠT−~ΠR−RT~ΠT−RTΛ+R)(W+W−)=:Q+W−(Λ−+RTΛ+R)W−≤0 by (???). (36)

We re-order :

 Q=(W+RW−)T(Λ++~Π+~ΠT−~Π−~ΠT−Λ+)(W+RW−)

We choose with scalar and get

 Q=(W+RW−)T(Λ+(1+2α)−Λ+−Λ+−Λ+)(W+RW−)=(W+RW−)T((1+2α)−1−1−1)⊗Λ+>0(W+RW−)

We choose such that the matrix is negative semidefinite. With , it is

 Q=−(W+RW−)T(1111)=:G⊗Λ+(W+RW−)

has the eigenvalues and and we obtain stability thanks to (36) and (31).

#### 4.1.3 Extension to 2D Symmetric Systems

Next, we will extend our investigation to the general hyperbolic system

 ∂U∂t+A∂U∂x+B∂U∂y =0, (x,y)∈Ω,t>0 (37) LnU =Gn (x,y)∈∂Ω,t>0

where are the Jacobian matrices of the system, the matrix and the vector are known, is the local outward unit vector, is the number of boundary conditions to satisfy. We assume that are constant and the system (37

) is symmetrizable. It exists a symmetric and invertible matrix

such that for any vector the matrix

 Bn=AnP

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

 P1/2∂V∂t+AP1/2∂V∂x+BP1/2∂V∂y=0. (38)

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

 ∂V∂t+P−1/2AP1/2∂V∂x+P−1/2BP1/2∂V∂y=0. (39)

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

 ∂U∂t+A∂U∂x+B∂U∂y=Πn(LnU−Gn)δ∂Ω, (40)

where is Dirac distribution on 555i.e. for any smooth enough, and is our boundary projection operator.

#### Energy balance

Again, we start by considering the energy balance for the weak formulation (40) in the continuous setting. We define the global energy of the solution of (37) by

 E=12∫ΩVTUdΩ, (41)

where we take into account the symmetrizability of the system (37). We multiply (40) and integrate over leads to

 ∫ΩVT(∂U∂t+A∂U∂x+B∂U∂y)dΩ=∫∂ΩVTΠn(Ln−Gn)dγ. (42)

We reformulate the left-hand side of (42).

 ∫ΩVT(∂U∂t+A∂U∂x+B∂U∂y)dΩ=∫ΩVT∂U∂tdΩ+∫ΩV