A Posteriori Error Estimates for Hughes Stabilized SUPG Technique and Adaptive Refinement for a Convection- Diffusion Problem

The motive of the present work is to propose an adaptive numerical technique for singularly perturbed convection-diffusion problem in two dimensions. It has been observed that for small singular perturbation parameter, the problem under consideration displays sharp interior or boundary layers in the solution which cannot be captured by standard numerical techniques. In the present work, Hughes stabilization strategy along with the streamline upwind/Petrov-Galerkin (SUPG) method has been proposed to capture these boundary layers. Reliable a posteriori error estimates in energy norm on anisotropic meshes have been developed for the proposed scheme. But these estimates prove to be dependent on the singular perturbation parameter. Therefore, to overcome the difficulty of oscillations in the solution, an efficient adaptive mesh refinement algorithm has been proposed. Numerical experiments have been performed to test the efficiency of the proposed algorithm.


Introduction
Singularly perturbed problems occur frequently in various branches of applied science and engineering, e.g., fluid dynamics, aerodynamics, oceanography, quantum mechanics, chemical reactor theory, reaction-diffusion processes, and radiating flows. In general, it has been observed that singularly perturbed problems exhibit singularities as the singular perturbation parameter ε ⟶ 0. Therefore, it becomes essential to implement some robust numerical technique to capture these singularities. In literature, there exist various numerical techniques to handle these singularities. Adaptive mesh refinement techniques are one of such techniques. Very few researchers have proposed adaptive refinement strategies for singularly perturbed convection-diffusion problems. Generally, adaptive refinement techniques are based on two types of error estimates, namely, a priori error estimates and a posteriori error estimates. Nicaise [1] developed a posteriori residual error estimates for convection-diffusionreaction problems using some cell-centered finite volume methods. Based on a posteriori error estimates, the author proposed an adaptive algorithm. John [2,3] did numerical study of various a posteriori error estimates and indicators for convection-diffusion problems. On the basis of error estimates, the author proposed numerical solution of singularly perturbed convection-dominated problems on adaptive refined grid. Repin and Nicaise [4] derived a posteriori error estimates for linear convection-diffusion-reaction problems using functional arguments. Verfurth [5] derived a posteriori error estimates for convection-dominated stationary convection-diffusion equation using locally refined isotropic meshes. Zhao et al. [6] proposed adaptive numerical technique for convection-diffusion equations based on semirobust residual a posteriori error estimates for lower order nonconforming finite element approximations of streamline diffusion method.
From literature, we know that the classical finite element methods [7] fail to provide satisfactory results for small values of singular perturbation parameter, i.e., when ε ⟶ 0. It has also been observed that streamline upwind/Petrov-Galerkin (SUPG) method provides good approximate solution in the region where there is no sharp change in the solution but fails in the small subregions of sharp boundary layers. It has been observed that occurrence of these nonphysical oscillations in the region of sharp boundary layers in the discrete solution of SUPG method is based on the fact that this scheme is not monotonicity preserving. To overcome this difficulty, in the present work, we have proposed Hughes stabilization strategy [8] along with the SUPG method. It involves suitable addition of one more term which is multiple of a function in the direction where spurious oscillations were seen in approximate SUPG solution. This additional term is added on the lefthand side of SUPG discretization of convection-diffusion problem. The a posteriori error estimates have been derived for the proposed scheme. Based on these estimates, an anisotropic mesh refinement strategy has been proposed for singularly perturbed problems.
The outline of the paper is as follows.
In Section 2, the continuous problem under consideration and its streamline upwind Petrov-Galerkin finite element approximation have been presented. In Section 3, some auxiliary tools which are required for deriving reliable error bounds have been presented. In Section 4, we have discussed residual-based a posteriori error estimates and derived error bounds on anisotropic meshes. An adaptive refinement algorithm based on derived a posteriori error estimates has been proposed in Section 5. Section 6 deals with some numerical experiments which have been performed to analyze the robustness and efficiency of the proposed adaptive refinement strategy. In the last section, conclusion has been presented.
Throughout the paper, we assume that −ð1/2Þ∇·a For any open bounded subset K ⊂ Ω, let H 1 ðKÞ be the standard Sobolev space. Further, we define Let be energy norm on bounded subset K ⊂ Ω. The weak formulation of equations (1a), (1b) and (1c) is given by the following.
Find u ∈ H 1 ðΩÞ such that The existence and uniqueness of the solution of the above weak formulation (4) are guaranteed using Lax Milgram Lemma together with condition (1c). Let Γ h be the admissible and shape-regular triangulation of domain Ω consisting of triangles. Let L be any two-dimensional element with edge E. Let n L,E = ðn x , n y Þ be unit outward normal vector to L along E (see Figure 1). Fixing one of the two normal vectors, let n E be the normal vector for each edge E.
It has been observed that the solution of singularly perturbed problem displays boundary layers if the Peclet number, as discussed below, is large. Define local mesh Peclet number as where h K min is the minimal length of element K as defined in the next section.
Let V h = fv h ∈ H 1 ðΩÞ: v h j K ∈ P 1 ðKÞg, where P 1 ðKÞ is the space of all linear polynomials over the element K and Next, we discuss the SUPG method along with the Hughes stabilization technique for approximating the solution of problem (1a), (1b) and (1c). The SUPG method [9] for problem (1a), (1b) and (1c) is defined as follows. Find and σ is a nonnegative stabilization parameter. This additional term increases the robustness of SUPG method in the boundary layer region by controlling oscillations. Using Hughes stabilization technique to SUPG finite element method, equation (1a), (1b) and (1c) is discretized as follows.
Let ρ K be the stabilization parameter over each element K. Ross et al. [7] showed that the approximate solution u h obtained using SUPG finite element discretization exists and is unique provided stabilization parameter ρ K is small and satisfies where h K min is the minimal length of element K and the constant ν satisfies the inequality From inequality (12), it can easily be observed that ν = 0 for piecewise linear functions in V h 0 . Therefore, the above bounds reduce to 0 ≤ ρ K ≤ b 0 /2kbk −2 ∞,K . In order to simplify the calculations, we introduce the notation c ≲ d which means that there exists a positive constant A independent of c, d, Γ h , and ε such that c ≤ Ad. Further, we assume that Also, for any mesh function v h ∈ V h 0 , using (12) and scaling arguments, we can get Using energy norm def. (3), Thus, we have Again, from energy norm, we have Using (16) and (17), we get

