Comparison of SUPG with Bubble Stabilization Parameters and the Standard SUPG

Abstract

We study a streamline upwind Petrov-Galerkin method (SUPG) with bubble stabilization coefficients on quasiuniform triangular meshes. The new algorithm is a consistent Petrov-Galerkin method and shows similar numerical performances as the standard SUPG when the mesh Péclet number is greater than 1. Relationship between the new algorithm and the standard SUPG will be explored. Numerical experiments support these results.

1. Introduction

We consider the reaction-convection-diffusion problem in 2D −𝜀Δ𝑢+𝐛⋅∇𝑢+𝑐𝑢=𝑓,𝑥,𝑦∈Ω,𝑢=0,𝑥,𝑦∈𝜕Ω, (1) where Ω is a bounded open domain in ℝ², with the boundary 𝜕Ω, 𝐛∈(𝑊^∞₁(Ω))²,    𝑐∈𝐿^∞(Ω) are given functions, and 𝜀≪|𝐛| is a small positive parameter. For simplicity, we only consider the case that Ω=(0,1)². In the following, suppose that there is a constant 𝜇₀>0 such that 1𝑐−2∇⋅𝐛≥𝜇₀>0,∀𝑥,𝑦∈Ω. (2) It is also assumed that 𝑓 is sufficiently smooth. In the case of 𝐛=(𝑏₁,𝑏₂)^𝑇>(0,0)^𝑇 (𝑏₁,𝑏₂≥𝛽, where 𝛽 is a positive constant), the solution of (1) typically has two exponential layers of width 𝑂(𝜀ln(1/𝜀)) at the sides 𝑥=1 and 𝑦=1 of Ω. In the case of 𝐛=(𝑏,0)^𝑇 (𝑏≥𝛽>0), the solution of (1) typically has an exponential layer of width 𝑂(𝜀) near the outflow boundary at 𝑥=1 and two characteristic (or parabolic) layers of width 𝑂(√𝜀) near the characteristic boundaries at 𝑦=0 and 𝑦=1.

When the mesh Péclet number Pe is greater than 1, there exist global unphysical oscillations in numerical solutions of standard discretization schemes on general meshes. Hence, stabilized methods and/or a priori adapted meshes are widely used in order to get discrete solutions with satisfactory accuracy. An overview on these methods can be found in the survey [1, 2].

One of the most famous stabilized finite element methods is the streamline upwind Petrov-Galerkin method (SUPG). The SUPG proposed by Hughes and Brooks [3] is known to provide good stability properties and high accuracy. However, there are several drawbacks in SUPG, such as lacking discrete maximum principle and involving second derivatives and difficulties in determining the stabilization coefficients. Driven by these problems, many researchers were devoted to improving the SUPG. A lot of numerical methods were proposed, such as SOLD [4], nonlinear residual [5], and LPS [6]. Also, relations of SUPG and other numerical formulation, like residual-free bubble method [7] and variational multiscale method [8], were studied to seek possible directions of improvement.

By means of residual, the SUPG adds to the original bilinear form a term which introduces a suitable amount of artificial viscosity in the direction of streamlines. Also, the SUPG can be viewed as an inconsistent Petrov-Galerkin method, since its modified weighting function cannot apply to the diffusive term (see the details in Section 4).

For this reason, we are to analyze the SUPG with bubble stabilization coefficients in 2D and compare its numerical performance with the SUPG’s. From theoretical analysis and numerical results, we find that the new scheme is classified into the consistent Petrov-Galerkin formulation (CPGF) and behaves as well as SUPG. Also, the standard finite element method (FEM), which shows excellent performances for Pe≤1, can be classified into the CPGF and viewed as a special “SUPG with bubble stabilization coefficients” by taking 𝛿_𝐾=0 in (7). Thus, the FEM and our new scheme could be viewed as two reference numerical methods in the CPGF. This provides possibilities of constructing new numerical schemes between them in the CPGF. In fact, in our forthcoming works, we obtain a linear maximum-principle-preserving stabilized method in the CPGF by means of the FEM and the SUPG with bubble stabilization coefficients, which shows better numerical performances than the standard SUPG. Thus, our results in this paper can be viewed as a starting point to construct new numerical schemes in the CPGF.

The remainder of this paper is organized as follows. In Section 2 we formulate the problem and introduce notation and mesh. Theoretical results including stability analysis and energy norm estimates can be found in Section 3. Section 4 is devoted to the relationship between SUPG, standard finite element method, and our method. Finally, numerical experiments that illustrate our theoretical results are presented in Section 5.

2. Mesh and Numerical Formulation

