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 section5 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
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
A linear combination of basis functions will be used to describe the numerical solution
As basis functions we are working either with Lagrange interpolation where the degrees of freedom are associated to points inor 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
Using the divergence theorem, we get
By choosing for any and inserting (3), we obtain a system of equations:
that in practice we compute using a quadrature rule:
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:
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
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.
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
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).
The scheme (9) is called strongly energy stable if
holds. The term is bounded for any finite and independent from , and the mesh.
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:
Let be an SBP operator with that fulfills
Let and , and . If with , and with then the scheme (9) is energy stable.
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
let us consider the variational formulation of (8) with test function and inserting the approximation yields
In matrix formulation (14) can be written
as it is also described in [nordstrom2006conservative]. Let us check (12). Therefore, we consider
If the boundaries are included in the set of degrees of freedom, then we obtain
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:
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:
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):
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
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.
We use the energy approach and multiply our discretization with instead of and add the transposed equation using . We obtain
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
a skew-symmetric formulation must be applied in the way as it is described in
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 forand 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.
4.1.1 The scalar case
Consider the initial boundary value problem
Without loss of generality, it is enough to consider the homogeneous boundary conditions and we consider the spatial operator
considered in the subspace of functions for which . This operator will be dissipative if . Using the Gauss-Green theorem, we obtain
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
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.
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
We denote by , the adjoint of , multiply with and obtain
We transpose (24) and add both equations together. This results in
The boundary terms (BT) correspond to with a properly chosen . The matrix
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:
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
The energy satisfies
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:
In order to enforce the boundary conditions weakly, we modify this relation by:
The operators and are selected in such a way that
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 .
If , then
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:
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
in the sequel.
Step 1: Strong Implementation
First, we consider again the strong implementation of the problem.
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
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
if the matrix in the bracket is negative semidefinite.
Step 2: Weak Implementation
Assume now that we have chosen an such that
The energy is given (remember )
Second, we consider everything at the boundary
Since , we focus only on the part. We define such that and get for the integrands
Collecting the terms, we obtain
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
We re-order :
We choose with scalar and get
We choose such that the matrix is negative semidefinite. With , it is
4.1.3 Extension to 2D Symmetric Systems
Next, we will extend our investigation to the general hyperbolic system
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 matrixsuch 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 (37) will become
Multiplying (38) 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 (37):
where is Dirac distribution on 555i.e. for any smooth enough, and is our boundary projection operator.
We reformulate the left-hand side of (42).