Some Important Notations and Tools
Since singularly perturbed convection-diffusion problems exhibit sharp boundary layers when Peclet number becomes large or the singular perturbation parameter becomes smaller, in such situations, elements with large aspect ratio (anisotropic meshes) are recommended. In this section, we will discuss some important results on anisotropic meshes.

Notations.
Consider an arbitrary triangle K ∈ Γ h with Q 0 Q 1 as the longest edge (see Figure 2). Denote two orthogonal vectors q i,K with length h i,K = jq i,K j, i = 1, 2, where q 1,K is taken along the largest edge Q 0 Q 1 . From Figure 2, it can be verified that h 1,K ≥ h 2,K . Define h K min = h 2,K . These q i,K 's correspond to two anisotropic directions. Further, define an orthogonal matrix C K = ðq 1,K , q 2,K Þ ∈ ℝ 2×2 . Let α K be the scaling factor defined as We represent triangles by K or K ′ or K i and its edges by E. Further define its height over edge E as where jKj represents the area of triangle K. Let w E be the bounded domain formed by using two triangles having common edge E. Further, define w K to be the domain consisting of triangle K and its edge neighboring triangles.

Advances in Mathematical Physics
Let be mesh Peclet number on the domain w K where Pe K is defined in (7). For an interior edge Since the mesh considered is assumed to be shape-regular and admissible, along with these requirements, we take and the number of triangles with node y j is bounded uniformly.

Interpolation.
In order to obtain reliable error upper bounds, we define a suitable matching function [11,12] to measure the alignment of anisotropic mesh Γ h and anisotropic function.
Definition 1 (matching function). Let u ∈ H 1 ðΩÞ and Γ h ∈ F be the triangulation of Ω. We define M 1 : where C K ∈ ℝ 2×2 as defined earlier.
We can easily verify that M 1 ðu, Γ h Þ~1 for isotropic meshes. Similarly, it can easily be observed that M 1 ðu, Γ h Þ~1 for anisotropic meshes suitably aligned with anisotropic function u. Therefore, M 1 ðu, Γ h Þ ≈ C for anisotropic meshes.
To propose reliable error estimates in energy norm, we will use Clément interpolation operator R C [13] for u ∈ H 1 ð ΩÞ as standard Lagrange interpolation cannot be defined for these functions.

Lemma 2.
Let u ∈ H 1 0 ðΩÞ and α K be the scaling factor defined by (19). Then, the Clément interpolation operator Proof. The proof is discussed in [14].

Residual Error Estimates
In this section, firstly we discuss exact and approximate residuals. Further, we will develop reliable error upper bounds for Hughes stabilized SUPG finite element solution on anisotropic meshes. It is shown that the error bounds obtained depend on anisotropic interpolation.

Exact Residuals.
We define exact element residual R K and exact edge residual R E as where n E ⊥E ⊂ Ω \ ∂Ω is the unitary normal vector and n⊥E ⊂ ∂Ω N is the outer unitary normal vector.