First we define a finite element space on triangular meshes 𝑉_ℎ∶=𝑣_ℎ∈𝐶Ω∶𝑣_ℎ𝜕Ω=0,𝑣_ℎ𝐾islinear,∀𝐾∈𝒯_ℎ, (3) where the term “linear” is to be understood in the usual isoparametric sense. Here we assume that the triangulations 𝒯_ℎ on Ω are quasiuniform: ℎ_𝐾>𝛼₀ℎ,𝛼_𝐾>𝛼₀, (4) for any 𝐾∈𝒯_ℎ, where 𝛼₀ is a positive constant and ℎ_𝐾,𝛼_𝐾,ℎ denote, respectively, the diameter of 𝐾, the smallest angle of 𝐾, and the maximum of all diameters of triangles in 𝒯_ℎ.

Using the linear finite element space 𝑉_ℎ, we can state the standard Galerkin discretisation of (1) which reads as follows.

Find 𝑢_ℎ∈𝑉_ℎ such that for all 𝑣_ℎ∈𝑉_ℎ, 𝑎_ℎ𝑢_ℎ,𝑣_ℎ=𝑓_ℎ𝑣_ℎ, (5) where 𝑎(𝑢,𝑣)=𝜀(∇𝑢,∇𝑣)+(𝐛⋅∇𝑢+𝑐𝑢,𝑣).

The SUPG consists in adding to the original bilinear form a term which introduces a suitable amount of artificial viscosity in the direction of streamlines. In this case, the SUPG reads as follows.

Find 𝑢_ℎ∈𝑉_ℎ such that for all 𝑣_ℎ∈𝑉_ℎ, 𝑎_ℎ𝑢_ℎ,𝑣_ℎ=𝑓_ℎ𝑣_ℎ, (6) where 𝑎_ℎ𝑢_ℎ,𝑣_ℎ=𝑎𝑢_ℎ,𝑣_ℎ+_𝐾−𝜀Δ𝑢_ℎ+𝐛⋅∇𝑢_ℎ+𝑐𝑢_ℎ,𝛿_𝐾𝐛⋅∇𝑣_ℎ𝐾,𝑓_ℎ𝑣_ℎ=𝑓,𝑣_ℎ+_𝐾𝑓,𝛿_𝐾𝐛⋅∇𝑣_ℎ𝐾. (7)

In (7) the term 𝛿_𝐾 is defined as 𝛿_𝐾(𝑥,𝑦)∶=𝛿_𝐾=𝐶^∗_𝐾ℎ𝜆₁𝜆₂𝜆₃,(𝑥,𝑦)∈𝐾 in which 𝐶^∗_𝐾 is a constant to be determined and 𝜆_𝑖   (𝑖=1,2,3) are the linear basis functions. Actually, 𝛿_𝐾 is a bubble function. Moreover, 𝜆₁𝜆₂𝜆₃≤127in𝐾,_𝐾𝜆₁𝜆₂𝜆₃=𝑆(𝐾),60 (8) where 𝑆(𝐾) represents the area of 𝐾.

Finally, we define a special energy norm (SD norm) associated with 𝑎_ℎ(⋅,⋅): 𝑢_ℎ²∶=𝜀𝑢_ℎ²₁+_𝐾𝛿_𝐾^1/2𝐛⋅∇𝑢_ℎ²_𝐾+𝑢_ℎ². (9)

We denote by ‖⋅‖_𝐷 the 𝐿² norm in 𝐿²(𝐷); that is, ‖𝑣‖²_𝐷=(𝑣,𝑣)_𝐷∀𝑣∈𝐿²(𝐷). (10) If 𝐷=Ω, then we drop Ω from the notation.

3. Stability and Energy Estimates

Throughout this subsection, we assume 𝐶 is some positive constant.

3.1. Stability Analysis

The stability properties are a consequence of the following.

Lemma 1. Let the parameter 𝐶^∗_𝐾 in 𝛿_𝐾 satisfies 𝐶^∗_𝐾≤27𝜇₀ℎ𝐶²_𝐾 (11) for each 𝐾∈𝒯_ℎ, where 𝐶_𝐾∶=‖𝑐‖_𝐿^∞_,𝐾. Then the discrete bilinear form is coercive; that is, 𝑎_ℎ𝑣_ℎ,𝑣_ℎ≥𝐴2𝑣_ℎ², (12) where 𝐴∶=max{𝜇₀,1}.

