ANALYSIS OF A DG METHOD FOR SINGULARLY PERTURBED CONVECTION-DIFFUSION PROBLEMS

In this article, we studied a discontinuous Galerkin finite element method for convection-diffusion-reaction problems with singular perturbation. Our approach is highly flexible by allowing the use of discontinuous approxi-mating function on polytopal mesh without imposing extra conditions on the convection coefficient. A priori error estimate is devised in a suitable energy norm on general polytopal mesh. Numerical examples are provided.

It is well-known that, when convection dominates diffusion (i.e. ϵ ≪ 1, which is also referred to as the singular perturbation parameter in this case), the solution of the boundary value problem typically possesses layers, which are thin regions where the solution and/or its derivatives change rapidly. Standard numerical methods will produce solutions with nonphysical oscillations in this case, unless the computational mesh has a size of the magnitude of the layers. Numerical stabilization techniques, including fitted mesh methods and fitted operator methods, have been developed to resolve the difficulty; cf. the books [8,9,11] and the references therein. In particular, DG methods with interior penalty have been proved to be an effective approach for solving convection-diffusion problems [3,6,7].
In this paper, we are interested in using discontinuous Galerkin methods with interior penalty to solve convection-diffusion-reaction problems with singular perturbation. In [2], Ayuso and Marini apply weighted-residual approach to recover discontinuous Galerkin formulations for advection-diffusion-reaction problems with singular perturbation. An optimal error estimation of the order O((ϵ 1/2 + h 1/2 )h k ) has been obtain on triangular mesh under a set of assumptions for the convection coefficient b. Analysis of finite element method on polygonal mesh is a new trend in the field of numerical PDE, such as the works [1,5,10,12] and the references therein. The goal of this article is to analyze an interior penalty discontinuous Galerkin method for the problem (1.1)-(1.2) on polygonal mesh. Compare to [2], our analysis is simple and allows the use of general mesh such as polytopal mesh, hybrid mesh and mesh with hanging nodes. In addition, a set of strict assumptions for the convection coefficient in [2] is removed in our analysis.
The DG methods in [2] and the DG methods proposed in this paper share the same finite element space. For the DG formulations in [2], integration by parts is used for the convection term and as a consequence, derivative of trial function turns to derivative of the test function. In this paper, we handle the convection term differently without applying integration by parts on it. As a result, we can relax the restriction of the convection coefficient.
The rest of this article is organized as follows. In Section 2, the DG FEM will be introduced. Analysis of the DG method is found in Section 3. Numerical experiments are presented in Section 4 to support the theoretical results.

Discontinuous Galerkin Finite Element Schemes
In this paper, standard definitions and notations of Sobolev spaces are used. For a polyhedron D ⊆ Ω, we denote H s (D) = W s,2 (D) the Hilbertian Sobolev space of index s ≥ 0 defined on D. The associated inner product, norm, and seminorms in H s (D) are denoted by (·, ·) s,D , ∥ · ∥ s,D , and | · | s,D , respectively. When s = 0, H 0 (D) coincides with the space of square integrable functions L 2 (D). In this case, the subscript s is suppressed from the notations of norm, semi-norm, and inner product. Furthermore, the subscript D is also suppressed when D = Ω. Throughout this article, we use C for generic constants independent of ϵ, mesh size, and the solution to equation (1.1)-(1.2), which may not necessarily be the same at each occurrence. We will use plain and bold fonts for scalars and vectors, respectively.
Let T h be a shape regular quasi-uniform triangulation of the domain Ω which consists of polygons in two dimension or polyhedra in three dimension satisfying a set of conditions specified in [13]. Denote by E h the set of all edges or flat faces in T h , and let E 0 h = E h \∂Ω be the set of all interior edges or flat faces. For every element T ∈ T h , we denote by h T its diameter. The mesh size for For a given integer k ≥ 1, let V h be a discontinuous Galerkin finite element space associated with T h defined as where P k is the space of polynomials of total degree up to k. Let T 1 and T 2 be elements sharing a common edge e, for which n 1 and n 2 are the unit outward normal vectors with respect to T 1 and T 2 , respectively. We define the jump and the average of a scalar valued function v on e as We introduce the following inner products for broken Sobolev spaces Define For any T ∈ T h , let v o represent the value of v at the element adjacent to T and let v o = 0 on e if e ⊂ ∂Ω ∩ ∂T . We may introduce some forms on V h as follows:

Analysis
First for any v ∈ V h , we define It is easy to verify that ||| · ||| is a norm in V h . For any function φ ∈ H 1 (T ), the following trace inequality holds true, The following lemmas are useful to establish the coercivity of the bilinear form a(·, ·).

3)
Proof. It follows from the definition of integration by parts that which implies (3.2). The equation (3.3) is a direct consequence of (3.2).
Proof. It follows from the definitions of ∂ + T and ∂ − T , which completes the proof of the lemma.
The following coercivity result holds.

Lemma 3.3.
For v ∈ V h , then for α large enough, Therefore, the DG formulation (2.2) has a unique solution.
Proof. It is well known that for α large enough
We next establish a priori error estimation of the DG method. Let Q h be a element-wise defined L 2 projections such that for each element T ∈ T h , Q h :

2). Then for any
Proof. Using the Cauchy-Schwarz inequality, the trace inequality and the definition of Q h , we have It follows from the definition of a c (·, ·) that We will bound the three terms on the right hand side of the above equation.
The inverse inequality and the definition of Q h imply whereb is the average of b on each element T ∈ T h . Obviously, we have It follows from the definitions of ∂ + T and ∂ − T , Using the equation above with w = u − Q h u and trace inequality (3.1), we have Combining the above estimates with (3.9) gives (3.8). We complete the proof of the lemma.

Theorem 3.1. Let u h ∈ V h be the DG finite element solution of the problem (1.1)-(1.2) arising from (2.2). Then there exists a constant C such that
(3.10) Proof. Let u be the solution of the problem (1.1)-(1.2). Then it is obvious that Subtracting (2.2) from the equation above yields Adding and subtracting Q h u give Letting v = Q h u − u h := e h , and using (3.5), (3.7) and (3.8), we have Using the triangle inequality, we have proved the theorem.

Numerical example
In this section, we present numerical results for the DG formulation (2.2). Examples of different types are used to confirm numerically the theoretical estimates in Section 3. P k elements (k = 1, 2, 3, 4) on either uniform rectangular meshes, or uniform pentagonal (hybrid with rectangles) meshes, shown in Figure 1, are used in all the examples. Note that as Q k elements fit rectangular meshes, the rectangular meshes are still considered as general polygonal meshes for the P k elements.  The source term f is determined such that the exact solution of (1.1) is We computed approximate solutions to the boundary value problem for ϵ = 10 −3 and 10 −9 to confirm the error estimates of the DG scheme. Here, α = 10 is used for (2.2) in all computations. Numerical error e h = Q h u−u h is measured in both the L 2 norm ||·|| and the energy norm |||·|||. The corresponding convergence rates for linear, quadratic, and cubic elements have been collected in Tables 1-3, respectively. Each convergence rate r is computed under the presumption that one has convergence of order O(h r ). The optimal L 2 estimate of order O(h k+1 ) is observed. On the other hand, numerical errors measured in the energy norm converges at O(h k+1/2 ), which matches the estimate in Theorem 3.1.  ( This problem has two characteristic boundary layers with the width of O( √ ϵ) near boundaries y = 0 and y = 1. P 1 finite element space is used on rectangular meshes. The numerical solution behaves as expected when ϵ = 10 −4 , (cf. Table 4), where it converges before resolving the layer and after resolving the layer. The numerical solutions for the boundary value problem with ϵ = 10 −4 and ϵ = 10 −10 are plotted in Figure 2. We can see the first solution resolves the boundary layer while the second does not. When ϵ = 10 −10 , the strongly convection-dominated case, optimal convergence rates in both L 2 and energy norms are obtained as shown in Table 4.   Due to the discontinuities in boundary conditions, the solution of this problem is not in H 1 (Ω). Numerical oscillations occur near the interior layer caused by the joints of the conflicting Dirichlet boundary conditions. This phenomenon has been reported in the literature for DG methods; see, e.g., [2,4]. We use P 2 and P 4 elements over rectangular meshes. The numerical solutions on a mesh of 1,024 uniform squares are plotted in Figure 3. The solutions do not oscillate outside of the internal layer.   where ϵ = 10 −9 , b = (1/2 − y, x − 1/2) T , and c = 10 −3 .
The numerical solutions for P 1 and P 4 elements over a mesh of 1,024 uniform squares are plotted in Figure 4. The numerical results indicate that the discontinuous Galerkin discretization is stable, and accurate.