Approximate Residuals.
Let Q be the approximation operator used to approximate the element residual and the face residual, i.e., where we have denoted (approximate) element residual by r K and the (approximate) face residual by r E . Since the finite element solution v h is linear, thus 4.3. Residual Error Estimator. Residual error estimator η K and the approximation term ζ K over triangle K are given by where α K is the scaling factor defined earlier and a 2 K = 3α 2 K . Further, we define global error estimators as  Advances in Mathematical Physics Next, we derive reliable upper error bounds on anisotropic meshes.
Theorem 1 (residual error estimation). Let v ∈ H 1 0 ðΩÞ be the exact solution and v h ∈ V h 0 be the finite element solution obtained by the proposed scheme. Then, the error in energy norm is bounded above globally by Proof. We know that Bðv, vÞ ≥ kjvjk 2 ∀v ∈ H 1 0 ðΩÞ.

Using this result, we get
where u = v − v h . Introducing Clément interpolation operator R C , we can write the bilinear form B(.,.) as Now, using the error equation and integration by parts, we have Using equation (33), the middle term of equation (32) can be written as Using Cauchy Schwarz inequality, we get Further, using Lemma 2, we get Therefore, the term Bðv − v h , u − R C uÞ is bounded above by From (18) and (19), we have Next, we will find the bounds on the second term Bðv − v h , R C uÞ of equation (32).
Using the standard Galerkin orthogonality condition and standard scaling results, the above equation reduces to 5 Advances in Mathematical Physics We know that for Clément operator [11] Thus, we have It may be noted that the effect of nonlinear term in the L ∞ norm will be bounded by that of the term kak ∞ as shown below, i.e., Since a, the convection coefficient, is assumed to be smooth in the domain under consideration, it is bounded above. Hence, the nonlinear term ka h k ∞,K is taken as bounded above by some constant and is absorbed in the constant term.
We know that Therefore, Since using equation (45) and equation (37) in equation (32), we get Using triangle inequalities,

Adaptive Refinement Algorithm
In this section, we propose an adaptive refinement strategy based on the a posteriori error estimates obtained in the last section. We propose the following adaptive refinement algorithm: (1) Discretize the domain using triangular elements. Triangulation is being carried out using red refinement (2) Solve the problem using the proposed scheme described in Section 2 (3) Over each element K, the residual error estimates η K have been calculated as defined in Section 4 (4) Mark the elements fK e i g M e i =1 satisfying η K e i > C max L′ η L ′ , where C is a user chosen constant from ð0, 1Þ, for refinement (5) Refine these marked elements using green refinement procedure 6 Advances in Mathematical Physics (6) Refine all those elements having hanging nodes also to avoid the discontinuity in the solution (7) Solve the problem again on the new adapted mesh (8) Repeat the process of grid refinement until the solution has been obtained up to a given desired accuracy

Numerical Results
In this section, numerical experiments have been carried out in order to test the efficiency and robustness of the proposed adaptive refinement technique based on the derived error estimates.
Example 4. Consider the following singularly perturbed convection-diffusion problem: The right-hand side function f is so chosen to satisfy the exact solution The solution of the above problem exhibits exponential boundary layers along the lines x = 1 and y = 1. For adaptive refinement, an anisotropic triangular mesh has been taken into consideration. In Figures 3(a) and 3(b), we present a portion of adaptive triangular mesh for ε = 2 −3 with different degrees of freedom. Figures 4(a) and 4(b) present adaptive refined meshes for ε = 2 −5 with different degrees of freedoms.
In Figures 5 and 6, the numerical solution obtained using the proposed refinement algorithm for different values of the singular perturbation parameter ε = 2 −3 and ε = 2 −5 has been plotted. It can be easily seen that even very sharp boundary layers have been efficiently captured using the proposed refinement algorithm. From the solutions, it can also be observed that the problem is very sensitive to the singular perturbation parameter, i.e., even for ε = 2 −5 , very sharp boundary layers appear in the solution. In Figure 7, energy norm errors for ε = 2 −5 have been presented. The behavior of effectivity index ψ = kjv − v h jk/ErrorEstimator which is used to measure reliability of the estimator is shown in Figure 8.

Conclusion
In the presented work, an adaptive numerical technique has been proposed for singularly perturbed convectiondiffusion problems in two dimensions. The singularly perturbed problem under consideration has been solved using Hughes stabilization technique under SUPG finite element framework. Anisotropic meshes have been considered for the domain discretization. Reliable a posteriori error estimates have been developed in energy norm on anisotropic meshes. Based on these a posteriori error estimates, an  7 Advances in Mathematical Physics adaptive mesh refinement strategy has been proposed. It has been depicted through numerical experiments that the proposed adaptive refinement strategy is very much efficient in capturing sharp boundary layers as the singular perturbation parameter ε approaches to 0.

Data Availability
There is no data were used to support this study    Advances in Mathematical Physics