Proof. By divergence theorem we obtain 𝑎_ℎ𝑣_ℎ,𝑣_ℎ=𝜀𝑣_ℎ²₁+_Ω1𝑐−2∇⋅𝐛𝑣²_ℎ+_𝐾𝛿_𝐾^1/2𝐛⋅∇𝑣_ℎ²_𝐾+_𝐾𝐾𝑐𝑣_ℎ𝛿_𝐾𝐛⋅∇𝑣_ℎ. (13) Obviously the sum of first three terms is greater than 𝐴|||𝑣_ℎ|||² from the condition of (2), and we only need to estimate the last term. Using Hölder inequality, we have _𝐾𝐾𝑐𝑣_ℎ𝛿_𝐾𝐛⋅∇𝑣_ℎ≤_𝐾𝐾𝑐²𝑣²_ℎ𝛿_𝐾^1/2_𝐾𝐛⋅∇𝑣_ℎ²𝛿_𝐾^1/2≤12_𝐾𝐾𝑐²𝑣²_ℎ𝛿_𝐾+_𝐾𝐛⋅∇𝑣_ℎ²𝛿_𝐾≤𝜇₀2𝑣_ℎ²+12_𝐾𝛿_𝐾^1/2𝐛⋅∇𝑣_ℎ²_𝐾≤𝐴2𝑣_ℎ², (14) where (14) is based on the definition of 𝛿_𝐾, (8), and (11): _𝐾𝐾𝑐²𝑣²_ℎ𝛿_𝐾≤𝐶²_𝐾𝐾𝐾𝐶^∗_𝐾ℎ𝜆₁𝜆₂𝜆₃𝑣²_ℎ≤_𝐾𝐾𝜇₀𝑣²_ℎ=𝜇₀𝑣_ℎ². (15) Then the proof of the lemma is finished.

3.2. Energy Norm Estimate

We denote by 𝑢^𝐼 the nodal piecewise linear interpolant to 𝑢 over 𝒯_ℎ. From Lemma 1 and the fact of 𝑎_ℎ(𝑢−𝑢_ℎ,𝑢^𝐼−𝑢_ℎ)=0 one gets 𝜇₀2𝑢^𝐼−𝑢_ℎ²≤𝑎_ℎ𝑢^𝐼−𝑢_ℎ,𝑢^𝐼−𝑢_ℎ=𝑎_ℎ𝑢^𝐼−𝑢,𝑢^𝐼−𝑢_ℎ. (16) Denote 𝜂∶=𝑢^𝐼−𝑢,𝑒∶=𝑢^𝐼−𝑢_ℎ and estimate the right-hand side of (16) term by term: 𝜀∇𝜂,∇𝑒≤𝜀^1/2𝜂₁|||𝑒|||,=𝐛⋅∇𝜂+𝑐𝜂,𝑒_𝜕Ω𝐛⋅𝐧𝜂𝑒−_Ω∇⋅𝐛𝜂𝑒−_Ω(𝐛⋅∇𝑒)𝜂+_Ω=𝑐𝜂𝑒_Ω(𝑐−∇⋅𝐛)𝜂𝑒−_Ω(𝐛⋅∇𝑒)𝜂≤𝐶_Ω𝜂^21/2_Ω𝑒^21/2+_Ω(𝐛⋅∇𝑒)^21/2_Ω𝜂^21/2≤𝐶ℎ²|||𝑒|||+𝐶ℎ^−1/2|||𝑒|||ℎ²≤𝐶ℎ^3/2|||𝑒|||, (17)(18) where we have used the standard interpolation results ‖𝜂‖≤𝐶ℎ² (see [9]) and the first inequality of (18) is obtained by _Ω(𝐛⋅∇𝑒)²≤𝐶𝛽_Ω(∇𝑒)²=𝐶ℎ⁻¹𝛽(∇𝑒)²_𝐾𝐾𝛿_𝐾≤𝐶ℎ⁻¹_𝐾𝐾𝛿_𝐾(𝐛⋅∇𝑒)²≤𝐶ℎ⁻¹|||𝑒|||²,_𝐾𝜀Δ𝑢+𝐛⋅∇𝜂+𝑐𝜂,𝛿_𝐾𝐛⋅∇𝑒≤_𝐾𝐾𝐶𝜀+𝐛⋅∇𝜂+𝑐𝜂𝛿_𝐾^1/2𝛿_𝐾^1/2𝐛⋅∇𝑒≤_𝐾(𝐶𝜀+𝐛⋅∇𝜂+𝑐𝜂)𝛿_𝐾^1/2_𝐾𝛿_𝐾^1/2𝐛⋅∇𝑒_𝐾≤_𝐾𝛿_𝐾^1/2𝐶𝜀+𝐛⋅∇𝜂+𝑐𝜂²_𝐾^1/2_𝐾𝛿_𝐾^1/2𝐛⋅∇𝑒²_𝐾^1/2≤𝐶ℎ²|||𝑒|||. (19) Combining all of these estimates, one gets 𝑢−𝑢_ℎ≤𝐶ℎ^3/2. (20)

4. Comparison of SUPG with Bubble Stabilization Coefficients and the Standard SUPG

Consider the bilinear form of SUPG: 𝜀∇𝑢_ℎ,∇𝑣_ℎ+𝐛⋅∇𝑢_ℎ+𝑐𝑢_ℎ,𝑣_ℎ+_𝐾𝛿_𝐾−𝜀Δ𝑢_ℎ+𝐛⋅∇𝑢_ℎ+𝑐𝑢_ℎ,𝐛⋅∇𝑣_ℎ𝐾=𝑓,𝑣_ℎ+_𝐾𝑓,𝛿_𝐾𝐛⋅∇𝑣_ℎ𝐾, (21) where 𝛿_𝐾 is constant in 𝐾.

It can be rewritten in the form 𝜀∇𝑢_ℎ,∇𝑣_ℎ+𝐛⋅∇𝑢_ℎ+𝑐𝑢_ℎ,𝑣_ℎ+𝛿𝐛⋅∇𝑣_ℎ=𝑓,𝑣_ℎ+𝛿𝐛⋅∇𝑣_ℎ, (22) where 𝛿|_𝐾=𝛿_𝐾.

Notice that (22) does not correspond to a consistent Petrov-Galerkin formulation in general except in the case when 𝐛≡𝐂 and 𝑣_ℎ is linear. Clearly, when 𝛿_𝐾 is a bubble function, 𝜀∇𝑢_ℎ,∇𝛿𝐛⋅∇𝑣_ℎ=_𝐾𝜀∇𝑢_ℎ,∇𝛿_𝐾𝐛⋅∇𝑣_ℎ𝐾=𝜀_𝐾𝜕𝐾∇𝑢_ℎ⋅𝐧𝛿_𝐾𝐛⋅∇𝑣_ℎ−𝜀_𝐾𝐾∇⋅∇𝑢_ℎ𝛿_𝐾𝐛⋅∇𝑣_ℎ=𝜀_𝐾𝜕𝐾∇𝑢_ℎ⋅𝐧𝛿_𝐾𝐛⋅∇𝑣_ℎ≡0, (23) since 𝛿_𝐾 vanishes on the boundary of 𝐾. Thus (22) can be written as 𝜀∇𝑢_ℎ,∇𝑣_ℎ+𝛿𝐛⋅∇𝑣_ℎ+𝐛⋅∇𝑢_ℎ+𝑐𝑢_ℎ,𝑣_ℎ+𝛿𝐛⋅∇𝑣_ℎ=𝑓,𝑣_ℎ+𝛿𝐛⋅∇𝑣_ℎ. (24) Then SUPG with bubble stabilization coefficients is classified into the consistent Petrov-Galerkin formulation.

Moreover, in general, piecewise constants 𝛿_𝐾 in SUPG make the test function 𝑣_ℎ+𝛿_𝐾𝐛⋅∇𝑣_ℎ discontinuous. However, the test functions are continuous in the case of bubble stabilization coefficients and the consequent test space 𝑊_ℎ is contained in 𝐻¹₀(Ω). In this case, SUPG with bubble stabilization coefficients reads as follows.

Find 𝑢_ℎ∈𝑉_ℎ, such that for all 𝑤_ℎ∈𝑊_ℎ, 𝑎𝑢_ℎ,𝑤_ℎ=𝑓,𝑤_ℎ, (25) where 𝑊_ℎ∶=𝑤_ℎ∈𝐶Ω∶𝑤_ℎ𝐾=𝑣_ℎ𝐾+𝛿_𝐾𝐛⋅∇𝑣_ℎ𝐾,𝑣_ℎ∈𝑉_ℎ. (26)

On the other hand, SUPG with bubble stabilization coefficients is gradually close to the FEM in the same space 𝐻¹₀(Ω) as 𝐶^∗_𝐾→0 for any ∈𝒯_ℎ.

In a word, SUPG with bubble stabilization coefficients inherits the advantages of SUPG and constructs the relation between FEM and SUPG in consistent Petrov-Galerkin formulations.

5. Numerical Experiments

In this section we give numerical results that appear to support our theoretical results. Errors and convergence rates for our numerical scheme are presented. All calculations are carried out by using Intel visual Fortran 11. The discrete problems are solved by using a version of Pardiso solver (see [10, 11]).

For the computations we have chosen 𝛿_𝐾=60.0ℎ𝜆₁𝜆₂𝜆₃ in SUPG with bubble stabilization coefficients and 𝛿_𝐾=1.0ℎ in SUPG. We set Ω=[0,1]² and calculate the errors and convergence rates in the subdomain Ω_𝐶 away from layers for Tables 1–16. For parabolic Problems 1 and 2, Ω_𝐶=[0,1/2]×[1/3,2/3]. For exponential Problems 3 and 4, Ω_𝐶=[0,1/2]×[0,1/2]. The errors in Tables 17–20 are calculated in the whole domain Ω. In the following we only list the results of the case of 𝜀=10⁻² and 𝜀=10⁻⁸ since the comparison results between the two methods of other cases like 𝜀=10⁻⁴, 𝜀=10⁻⁶, and 𝜀=10⁻¹⁰ are similar.

Table 1 

𝜀=10⁻⁸, Problem 1 (comparison of error).

Table 2 

𝜀=10⁻⁸, Problem 1 (comparison of convergence rate).

Table 3 

𝜀=10⁻², Problem 1 (comparison of error).

Table 4 

𝜀=10⁻², Problem 1 (comparison of convergence rate).

Table 5 

𝜀=10⁻⁸, Problem 2 (comparison of error).

Table 6 

𝜀=10⁻⁸, Problem 2 (comparison of convergence rate).

Table 7 

𝜀=10⁻², Problem 2 (comparison of error).

Table 8 

𝜀=10⁻², Problem 2 (comparison of convergence rate).

Table 9 

𝜀=10⁻⁸, Problem 3 (comparison of error).

Table 10 

𝜀=10⁻⁸, Problem 3 (comparison of convergence rate).

Table 11 

𝜀=10⁻², Problem 3 (comparison of error).

Table 12 

𝜀=10⁻², Problem 3 (comparison of convergence rate).

Table 13 

𝜀=10⁻⁸, Problem 4 (comparison of error).

Table 14 

𝜀=10⁻⁸, Problem 4 (comparison of convergence rate).

Table 15 

𝜀=10⁻², Problem 4 (comparison of error).

Table 16 

𝜀=10⁻², Problem 4 (comparison of convergence rate).

Table 17 

𝜀=10⁻⁸, Problem 1 (comparison of global error).

Table 18 

𝜀=10⁻⁸, Problem 2 (comparison of global error).

Table 19 

𝜀=10⁻⁸, Problem 3 (comparison of global error).

Table 20 

𝜀=10⁻⁸, Problem 4 (comparison of global error).

Problem 1. One has −𝜀Δ𝑢+(2−𝑥)𝑢_𝑥+32𝑢=𝑓𝑥,𝑦∈Ω,𝑢=0𝑥,𝑦∈𝜕Ω, (27) where 𝑓 is chosen such that the solution 𝑢 is 𝑢=𝜋𝑥,𝑦sin2𝑒𝑥−^{−((1−𝑥)/𝜀)}−𝑒^−(1/𝜀)1−𝑒^−(1/𝜀)×1−𝑒^{√−(𝑦/𝜀)}1−𝑒^{√−((1−𝑦)/𝜀)}1−𝑒^√−1/𝜀. (28)

Problem 2. One has −𝜀Δ𝑢+𝑢_𝑥+𝑢=𝑥(1−𝑥)+𝑦1−𝑦,𝑥,𝑦∈Ω,𝑢=0𝑥,𝑦∈𝜕Ω. (29)

Problem 3. One has −𝜀Δ𝑢+2𝑢_𝑥+3𝑢_𝑦+𝑢=𝑓𝑥,𝑦∈Ω,𝑢=0𝑥,𝑦∈𝜕Ω, (30) where 𝑓 is chosen such that the solution 𝑢 is 𝑢()𝑥,𝑦=2sin𝑥1−exp−21−𝑥𝜀×𝑦²−1−exp1−𝑦𝜀. (31)

Problem 4. One has −𝜀Δ𝑢+2𝑢_𝑥+3𝑢_𝑦+𝑢=𝑥(1−𝑥)+𝑦1−𝑦𝑥,𝑦∈Ω,𝑢=0𝑥,𝑦∈𝜕Ω. (32)

From the above tables it is shown that the errors and convergence rates of SUPG with bubble stabilization coefficients and standard SUPG in the maximum norm, in the 𝐿² norm, and in the SD norm are similar not only in the subdomain away from layers but also in the global sense. These results illustrate that SUPG with bubble stabilization coefficients also has good stability properties and high accuracy as standard SUPG.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was partially supported by Promotive Research Fund for Excellent Young and Middle-aged Scientists of Shandong Province (Grant no. BS2012SF012) and Independent Innovation Foundation of Shandong University (Grant no. 2012TS194).