A modified diffusion coefficient technique for the ... - Semantic Scholar
Applied Mathematics and Computation 219 (2013) 9317–9330
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A modified diffusion coefficient technique for the convection diffusion equation S.A. Mohamed a,⇑, N.A. Mohamed a, A.F. Abdel Gawad b, M.S. Matbuly a a b
Department of Engineering Mathematics, Faculty of Engineering, Zagazig University, Egypt Department of Mechanical Power Engineering, Faculty of Engineering, Zagazig University, Egypt
a r t i c l e
i n f o
Keywords: Convection–diffusion Singularly-perturbed equations Multigrid Modified diffusion coefficient
a b s t r a c t A new modified diffusion coefficient (MDC) technique for solving convection diffusion equation is proposed. The Galerkin finite-element discretization process is applied on the modified equation rather than the original one. For a class of one-dimensional convection–diffusion equations, we derive the modified diffusion coefficient analytically as a function of the equation coefficients and mesh size, then, prove that the discrete solution of this method coincides with the exact solution of the original equation for every mesh size and/or equation coefficients. The application of the derived analytic formula of MDC is extended for other classes of convection–diffusion equations, where the analytic formula is computed locally within each element according to the mesh size and the values of associated coefficients in each direction. The numerical results of the proposed approach for two-dimensional, variable coefficients, with boundary layers, convection-dominated problems show stable and accurate solutions even on coarse grids. Accordingly, multigrid based solvers retain their efficient convergence rates. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction The convection–diffusion equations play an important role in many engineering and physics phenomena. In most practical problems, the magnitude of convection coefficient is much greater than that of diffusion coefficient. So, these problems are called convection-dominated or singularly-perturbed. The numerical solution of these problems represents serious difficulties because the solution of these diffusion–convection problems possess boundary layers that are small sub-regions, where derivatives of the solution are very large. These boundary layers make standard finite-element or finite-difference methods unsuitable for solving these problems. This is because the numerical solutions produce non-physical oscillations and low-order of accuracy unless refined meshes are introduced in the boundary regions using an adaptive mesh-refinement strategy. For this strategy to be effective, it is important that the error does not propagate into regions where refinement is not needed. So, the computational costs increase to obtain satisfied numerical results. Another difficulty occurs when multigrid is used for solving convection-dominated problems using classical discretization methods. Even if the grid where the solution is computed provides suitable accuracy, the multigrid algorithm requires a sequence of coarser grids, and it is important that the discretizations on these grids capture the character of the solution with a reasonable degree of accuracy. For these reasons, it is necessary to have a discretization strategy that does not have the deficiencies exhibited by the classical discretization method.
⇑ Corresponding author. E-mail address: [email protected] (S.A. Mohamed). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.03.014
9318
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
These difficulties have been the researchers’ main motivation, during the last two decades, to obtain discrete formulations that are, at the same time, accurate for smooth problems and stable for problems with boundary layers. High-order compact finite-difference schemes have been formulated including polynomial [1,2] or exponential [3,4] schemes. The high-order compact difference techniques have also been generalized to three-dimensions [5,6]. However, those high-order compact difference schemes are all based on uniform discretization. One existing strategy to deal with boundary layers is to first discretize the equation on a non-uniform (stretched) grid then a grid transformation technique is used to map the non-uniform grid to a uniform one, on which the fourth-order compact scheme with uniform mesh size is applied [7]. This strategy has been proven to be successful to recover fourth-order convergence rate in the presence of boundary layers. Two potential problems with the grid-transformation strategy are that the transformed equation becomes more complex, containing mixed terms in addition to the difficulty in constructing suitable grid-transformation functions for general applications. Another strategy is to use high-order compact difference scheme that can accommodate non-uniform mesh size [8,9]. Within the finite-element context, the same difficulties arise when the standard Galerkin method is applied to convection–diffusion equations. To overcome these difficulties, various formulations have been introduced including: upwinding techniques, see e.g. [10, Chapter 8] for a description of some of the most common of these schemes, streamline-upwind Petrov/Galerkin finite-element method (SUPG) [11,12], stream-line diffusion methods [13,14] and other special finite-element formulations, see e.g. [15,16]. As it is well-known, streamline-upwind Petrov/Galerkin method (SUPG) corresponds to modifying a weighting function in the Galerkin formulation to produce a small additional diffusion in the streamline direction. The amount of such additional diffusion is tuned by a parameter s that must be chosen in a suitable way. The choice of s is still considered as a major drawback of the method; a lot of numerical tests and several recipes have been proposed for the choice of s. The main contribution of the present work is a new technique for solving convection–diffusion equation by predicting a modified diffusion coefficient (MDC) such that the discretization process applies on the modified equation rather than the original one. For a class of one-dimensional convection–diffusion equation, we derive the modified diffusion coefficient analytically as a function of the equation coefficients and mesh size, then, prove that the discrete solution of this method coincides with the exact solution of the original equation for every mesh size and/or equation coefficients. Extending the same technique to obtain analytic MDC for other classes of convection–diffusion equations is not always straight forward especially for higher dimensions. However, we have extended the derived analytic formula of MDC (of the studied class) to general convection–diffusion problems. The analytic formula is computed locally within each element according to the mesh size and the values of the associated coefficients in each direction. The numerical results for twodimensional, variable coefficients, convection-dominated problems show that although the discrete solution does not coincide with the exact one, it provides stable and accurate solution even on coarse grids. As a result, multigrid based solvers benefit from these accurate coarse grid solutions and retain its efficiency when applied for convection–diffusion equations. Many numerical results are presented to investigate the convergence of classical algebraic and geometric multigrid solvers as well as Krylov-subspace methods preconditioned by multigrid. This paper is organized as follows; after the Introduction, Section 2 discusses the oscillatory behavior of classical finiteelement methods for convection-dominated convection–diffusion problems and interprets this property by the analysis and the numerical results in 1-D example. The modified diffusion method is introduced in Section 3 and is computed analytically for a class of 1-D convection–diffusion equation. In Section 4, numerical results are presented to demonstrate the validity of the computed modified coefficient and extend its application to 2-D convection–diffusion problems with variable coefficient. Finally in Section 5, we present the conclusions. 2. Performance of Galerkin method for convection–diffusion problems We consider steady convection–diffusion equation with suitable boundary conditions
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð1Þ
where X Rd is a smooth convex domain, dimension d = 2 or 3, w is the wind vector providing the convection coefficients, scalar e is the diffusion coefficient and f, g are the forcing function in X and the boundary function on C, respectively. As it is well-known, the classical Galerkin finite-element method fails when applied to convection-dominated, convection–diffusion problems. There are several ways to show this fact; we will give an explanation by considering a simple one-dimensional problem in the next subsection, and later on, a numerical experiment in 2-D is given in Section 4.2. 2.1. A class of convection–diffusion problem Consider the 1-D steady convection–diffusion equation
u00 þ wu0 ¼ 0 in Xð0; 1Þ uð0Þ ¼ 0; uð1Þ ¼ 1
Assuming that both
and w are constants and using m ¼ w , Eq. (2) becomes
ð2Þ
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
mu00 þ u0 ¼ 0 in Xð0; 1Þ
9319
ð3Þ
uð0Þ ¼ 0 uð1Þ ¼ 1 The exact solution of Eq. (3) is given by
u¼
ex=m 1 e1=m 1
ð4Þ
Discretiz ation of Eq. (3) with continuous, piecewise-linear finite-elements on a uniform grid, yields the scheme (
m
ui1 2ui þuiþ1 h
u0 ¼ 0;
þ
uiþ1 ui1 2
¼ 0;
i ¼ 1; . . . ; n 1
ð5Þ
un ¼ 1
where n is the number of elements, h = 1/n, xi = ih, and ui is the finite-element solution at node i. In this particular case, the scheme given by Eq. (5) is identical to the finite-difference central scheme. Introducing the mesh Peclet number defined as P e ¼ 2hm 8 into Eq. (5), the discrete algebraic system is given by:
> < u0 ¼ 0 ðPe 1Þui1 þ 2ui þ ðPe 1Þuiþ1 ¼ 0; > : un ¼ 1
i ¼ 1; . . . ; n 1
ð6Þ
Fig. 1 shows some characteristics of the general performance of the convection–diffusion equations where two values for
t = 0.1, t = 0.01 are considered. In this figure, both the discrete solution obtained by the classical Galerkin finite-element method (Eq. (6)) for different values of mesh sizes h = 0.1, h = 0.05, h = 0.01 and the exact solutions are plotted. It is easy to conclude the following observations which are well-known as the general performance of the convection– diffusion equations. As t ¼ w decreases, the convection term becomes dominant and a boundary layer appears where the solution becomes nearly singular. When t ¼ w ¼ 0:1, the numerical solution converges to the exact one even on a coarse grid. However, for a smaller value (t = 0.01), the discrete solution oscillates on coarse grids. On the coarsest grid (h = 0.1), the oscillations spread away from the boundary layer and the oscillations become local in the boundary layer when h = 0.05 and then converges to the exact solution as h becomes sufficiently small. Also, it can be noticed that accurate solutions can be obtained by classical Galerkin method on refined grids only when h is small enough such that the problem Peclet number Pe ¼ 2hm < 1. In the next section, we try to interpret this performance. 2.2. Interpretation of the oscillatory performance The exact solution of the algebraic system (6) can be written explicitly in terms of Pe, however, this solution is much simplified if we introduce the new parameter:
r¼
1 þ Pe 1 Pe
ð7Þ
It is easy to prove that the exact solution of (6) is
(
i
1r ui ¼ 1r n
i ¼ 0; . . . ; n
ui ¼ 0
i ¼ 0; . . . ; n 1; un ¼ 1 if Pe ¼ 1
if Pe–1
ð8Þ
It is important to notice that: r is undefined if Pe = 1 and the discrete solution is given by the second line in (8). for 0 < Pe < 1 ! r > 1; and hence the discrete solution is free of oscillations. Moreover, the discrete solution converges to the exact solution as h ? 0. This can be proven directly by taking the limit of ui as h ? 0
limui ¼ lim h!0
h!0
1 ri 1 r xi =h ¼ lim 1 rn h!0 1 r 1=h
But since,
limr h!0
1=h
1 þ 2hm ¼ lim h!0 1 h 2m
then,
limui ¼ h!0
1 exi =m 1 e1=m
!1=h
¼ lim 1 þ h!0
2h 2m h
1=h
¼ e1=m
9320
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
h=0.1 1
Solution u
0.5
Galerkin (Pe=0.5, r=3, exact (Pe=0.5, r=3, exact (Pe= 5, r=-1.5, Gaerkin (Pe=5, r=-1.5,
nu=0.1) nu=0.1) nu=0.01) nu=0.01)
0
-0.5
-1 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x h=0.05
Solution u
1
Galerkin (Pe=0.25, r = 1.667, Exact (Pe=0.25, r = 1.667, Galerkin (Pe=2.5, r = -2.333, Exact (Pe=2.5, r = -2.333,
0.5
nu=0.1) nu=0.1) nu=0.01) nu=0.01)
0
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.6
0.7
0.8
0.9
1
x h=0.01 1
Solution u
0.8
Galerkin (Pe=0.05, r = 1.1053, nu=0.1 ) Exact (Pe=0.05, r = 1.1053, nu=0.1) Galerkin (Pe=0.5, r = 3, nu=0.01) Exact (Pe=0.5, r = 3, nu=0.01)
0.6 0.4 0.2 0 0
0.1
0.2
0.3
0.4
0.5
x Fig. 1. Performance of Galerkin FEM, (¼ 2ht, r ¼ 1þPe ). 1Pe
which coincides with the exact solution (4) for P e > 1 ! r < 1; and hence the solution is accomplished with oscillations. 2.3. Investigation of numerical solution for high peclet number We are interested in the behavior of the discrete solution, Eq. (8), for h fixed and m going to zero (i.e., Pe going to +1). First note that if Pe ¼ 2hm 1; then 0 < 2hm 1 and hence r is approximated by: r ffi 1 – 4m/h. Then, substitution in Eq. (8) and neglecting higher orders of 2hm, lead to
9321
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330 4
1
x 10
1.2
n=10, nu=0.0000001, Pe=500000
n=11, nu=0.0000001, Pe=454550
1
0
0.8
-1
u
u
0.6
-2
0.4
-3
0.2
-4
0
-5
0
0.2
0.4
0.6
0.8
1
-0.2
0
x
0.2
0.4
0.6
0.8
1
x
Fig. 2. Oscillation of the discrete solution for large values of Pe for even and odd numbers of elements (n).
ui ffi
1 ð1 4m=hÞi 1 ð1Þi ½1 þ ð4m=hÞi n ffi 1 ð1 4m=hÞ 1 ð1Þn ½1 þ ð4m=hÞn
The asymptotic analysis of ui as m tends to zero leads to Eq. (9) where two cases: N even and N odd are distinguished.
(
if n is even; if n is odd;
ui ffi ui ffi
i n
if i is even
ðh=2mÞ 1n ni if i is odd ð2m=hÞi
if i is even
ð9Þ
1 ð2m=hÞðn iÞ if i is odd
Direct substitution of m = 0 in Eq. (9) shows that ui oscillates between 0 and 1 when n is odd. For even values of n, it is clear that ui is oscillating between a bounded value ðni Þ; and a negative unbounded value. In Fig. 2, the discrete solutions are plotted for m = 107 for even (n = 10) and odd (n = 11) numbers of intervals. As it is clear, the behavior is in full agreement with the previous analysis (see Eq. (9)). When n is even, the left figure, the discrete solution oscillates between a finite value at even nodes ði ¼ 0; 2; 4; . . .Þ, and negative very large values at odd nodes ði ¼ 1; 3; 5; . . .Þ. On the other hand, the right figure shows that the discrete solution oscillates between 0 at even nodes and 1 at odd nodes. 3. The modified diffusion coefficient method (MDC) In this section we prove that the discrete solution (Eq. (8)) of the class problem Eq. (3) coincides with the exact one (given by Eq. (4)) for every Pe if Galerkin finite-element is applied to Eq. (3) with its original diffusion coefficient m replaced by a modified coefficient m0. 3.1. Estimation of the modified diffusion parameter Let m be replaced by m + d in Eq. (3) then, the discrete system becomes
(
ðm þ dÞ u0 ¼ 0;
ui1 2ui þuiþ1 h
þ
uiþ1 ui1 2
¼ 0;
i ¼ 1; . . . ; n 1
un ¼ 1
ð10Þ
Introducing a parameter R, defined as
R¼
h 2ðm þ dÞ
ð11Þ
we obtain the discrete algebraic system
8 > < u0 ¼ 0 ðR 1Þui1 þ 2ui þ ðR 1Þuiþ1 ¼ 0; > : un ¼ 1
i ¼ 1; . . . ; n 1
ð12Þ
9322
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Eq. (12) is identical to Eq. (6) (just Pe is replaced by R). Now, define
r¼
1þR 1R
ð13Þ
Analogous to Eq. (8), the exact solution of the discrete system (12) is
(
i
1r ui ¼ 1r n
i ¼ 0; . . . ; n
if R – 1
ui ¼ 0 i ¼ 0; . . . ; n 1; un ¼ 1 if R ¼ 1 Substituting = xi/h, where h = 1/n, then for R – 1, the discrete solution at xi is given by
ui ¼
1 r xi =h 1 qxi ¼ ; 1 r 1=h 1q
q ¼ r1=h
ð14Þ
Since the exact solution Eq. (3) can be written as
u¼
ex=m 1 1 lx ¼ ; e1=m 1 1l
l ¼ e1=m
ð15Þ
It is easy to notice that discrete solution (14) coincides with the exact one (15) if r1/h = e1/m or equivalently, ¼ eh=m . Then, substitution from Eq. (11) yields 22mmþ2dþh ¼ eh=m . Solving for d yields þ2dh simplified as
ð1þR Þ1=h ¼ e1=m whichleadsto 1þR 1R h=m 1R mhÞð2mþhÞ d ¼ e ð22ð1e , which can be h=m Þ
d ¼ m½1 þ Pe cot hðP e Þ
ð16Þ
Thus, we have proven the following theorem. Theorem 3.1. The Galerkin finite-element solution of Eq. (3) with the diffusion coefficient m replaced by m0 = m + d = mPe cot h(Pe) produces the exact solution of Eq. (3) at every node in the mesh. It is important to notice that this theorem is independent of the mesh size and more importantly is independent on the u
2u þu
value of m. In fact, Theorem 3.1 proves that the addition of the term d i1 h i iþ1 to the left-hand side of Eq. (3) compensates the truncation error introduced by the numerical scheme. The relation between m and m0 is shown in Fig. 3 for the value of h = 1/50. It is noticed that, as m ? 0, m0 ? h/2, however, for greater values of m, the curves of both m0 and m coincide. 3.2. Generalization of the diffusion parameter First, a slight generalization of the class problem Eq. (2) is now considered as:
u00 þ wu0 ¼ 0 in Xða; bÞ uðaÞ ¼ ua ;
ð17Þ
uðbÞ ¼ ub
h/2=0.01
0.12 0.1
Modified nu nu
0.08
Modified nu
0.06 0.04 0.02 0
0
0.02
0.04
0.06
0.08
0.1
nu Fig. 3. Relation between the diffusion coefficient m and the modified one m0.
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
9323
Defining m = e/w, the exact solution of Eq. (17) is given by
u ¼ ua þ ðub ua Þ
ð1 eðxaÞ=m Þ ð1 eðbaÞ=m Þ
ð18Þ
To obtain the exact solution of the discrete algebraic system produced by applying classical finite-element, Eq. (17) is first xa a ¼ uuu x ¼ ba normalized to return to the form of Eq. (3). Let u and then Eq. (17) reduces to ua
(
b
u 00 þ u 0 ¼ 0 in Xð0; 1Þ m ð1Þ ¼ 1 uð0Þ ¼ 0; u
ð19Þ
2
1 ¼ w ba 0 ¼ ddux ; u 00 ¼ ddx2u. The exact solution of the algebraic system resulting by Galerkin finite-element on a uniwhere m and u formly divided domain can be obtained as
ui ¼ ua þ ðub ua Þ
1 r i 1 ðr 1=h Þxi a ¼ u þ ðu u Þ a a b 1 rn 1 ðr1=h Þba
ð20Þ
e where ui is the discrete solution at xi ¼ a þ ih; where h ¼ ðb aÞ=n and r ¼ 1þP , P e ¼ 2hm – 1. Using procedures similar to that 1P e
adopted in the previous section, it could be proven that the discrete finite-element solution will coincides on the exact solu is replaced in Eq. (19) by m 0 ¼ m Pe coth ½P e ; where, P e ¼ 2hm ; Pe ¼ 2hm. Then, m 0 ¼ 2h coth ð2hmÞ and we have the tion Eq. (18) if m following theorem. Theorem 3.2. The Galerkin finite-element solution of Eq. (17), with the diffusion coefficient m = e/w replaced by m0 ¼ ðb aÞmPe coth ðPe Þ; produces the exact solution of Eq. (17) at every node in the mesh. Other generalizations of Theorem 3.1 are: 1. On non-uniform grids, m0 is computed locally within each element e using the element mesh size he. Using Theorem 3.2, one can prove that exact solution is produced for the considered problem class Eq. (3) on non-uniform meshes. 2. For higher dimensional problems, m0 is computed differently in each direction according to the coefficients component in that direction. Although there is no guarantee that this process produces the exact solution, it is expected that it gives good approximate solutions. 3. For variable coefficients problems, m0 is computed locally within each element by approximating the coefficients by their corresponding values at the center of that element. 4. Numerical results In this section, we consider the linear and stationary diffusion–convection problem:
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð21Þ
where e is the diffusive coefficient, w is the transport advective field, f is the volume source term, X R2 is a smooth convex domain and g is the boundary value prescribed for u on C. 4.1. Problem 1 [15] The domain X is the unit square 0 6 x 6 1; 0 6 y 6 1, with mixed boundary conditions. More precisely, u = 0 on x = 0, u = 1 on x = 1 and @u ¼ 0 on y = 0 and y = 1. The given transport field is w = (w, 0) and f = 0. The importance of this problem @n is that it is of one-dimensional nature and is reducible to
u00 þ wu0 ¼ 0 in Xð0; 1Þ uð0Þ ¼ 0;
uð1Þ ¼ 1
which is identical to Eq. (2) that was analyzed in the previous sections. 1 As a first experiment, the domain of the problem is discretized uniformly with h ¼ 16 and the coefficients are chosen such 1 h that m ¼ w ¼ 50 producing Peclet number P e ¼ 2m ¼ 1:5625. To examine the proposed modified diffusion coefficient (MDC), the problem is solved directly by classical Galerkin finite-element, using stream diffusion (SD) method as implemented in [17] and using the proposed MDC method. In Fig. 4, the solution of the three methods are compared with the exact solution (Eq. (3)). As expected, the MDC solution coincides on the exact one while the Galerkin solution oscillates near the boundary region and SD, although prevents oscillations, fails to approximate the solution.
9324
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Nu=1/50, h=1/16, Pe=1.5625
1 0.8
Galerkin Exact Solution SD Proposed MDC
Solution u
0.6 0.4 0.2 0 -0.2 -0.4 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x Fig. 4. Comparison of solutions of problem 1 on a uniform grid by different numerical techniques.
Table 1 1 Maximum discretization error of problem 1 on uniform grid (h ¼ 16 ). Peclet number
Galerkin
SD
Proposed MSD
1.5625 2 3 5
0.2634 0.3516 0.5025 0.6693
0.0439 0.0183 0.0025 4.5400e-005
4.0870e-015 3.8129e-015 3.9166e-015 3.1866e-015
1
MDC Nu=1/160 Exact Nu=1/160 MDC Nu=1/48 Exact Nu=1/48 MDC Nu=1/320 Exact Nu=1/320
0.8
u
0.6 0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
x Fig. 5. Full agreement of the computed solution by MDC and the exact solution for different values of m on a non-uniform grid.
The previous experiment is repeated, but with different values of Pe and the results of maximum discretization error computed as Maxi juexact ui j, where ui is the computed solutions at node i are summarized in Table 1. The results of MDC are in full agreement with the analysis of the previous section. It produces the exact solution (up to round-off error 1015) for different values of Pe. In the previous results, a uniform grid was used. In the next experiment, a non-uniform grid is used where the domain 0 < x < 1 is subdivided into 16 intervals using a Shishkin grid [18]. Shishkin grids are specially designed to fit boundary layers where a transition line determines the boundary between the fine discretization that resolves the layer, and the coarse discretization outside of the layer. Using MDC on a piecewise uniform Shishkin grid, the results for different values of m are plotted in Fig. 5 showing that the computed discrete solution coincides on the curve representing the exact solution. On a nonuniform mesh, element Peclet number is computed locally within each element i as P ie ¼ 2hmi , where hi is the element length. Table 2 summarizes the results. As expected by the analysis of Section 3.2, MDC method produces the exact solution for the 1-D class of the convection– diffusion equation (Eq. (2)) even on non-uniform grids.
9325
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330 Table 2 Maximum discretization error on a non-uniform grid for problem 1.
m
P ie
Maxi juexact ui j
1/48 1/160 1/320
2.8336 9.4455 18.8910
2.9976e-015 3.0531e-015 3.2058e-015
4.2. Problem 2 [7] We next consider a constant coefficient convection–diffusion equation
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð22Þ
defined on the unit square 0 < x < 1; 0 < y < 1 with w = (1,0). The boundary condition is prescribed as
uð0; yÞ ¼ sin py;
uðx; 0Þ ¼ uðx; 1Þ ¼ 0;
uð1; yÞ ¼ 2 sin py
The exact solution is [19]
u¼
1 ex=2 sinh r
sin pyð2ex=2
sinh rx þ sinh rð1 xÞÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ p2 þ 1=42 . This problem represents a convection-dominated flow and was used as one of the test problems by Gupta et al. [19]. The coefficient of the convective term is constant. Fig. 6 shows the exact solution with = 0.05. For most part of the domain, the exact solution is smooth. But it has a steep boundary layer along the downstream edge at x = 1. Problem 2 is solved using classical Galerkin finite-element on a uniform grid (h = 1/16) for e = 0.001 and e = 0.05. The solutions are plotted in Fig. 7 showing nonphysical oscillations that are spread along the whole domain for the lower value of e = 0.001 and is localized near the boundary region as e is increased to 0.05. This behavior is similar to that shown in Fig. 1 for the 1-D case.
2 1 0 1 0.5 0
0.2
0
0.8
1
0.2
0.4
0.6
0.4
Fig. 6. Exact solution of problem 2 (e = 0.05).
2
2 1.5
1
1 0.5
0
0 1
-1 1
0.5
0.5
0
0
0.2
0.4
0.6
0.8
1
0
0
0.6
0.8
(b) Fig. 7. Galerkin finite-element solution of problem 2 with uniform grid size h = 1/16 for different values of e.
1
9326
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
2 1 0 1 0.5 0
0
0.6
0.4
0.2
0.8
1
Fig. 8. Finite-element solution using MDC of problem 2 ( = 0.001, h = 1/16).
x 10
Maximum Discretization error norm
-4
2 1 0 1
0.5
0
1
0.5
0
Fig. 9. Distribution of the discretization error norm of the finite-element solution using MDC of problem 2 ( = 0.001, h = 1/16).
Table 3 Convergence of MDC finite-element method for problem 2. h
1/16
1/32
1/64
1/128
Maximum discretization error
1.3084e-004
5.8178e-005
2.7241e-005
1.3175e-005
GMG Convergence using MDC, nu=1/200
0
10
h=1/64 h=1/32 h=1/16
-2
10
-4
Residual norm
10
-6
10
-8
10
-10
10
-12
10
-14
10
-16
10
0
2
4
6
8
10
12
14
16
18
20
V-cycles Fig. 10. Convergence of GMG for different h for problem 2 using MDC.
Next, problem 2 is solved on the same uniform grid using MDC method and the solution for = 0.001is plotted in Fig. 8. The absolute discretization error distribution is shown in Fig. 9. It is noticed that although extending MDC to this 2-D convection–diffusion problem does not produce the exact solution, it provides non-oscillatory accurate solution even on a coarse uniform grid. Table 3 reports the maximum discretization error norm of MDC finite-element method for different mesh sizes. As expected, multigrid methods do not converge if classical Galerkin finite-element is used to discretize the convection– diffusion problems since it fails to provide good approximations on coarser grids. However, introducing the MDC-Galerkin
9327
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
AMG Convergence using MDC, nu=1/200
2
10
0
h=1/16 h=1/32 h=1/64 h=1/128
10
-2
10
Residual norm
-4
10
-6
10
-8
10
-10
10
-12
10
-14
10
-16
10
2
4
6
8
10
12
14
16
18
20
V-Cycles Fig. 11. Convergence of AMG for different h for problem 2 using MDC.
discretization succeeded to retain the convergence of both geometric (GMG) and algebraic (AMG) multigrid methods as shown in Figs. 10 and 11, respectively. Convergence of AMG is slightly better than that of GMG but convergence rates are dependent on the mesh size h for both methods. 4.3. Problem 3 [7] In this test problem, we chose in Eq. (22) the following convection coefficients
w ¼ ðxðx 1Þð1 2yÞ; yðy 1Þð1 2xÞÞ It is obvious that this test problem has a stagnation point at (0.5, 0.5). A stagnation point (x0, y0) is the one inside the computational domain where both convection coefficients vanish. Convection–diffusion equations with stagnation points in their domains are usually used to model recirculation flow problems. For small values of e, this type of problem is very hard to solve, especially when standard multigrid method is used as the solution technique. The exact solution is chosen similar to the one used [7] as
u ¼ ð1 aðxÞÞð1 aðyÞÞ where,
aðxÞ ¼
eð1xÞ= 1 eð1yÞ= 1 ; aðyÞ ¼ 1= 1= e 1 e 1
Computed solution epslon=0.1
Computed solution epslon=0.01
0
0
-0.2
-0.2 -0.4
-0.4
-0.6
-0.6
-0.8
-0.8 -1 0
0 0.2
0.5 0.4
0.6
0.8
1
1
-1 0
0 0.5
0.5 1
Fig. 12. Computed solution on a stretched grid of problem 3 with b = 1/20, n = 32 and different values of
1
= 0.1, = 0.01.
9328
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
An exponential stretched grid in both x- and y-directions is used. The interval (0,1) is divided into n-sub-intervals by n + 1 iQ
n ; nodes such that xi ¼ 1e 1eQ
i ¼ 0; . . . ; n. As i changes from 0 to n, xi gives the coordinates of nodes along the x-axis from 0 to
firstinterval ln b 1. If the stretching parameter Q is chosen such that the ratio b ¼ lengthof , then Q ¼ nn1 . lengthof lastinterval
The computed solution of problem 3 for = 0.1, = 0.01 on a stretched grid is shown in Fig. 12 where the solution has steep boundary layers near x = 0 and y = 0. The width of these boundary layers becomes thinner as e decreases and sharp variation in the solution is observed. It must be mentioned that the exact solutions are graphically indistinguishable from
Table 4 Convergence of MDC finite-element method for problem 3. e
0.10
0.05
0.01
Maximum discretization error
3.654 104
5.431 104
3.821 103
10 10
Residual error norm
10 10 10 10 10 10
GMRES with different preconditioners
1
AMG GMG ILU non
0
-1
-2
-3
-4
-5
-6
0
5
10
15
20
25
30
iterations Fig. 13. Convergence behavior of GMRES for different preconditioners for problem 3.
AMG Convergene for problem 3 with different stretching ratios B 1
10
B=0.01 B=0.1 B=0.05
0
Residual Error Norm
10
-1
10
-2
10
-3
10
-4
10
-5
10
-6
10
2
4
6
8
10
12
14
16
18
20
AMG Cycles Fig. 14. Convergence behavior of AMG for different stretching factors b for problem 3.
9329
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
the shown computed solutions. It is clear that although the solution seems smooth for e = 0.1, it becomes more steep as e decreases and hence obtaining accurate solutions becomes more difficult. This is clear in Table 4 which reports the maximum discretization error for different values of e on a grid with mesh size h = 1/32. The resulting algebraic linear system produced by applying the finite-element method with the modified diffusion coefficient of problem 3 is solved by the generalized minimum residual GMRES method [20] with and without preconditioners. The results in Fig. 13 show the effect of preconditioning on the convergence of the GMRES. For this problem, the best performance is attained by the algebraic multigrid (AMG) preconditioner; only 7 GMRES iterations were sufficient to reduce the residual error norm below 106 of its initial value. The convergence of the GMRES with geometric multigrid (GMG) and incomplete LU (ILU) preconditioners are slightly slower but are significantly better than GMRES without any preconditioner. It is important to mention that classical multigrid does not converge when used to solve the resulting system. However, AMG works efficiently as observed from Fig. 14 where the convergence behavior is shown for different stretching ratios b. 4.4. Problem 4 [2] Although the previous two problems are singularly perturbed, they seem to have quasi-1D character. Therefore, in this Section we consider a problem that has different variability in each direction and hence represents a real 2D convection diffusion problem. This problem is similar to problem 3 but the boundary conditions and right hand side f(x, y) are computed such that the exact solution is u(x, y) = sin px + cos 3py + exy.
Computed solution h=1/16
Computed solution h=1/64 3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0 1
1 0.5
1
0.8
0.6
0.5
0.4
0.5
0.2
0 0
0
0
Fig. 15. Computed solution of problem 4 on a uniform grid with different mesh sizes h = 1/16, h = 1/64.
AMG convergence h=1/16
0
10
h=1/32 h=1/64 h=1/128
Maximum Error Norm
0 1
-5
10
-10
10
-15
10
2
4
6
8
10
12
14
AMG V-cycles Fig. 16. Convergence behavior of AMG for different mesh sizes.
9330
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Table 5 Convergence of MDC finite-element method for problem 4. h
1/16
1/32
1/64
1/128
Maximum discretization error Error reduction factor
5.11 104 –
1.27 104 4.02
3.16 105 4.02
7.91 106 3.99
The computed solutions using MCD on uniform grids with different mesh sizes h = 1/16, h = 1/64 are plotted in Fig. 15. The algebraic systems produced by MCD were solved efficiently by AMG whose convergence rates are presented in Fig. 16. The discretization error versus mesh size are given in Table 5 which show that a reduction of the discretization error by a factor of 4 is attained when the mesh size is reduced by a factor of 2. 5. Conclusions In this paper, we investigate the idea that the discretization error introduced by a numerical scheme to solve a boundary value problem can be compensated. A one-dimensional convection–diffusion equation with known exact solution is considered. The problem is discretized by a classical linear finite-element technique (equivalent to central finite-difference) and the analytic solution of the discrete system was computed and compared to the exact one to get the condition for being identical. We choose this condition in the form of a change in the normalized diffusion coefficient and get an analytic formula for the modified diffusion coefficient as a function of the original coefficients and the mesh size. Numerical results are presented to demonstrate that implementation of the proposed technique produces the exact solution for this class of 1-D singularly-perturbed problems even on coarse grids with uniform or non-uniform mesh sizes. As an extension to 2-D and variable coefficients convection–diffusion problems, the obtained formula is computed locally within each element and each space direction to approximate the solution. To demonstrate this approach and investigate the convergence behavior of the multigrid method, several numerical experiments are presented and accurate solutions were obtained even on coarser grids. As a result, multigrid-based solvers retain its efficient convergence rates. References [1] M.M. Gupta, R.P. Manohar, J.W. Stephenson, A single cell high order scheme for the convection–diffusion equation with variable coefficients, Int. J. Numer. Methods Fluids 4 (1984) 641–651. [2] S. Karaa, J. Zhang, Convergence and performance of iterative methods for solving variable coefficient convection–diffusion equation with a fourth-order compact difference scheme, Comput. Math. Appl. 44 (2002) 457–479. [3] G.Q. Chen, Z. Gao, Z.F. Yang, A perturbational, h4 exponential finite difference scheme for convection diffusion equation, J. Comput. Phys. 104 (1993) 129–139. [4] Z.F. Tian, S.Q. Dai, High-order compact exponential finite difference methods for convection diffusion type problems, J. Comput. Phys. 220 (2007) 952– 974. [5] M.M. Gupta, J. Zhang, High accuracy multigrid solution of the 3D convection–diffusion equation, Appl. Math. Comput. 113 (2000) 249–274. [6] Y. Wang, J. Zhang, Fast and robust sixth-order multigrid computation for the three-dimensional convection-diffusion equation, J. Comput. Appl. Math. 234 (2010) 3496–3506. [7] L. Ge, J. Zhang, High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids, J. Comput. Phys. 171 (2001) 560–578. [8] Y. Ge, F. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems, J. Comput. Phys. 230 (2011) 4051–4070. [9] J. Zhang, H. Sun, J.J. Zhao, High order compact scheme with multigrid local mesh refinement procedure for convection diffusion problems, Comput. Methods Appl. Mech. Eng. 191 (2002) 4661–4674. [10] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, Berlin, Heidelberg, 2008. [11] A. Russo, Streamline-upwind Petrov/Galerkin method (SUPG) vs residual-free bubbles (RFB), Comput. Methods Appl. Mech. Eng. 195 (2006) 1608– 1620. [12] E.G. Carmo, G.B. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the streamline and approximate upwind/Petrov–Galerkin method, Comput. Methods Appl. Mech. Eng. 192 (2003) 3379–3396. [13] N. Madden, M. Stynes, Linear enhancements of the streamline diffusion method for convection–diffusion problems, Comput. Math. Appl. 32 (10) (1996) 29–42. [14] N. Kopteva, How accurate is the streamline-diffusion FEM inside characteristic (boundary and interior) layers?, Comput Methods Appl. Mech. Eng. 193 (2004) 4875–4889. [15] M. Farhloul, A.S. Mounim, A mixed-hybrid finite element method for convection–diffusion problems, Appl. Math. Comput. 171 (2005) 1037–1047. [16] X.-G. Li, C.K. Chan, S. Wang, The finite element method with weighted basis functions for singularly perturbed convection–diffusion problems, J. Comput. Phys. 195 (2004) 773–789. [17] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics, Oxford, New York, 2005. [18] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [19] M.M. Gupta, R.P. Manohar, J.W. Stephenson, High-order difference schemes for two-dimensional elliptic equations, Numer. Methods Partial Differ. Equ. 1 (1985) 71. [20] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, 2003.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A modified diffusion coefficient technique for the convection diffusion equation S.A. Mohamed a,⇑, N.A. Mohamed a, A.F. Abdel Gawad b, M.S. Matbuly a a b
Department of Engineering Mathematics, Faculty of Engineering, Zagazig University, Egypt Department of Mechanical Power Engineering, Faculty of Engineering, Zagazig University, Egypt
a r t i c l e
i n f o
Keywords: Convection–diffusion Singularly-perturbed equations Multigrid Modified diffusion coefficient
a b s t r a c t A new modified diffusion coefficient (MDC) technique for solving convection diffusion equation is proposed. The Galerkin finite-element discretization process is applied on the modified equation rather than the original one. For a class of one-dimensional convection–diffusion equations, we derive the modified diffusion coefficient analytically as a function of the equation coefficients and mesh size, then, prove that the discrete solution of this method coincides with the exact solution of the original equation for every mesh size and/or equation coefficients. The application of the derived analytic formula of MDC is extended for other classes of convection–diffusion equations, where the analytic formula is computed locally within each element according to the mesh size and the values of associated coefficients in each direction. The numerical results of the proposed approach for two-dimensional, variable coefficients, with boundary layers, convection-dominated problems show stable and accurate solutions even on coarse grids. Accordingly, multigrid based solvers retain their efficient convergence rates. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction The convection–diffusion equations play an important role in many engineering and physics phenomena. In most practical problems, the magnitude of convection coefficient is much greater than that of diffusion coefficient. So, these problems are called convection-dominated or singularly-perturbed. The numerical solution of these problems represents serious difficulties because the solution of these diffusion–convection problems possess boundary layers that are small sub-regions, where derivatives of the solution are very large. These boundary layers make standard finite-element or finite-difference methods unsuitable for solving these problems. This is because the numerical solutions produce non-physical oscillations and low-order of accuracy unless refined meshes are introduced in the boundary regions using an adaptive mesh-refinement strategy. For this strategy to be effective, it is important that the error does not propagate into regions where refinement is not needed. So, the computational costs increase to obtain satisfied numerical results. Another difficulty occurs when multigrid is used for solving convection-dominated problems using classical discretization methods. Even if the grid where the solution is computed provides suitable accuracy, the multigrid algorithm requires a sequence of coarser grids, and it is important that the discretizations on these grids capture the character of the solution with a reasonable degree of accuracy. For these reasons, it is necessary to have a discretization strategy that does not have the deficiencies exhibited by the classical discretization method.
⇑ Corresponding author. E-mail address: [email protected] (S.A. Mohamed). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.03.014
9318
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
These difficulties have been the researchers’ main motivation, during the last two decades, to obtain discrete formulations that are, at the same time, accurate for smooth problems and stable for problems with boundary layers. High-order compact finite-difference schemes have been formulated including polynomial [1,2] or exponential [3,4] schemes. The high-order compact difference techniques have also been generalized to three-dimensions [5,6]. However, those high-order compact difference schemes are all based on uniform discretization. One existing strategy to deal with boundary layers is to first discretize the equation on a non-uniform (stretched) grid then a grid transformation technique is used to map the non-uniform grid to a uniform one, on which the fourth-order compact scheme with uniform mesh size is applied [7]. This strategy has been proven to be successful to recover fourth-order convergence rate in the presence of boundary layers. Two potential problems with the grid-transformation strategy are that the transformed equation becomes more complex, containing mixed terms in addition to the difficulty in constructing suitable grid-transformation functions for general applications. Another strategy is to use high-order compact difference scheme that can accommodate non-uniform mesh size [8,9]. Within the finite-element context, the same difficulties arise when the standard Galerkin method is applied to convection–diffusion equations. To overcome these difficulties, various formulations have been introduced including: upwinding techniques, see e.g. [10, Chapter 8] for a description of some of the most common of these schemes, streamline-upwind Petrov/Galerkin finite-element method (SUPG) [11,12], stream-line diffusion methods [13,14] and other special finite-element formulations, see e.g. [15,16]. As it is well-known, streamline-upwind Petrov/Galerkin method (SUPG) corresponds to modifying a weighting function in the Galerkin formulation to produce a small additional diffusion in the streamline direction. The amount of such additional diffusion is tuned by a parameter s that must be chosen in a suitable way. The choice of s is still considered as a major drawback of the method; a lot of numerical tests and several recipes have been proposed for the choice of s. The main contribution of the present work is a new technique for solving convection–diffusion equation by predicting a modified diffusion coefficient (MDC) such that the discretization process applies on the modified equation rather than the original one. For a class of one-dimensional convection–diffusion equation, we derive the modified diffusion coefficient analytically as a function of the equation coefficients and mesh size, then, prove that the discrete solution of this method coincides with the exact solution of the original equation for every mesh size and/or equation coefficients. Extending the same technique to obtain analytic MDC for other classes of convection–diffusion equations is not always straight forward especially for higher dimensions. However, we have extended the derived analytic formula of MDC (of the studied class) to general convection–diffusion problems. The analytic formula is computed locally within each element according to the mesh size and the values of the associated coefficients in each direction. The numerical results for twodimensional, variable coefficients, convection-dominated problems show that although the discrete solution does not coincide with the exact one, it provides stable and accurate solution even on coarse grids. As a result, multigrid based solvers benefit from these accurate coarse grid solutions and retain its efficiency when applied for convection–diffusion equations. Many numerical results are presented to investigate the convergence of classical algebraic and geometric multigrid solvers as well as Krylov-subspace methods preconditioned by multigrid. This paper is organized as follows; after the Introduction, Section 2 discusses the oscillatory behavior of classical finiteelement methods for convection-dominated convection–diffusion problems and interprets this property by the analysis and the numerical results in 1-D example. The modified diffusion method is introduced in Section 3 and is computed analytically for a class of 1-D convection–diffusion equation. In Section 4, numerical results are presented to demonstrate the validity of the computed modified coefficient and extend its application to 2-D convection–diffusion problems with variable coefficient. Finally in Section 5, we present the conclusions. 2. Performance of Galerkin method for convection–diffusion problems We consider steady convection–diffusion equation with suitable boundary conditions
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð1Þ
where X Rd is a smooth convex domain, dimension d = 2 or 3, w is the wind vector providing the convection coefficients, scalar e is the diffusion coefficient and f, g are the forcing function in X and the boundary function on C, respectively. As it is well-known, the classical Galerkin finite-element method fails when applied to convection-dominated, convection–diffusion problems. There are several ways to show this fact; we will give an explanation by considering a simple one-dimensional problem in the next subsection, and later on, a numerical experiment in 2-D is given in Section 4.2. 2.1. A class of convection–diffusion problem Consider the 1-D steady convection–diffusion equation
u00 þ wu0 ¼ 0 in Xð0; 1Þ uð0Þ ¼ 0; uð1Þ ¼ 1
Assuming that both
and w are constants and using m ¼ w , Eq. (2) becomes
ð2Þ
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
mu00 þ u0 ¼ 0 in Xð0; 1Þ
9319
ð3Þ
uð0Þ ¼ 0 uð1Þ ¼ 1 The exact solution of Eq. (3) is given by
u¼
ex=m 1 e1=m 1
ð4Þ
Discretiz ation of Eq. (3) with continuous, piecewise-linear finite-elements on a uniform grid, yields the scheme (
m
ui1 2ui þuiþ1 h
u0 ¼ 0;
þ
uiþ1 ui1 2
¼ 0;
i ¼ 1; . . . ; n 1
ð5Þ
un ¼ 1
where n is the number of elements, h = 1/n, xi = ih, and ui is the finite-element solution at node i. In this particular case, the scheme given by Eq. (5) is identical to the finite-difference central scheme. Introducing the mesh Peclet number defined as P e ¼ 2hm 8 into Eq. (5), the discrete algebraic system is given by:
> < u0 ¼ 0 ðPe 1Þui1 þ 2ui þ ðPe 1Þuiþ1 ¼ 0; > : un ¼ 1
i ¼ 1; . . . ; n 1
ð6Þ
Fig. 1 shows some characteristics of the general performance of the convection–diffusion equations where two values for
t = 0.1, t = 0.01 are considered. In this figure, both the discrete solution obtained by the classical Galerkin finite-element method (Eq. (6)) for different values of mesh sizes h = 0.1, h = 0.05, h = 0.01 and the exact solutions are plotted. It is easy to conclude the following observations which are well-known as the general performance of the convection– diffusion equations. As t ¼ w decreases, the convection term becomes dominant and a boundary layer appears where the solution becomes nearly singular. When t ¼ w ¼ 0:1, the numerical solution converges to the exact one even on a coarse grid. However, for a smaller value (t = 0.01), the discrete solution oscillates on coarse grids. On the coarsest grid (h = 0.1), the oscillations spread away from the boundary layer and the oscillations become local in the boundary layer when h = 0.05 and then converges to the exact solution as h becomes sufficiently small. Also, it can be noticed that accurate solutions can be obtained by classical Galerkin method on refined grids only when h is small enough such that the problem Peclet number Pe ¼ 2hm < 1. In the next section, we try to interpret this performance. 2.2. Interpretation of the oscillatory performance The exact solution of the algebraic system (6) can be written explicitly in terms of Pe, however, this solution is much simplified if we introduce the new parameter:
r¼
1 þ Pe 1 Pe
ð7Þ
It is easy to prove that the exact solution of (6) is
(
i
1r ui ¼ 1r n
i ¼ 0; . . . ; n
ui ¼ 0
i ¼ 0; . . . ; n 1; un ¼ 1 if Pe ¼ 1
if Pe–1
ð8Þ
It is important to notice that: r is undefined if Pe = 1 and the discrete solution is given by the second line in (8). for 0 < Pe < 1 ! r > 1; and hence the discrete solution is free of oscillations. Moreover, the discrete solution converges to the exact solution as h ? 0. This can be proven directly by taking the limit of ui as h ? 0
limui ¼ lim h!0
h!0
1 ri 1 r xi =h ¼ lim 1 rn h!0 1 r 1=h
But since,
limr h!0
1=h
1 þ 2hm ¼ lim h!0 1 h 2m
then,
limui ¼ h!0
1 exi =m 1 e1=m
!1=h
¼ lim 1 þ h!0
2h 2m h
1=h
¼ e1=m
9320
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
h=0.1 1
Solution u
0.5
Galerkin (Pe=0.5, r=3, exact (Pe=0.5, r=3, exact (Pe= 5, r=-1.5, Gaerkin (Pe=5, r=-1.5,
nu=0.1) nu=0.1) nu=0.01) nu=0.01)
0
-0.5
-1 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x h=0.05
Solution u
1
Galerkin (Pe=0.25, r = 1.667, Exact (Pe=0.25, r = 1.667, Galerkin (Pe=2.5, r = -2.333, Exact (Pe=2.5, r = -2.333,
0.5
nu=0.1) nu=0.1) nu=0.01) nu=0.01)
0
-0.5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0.6
0.7
0.8
0.9
1
x h=0.01 1
Solution u
0.8
Galerkin (Pe=0.05, r = 1.1053, nu=0.1 ) Exact (Pe=0.05, r = 1.1053, nu=0.1) Galerkin (Pe=0.5, r = 3, nu=0.01) Exact (Pe=0.5, r = 3, nu=0.01)
0.6 0.4 0.2 0 0
0.1
0.2
0.3
0.4
0.5
x Fig. 1. Performance of Galerkin FEM, (¼ 2ht, r ¼ 1þPe ). 1Pe
which coincides with the exact solution (4) for P e > 1 ! r < 1; and hence the solution is accomplished with oscillations. 2.3. Investigation of numerical solution for high peclet number We are interested in the behavior of the discrete solution, Eq. (8), for h fixed and m going to zero (i.e., Pe going to +1). First note that if Pe ¼ 2hm 1; then 0 < 2hm 1 and hence r is approximated by: r ffi 1 – 4m/h. Then, substitution in Eq. (8) and neglecting higher orders of 2hm, lead to
9321
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330 4
1
x 10
1.2
n=10, nu=0.0000001, Pe=500000
n=11, nu=0.0000001, Pe=454550
1
0
0.8
-1
u
u
0.6
-2
0.4
-3
0.2
-4
0
-5
0
0.2
0.4
0.6
0.8
1
-0.2
0
x
0.2
0.4
0.6
0.8
1
x
Fig. 2. Oscillation of the discrete solution for large values of Pe for even and odd numbers of elements (n).
ui ffi
1 ð1 4m=hÞi 1 ð1Þi ½1 þ ð4m=hÞi n ffi 1 ð1 4m=hÞ 1 ð1Þn ½1 þ ð4m=hÞn
The asymptotic analysis of ui as m tends to zero leads to Eq. (9) where two cases: N even and N odd are distinguished.
(
if n is even; if n is odd;
ui ffi ui ffi
i n
if i is even
ðh=2mÞ 1n ni if i is odd ð2m=hÞi
if i is even
ð9Þ
1 ð2m=hÞðn iÞ if i is odd
Direct substitution of m = 0 in Eq. (9) shows that ui oscillates between 0 and 1 when n is odd. For even values of n, it is clear that ui is oscillating between a bounded value ðni Þ; and a negative unbounded value. In Fig. 2, the discrete solutions are plotted for m = 107 for even (n = 10) and odd (n = 11) numbers of intervals. As it is clear, the behavior is in full agreement with the previous analysis (see Eq. (9)). When n is even, the left figure, the discrete solution oscillates between a finite value at even nodes ði ¼ 0; 2; 4; . . .Þ, and negative very large values at odd nodes ði ¼ 1; 3; 5; . . .Þ. On the other hand, the right figure shows that the discrete solution oscillates between 0 at even nodes and 1 at odd nodes. 3. The modified diffusion coefficient method (MDC) In this section we prove that the discrete solution (Eq. (8)) of the class problem Eq. (3) coincides with the exact one (given by Eq. (4)) for every Pe if Galerkin finite-element is applied to Eq. (3) with its original diffusion coefficient m replaced by a modified coefficient m0. 3.1. Estimation of the modified diffusion parameter Let m be replaced by m + d in Eq. (3) then, the discrete system becomes
(
ðm þ dÞ u0 ¼ 0;
ui1 2ui þuiþ1 h
þ
uiþ1 ui1 2
¼ 0;
i ¼ 1; . . . ; n 1
un ¼ 1
ð10Þ
Introducing a parameter R, defined as
R¼
h 2ðm þ dÞ
ð11Þ
we obtain the discrete algebraic system
8 > < u0 ¼ 0 ðR 1Þui1 þ 2ui þ ðR 1Þuiþ1 ¼ 0; > : un ¼ 1
i ¼ 1; . . . ; n 1
ð12Þ
9322
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Eq. (12) is identical to Eq. (6) (just Pe is replaced by R). Now, define
r¼
1þR 1R
ð13Þ
Analogous to Eq. (8), the exact solution of the discrete system (12) is
(
i
1r ui ¼ 1r n
i ¼ 0; . . . ; n
if R – 1
ui ¼ 0 i ¼ 0; . . . ; n 1; un ¼ 1 if R ¼ 1 Substituting = xi/h, where h = 1/n, then for R – 1, the discrete solution at xi is given by
ui ¼
1 r xi =h 1 qxi ¼ ; 1 r 1=h 1q
q ¼ r1=h
ð14Þ
Since the exact solution Eq. (3) can be written as
u¼
ex=m 1 1 lx ¼ ; e1=m 1 1l
l ¼ e1=m
ð15Þ
It is easy to notice that discrete solution (14) coincides with the exact one (15) if r1/h = e1/m or equivalently, ¼ eh=m . Then, substitution from Eq. (11) yields 22mmþ2dþh ¼ eh=m . Solving for d yields þ2dh simplified as
ð1þR Þ1=h ¼ e1=m whichleadsto 1þR 1R h=m 1R mhÞð2mþhÞ d ¼ e ð22ð1e , which can be h=m Þ
d ¼ m½1 þ Pe cot hðP e Þ
ð16Þ
Thus, we have proven the following theorem. Theorem 3.1. The Galerkin finite-element solution of Eq. (3) with the diffusion coefficient m replaced by m0 = m + d = mPe cot h(Pe) produces the exact solution of Eq. (3) at every node in the mesh. It is important to notice that this theorem is independent of the mesh size and more importantly is independent on the u
2u þu
value of m. In fact, Theorem 3.1 proves that the addition of the term d i1 h i iþ1 to the left-hand side of Eq. (3) compensates the truncation error introduced by the numerical scheme. The relation between m and m0 is shown in Fig. 3 for the value of h = 1/50. It is noticed that, as m ? 0, m0 ? h/2, however, for greater values of m, the curves of both m0 and m coincide. 3.2. Generalization of the diffusion parameter First, a slight generalization of the class problem Eq. (2) is now considered as:
u00 þ wu0 ¼ 0 in Xða; bÞ uðaÞ ¼ ua ;
ð17Þ
uðbÞ ¼ ub
h/2=0.01
0.12 0.1
Modified nu nu
0.08
Modified nu
0.06 0.04 0.02 0
0
0.02
0.04
0.06
0.08
0.1
nu Fig. 3. Relation between the diffusion coefficient m and the modified one m0.
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
9323
Defining m = e/w, the exact solution of Eq. (17) is given by
u ¼ ua þ ðub ua Þ
ð1 eðxaÞ=m Þ ð1 eðbaÞ=m Þ
ð18Þ
To obtain the exact solution of the discrete algebraic system produced by applying classical finite-element, Eq. (17) is first xa a ¼ uuu x ¼ ba normalized to return to the form of Eq. (3). Let u and then Eq. (17) reduces to ua
(
b
u 00 þ u 0 ¼ 0 in Xð0; 1Þ m ð1Þ ¼ 1 uð0Þ ¼ 0; u
ð19Þ
2
1 ¼ w ba 0 ¼ ddux ; u 00 ¼ ddx2u. The exact solution of the algebraic system resulting by Galerkin finite-element on a uniwhere m and u formly divided domain can be obtained as
ui ¼ ua þ ðub ua Þ
1 r i 1 ðr 1=h Þxi a ¼ u þ ðu u Þ a a b 1 rn 1 ðr1=h Þba
ð20Þ
e where ui is the discrete solution at xi ¼ a þ ih; where h ¼ ðb aÞ=n and r ¼ 1þP , P e ¼ 2hm – 1. Using procedures similar to that 1P e
adopted in the previous section, it could be proven that the discrete finite-element solution will coincides on the exact solu is replaced in Eq. (19) by m 0 ¼ m Pe coth ½P e ; where, P e ¼ 2hm ; Pe ¼ 2hm. Then, m 0 ¼ 2h coth ð2hmÞ and we have the tion Eq. (18) if m following theorem. Theorem 3.2. The Galerkin finite-element solution of Eq. (17), with the diffusion coefficient m = e/w replaced by m0 ¼ ðb aÞmPe coth ðPe Þ; produces the exact solution of Eq. (17) at every node in the mesh. Other generalizations of Theorem 3.1 are: 1. On non-uniform grids, m0 is computed locally within each element e using the element mesh size he. Using Theorem 3.2, one can prove that exact solution is produced for the considered problem class Eq. (3) on non-uniform meshes. 2. For higher dimensional problems, m0 is computed differently in each direction according to the coefficients component in that direction. Although there is no guarantee that this process produces the exact solution, it is expected that it gives good approximate solutions. 3. For variable coefficients problems, m0 is computed locally within each element by approximating the coefficients by their corresponding values at the center of that element. 4. Numerical results In this section, we consider the linear and stationary diffusion–convection problem:
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð21Þ
where e is the diffusive coefficient, w is the transport advective field, f is the volume source term, X R2 is a smooth convex domain and g is the boundary value prescribed for u on C. 4.1. Problem 1 [15] The domain X is the unit square 0 6 x 6 1; 0 6 y 6 1, with mixed boundary conditions. More precisely, u = 0 on x = 0, u = 1 on x = 1 and @u ¼ 0 on y = 0 and y = 1. The given transport field is w = (w, 0) and f = 0. The importance of this problem @n is that it is of one-dimensional nature and is reducible to
u00 þ wu0 ¼ 0 in Xð0; 1Þ uð0Þ ¼ 0;
uð1Þ ¼ 1
which is identical to Eq. (2) that was analyzed in the previous sections. 1 As a first experiment, the domain of the problem is discretized uniformly with h ¼ 16 and the coefficients are chosen such 1 h that m ¼ w ¼ 50 producing Peclet number P e ¼ 2m ¼ 1:5625. To examine the proposed modified diffusion coefficient (MDC), the problem is solved directly by classical Galerkin finite-element, using stream diffusion (SD) method as implemented in [17] and using the proposed MDC method. In Fig. 4, the solution of the three methods are compared with the exact solution (Eq. (3)). As expected, the MDC solution coincides on the exact one while the Galerkin solution oscillates near the boundary region and SD, although prevents oscillations, fails to approximate the solution.
9324
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Nu=1/50, h=1/16, Pe=1.5625
1 0.8
Galerkin Exact Solution SD Proposed MDC
Solution u
0.6 0.4 0.2 0 -0.2 -0.4 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x Fig. 4. Comparison of solutions of problem 1 on a uniform grid by different numerical techniques.
Table 1 1 Maximum discretization error of problem 1 on uniform grid (h ¼ 16 ). Peclet number
Galerkin
SD
Proposed MSD
1.5625 2 3 5
0.2634 0.3516 0.5025 0.6693
0.0439 0.0183 0.0025 4.5400e-005
4.0870e-015 3.8129e-015 3.9166e-015 3.1866e-015
1
MDC Nu=1/160 Exact Nu=1/160 MDC Nu=1/48 Exact Nu=1/48 MDC Nu=1/320 Exact Nu=1/320
0.8
u
0.6 0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
x Fig. 5. Full agreement of the computed solution by MDC and the exact solution for different values of m on a non-uniform grid.
The previous experiment is repeated, but with different values of Pe and the results of maximum discretization error computed as Maxi juexact ui j, where ui is the computed solutions at node i are summarized in Table 1. The results of MDC are in full agreement with the analysis of the previous section. It produces the exact solution (up to round-off error 1015) for different values of Pe. In the previous results, a uniform grid was used. In the next experiment, a non-uniform grid is used where the domain 0 < x < 1 is subdivided into 16 intervals using a Shishkin grid [18]. Shishkin grids are specially designed to fit boundary layers where a transition line determines the boundary between the fine discretization that resolves the layer, and the coarse discretization outside of the layer. Using MDC on a piecewise uniform Shishkin grid, the results for different values of m are plotted in Fig. 5 showing that the computed discrete solution coincides on the curve representing the exact solution. On a nonuniform mesh, element Peclet number is computed locally within each element i as P ie ¼ 2hmi , where hi is the element length. Table 2 summarizes the results. As expected by the analysis of Section 3.2, MDC method produces the exact solution for the 1-D class of the convection– diffusion equation (Eq. (2)) even on non-uniform grids.
9325
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330 Table 2 Maximum discretization error on a non-uniform grid for problem 1.
m
P ie
Maxi juexact ui j
1/48 1/160 1/320
2.8336 9.4455 18.8910
2.9976e-015 3.0531e-015 3.2058e-015
4.2. Problem 2 [7] We next consider a constant coefficient convection–diffusion equation
r ðruÞ þ w ru ¼ f
in X
u¼g
on C
ð22Þ
defined on the unit square 0 < x < 1; 0 < y < 1 with w = (1,0). The boundary condition is prescribed as
uð0; yÞ ¼ sin py;
uðx; 0Þ ¼ uðx; 1Þ ¼ 0;
uð1; yÞ ¼ 2 sin py
The exact solution is [19]
u¼
1 ex=2 sinh r
sin pyð2ex=2
sinh rx þ sinh rð1 xÞÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ p2 þ 1=42 . This problem represents a convection-dominated flow and was used as one of the test problems by Gupta et al. [19]. The coefficient of the convective term is constant. Fig. 6 shows the exact solution with = 0.05. For most part of the domain, the exact solution is smooth. But it has a steep boundary layer along the downstream edge at x = 1. Problem 2 is solved using classical Galerkin finite-element on a uniform grid (h = 1/16) for e = 0.001 and e = 0.05. The solutions are plotted in Fig. 7 showing nonphysical oscillations that are spread along the whole domain for the lower value of e = 0.001 and is localized near the boundary region as e is increased to 0.05. This behavior is similar to that shown in Fig. 1 for the 1-D case.
2 1 0 1 0.5 0
0.2
0
0.8
1
0.2
0.4
0.6
0.4
Fig. 6. Exact solution of problem 2 (e = 0.05).
2
2 1.5
1
1 0.5
0
0 1
-1 1
0.5
0.5
0
0
0.2
0.4
0.6
0.8
1
0
0
0.6
0.8
(b) Fig. 7. Galerkin finite-element solution of problem 2 with uniform grid size h = 1/16 for different values of e.
1
9326
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
2 1 0 1 0.5 0
0
0.6
0.4
0.2
0.8
1
Fig. 8. Finite-element solution using MDC of problem 2 ( = 0.001, h = 1/16).
x 10
Maximum Discretization error norm
-4
2 1 0 1
0.5
0
1
0.5
0
Fig. 9. Distribution of the discretization error norm of the finite-element solution using MDC of problem 2 ( = 0.001, h = 1/16).
Table 3 Convergence of MDC finite-element method for problem 2. h
1/16
1/32
1/64
1/128
Maximum discretization error
1.3084e-004
5.8178e-005
2.7241e-005
1.3175e-005
GMG Convergence using MDC, nu=1/200
0
10
h=1/64 h=1/32 h=1/16
-2
10
-4
Residual norm
10
-6
10
-8
10
-10
10
-12
10
-14
10
-16
10
0
2
4
6
8
10
12
14
16
18
20
V-cycles Fig. 10. Convergence of GMG for different h for problem 2 using MDC.
Next, problem 2 is solved on the same uniform grid using MDC method and the solution for = 0.001is plotted in Fig. 8. The absolute discretization error distribution is shown in Fig. 9. It is noticed that although extending MDC to this 2-D convection–diffusion problem does not produce the exact solution, it provides non-oscillatory accurate solution even on a coarse uniform grid. Table 3 reports the maximum discretization error norm of MDC finite-element method for different mesh sizes. As expected, multigrid methods do not converge if classical Galerkin finite-element is used to discretize the convection– diffusion problems since it fails to provide good approximations on coarser grids. However, introducing the MDC-Galerkin
9327
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
AMG Convergence using MDC, nu=1/200
2
10
0
h=1/16 h=1/32 h=1/64 h=1/128
10
-2
10
Residual norm
-4
10
-6
10
-8
10
-10
10
-12
10
-14
10
-16
10
2
4
6
8
10
12
14
16
18
20
V-Cycles Fig. 11. Convergence of AMG for different h for problem 2 using MDC.
discretization succeeded to retain the convergence of both geometric (GMG) and algebraic (AMG) multigrid methods as shown in Figs. 10 and 11, respectively. Convergence of AMG is slightly better than that of GMG but convergence rates are dependent on the mesh size h for both methods. 4.3. Problem 3 [7] In this test problem, we chose in Eq. (22) the following convection coefficients
w ¼ ðxðx 1Þð1 2yÞ; yðy 1Þð1 2xÞÞ It is obvious that this test problem has a stagnation point at (0.5, 0.5). A stagnation point (x0, y0) is the one inside the computational domain where both convection coefficients vanish. Convection–diffusion equations with stagnation points in their domains are usually used to model recirculation flow problems. For small values of e, this type of problem is very hard to solve, especially when standard multigrid method is used as the solution technique. The exact solution is chosen similar to the one used [7] as
u ¼ ð1 aðxÞÞð1 aðyÞÞ where,
aðxÞ ¼
eð1xÞ= 1 eð1yÞ= 1 ; aðyÞ ¼ 1= 1= e 1 e 1
Computed solution epslon=0.1
Computed solution epslon=0.01
0
0
-0.2
-0.2 -0.4
-0.4
-0.6
-0.6
-0.8
-0.8 -1 0
0 0.2
0.5 0.4
0.6
0.8
1
1
-1 0
0 0.5
0.5 1
Fig. 12. Computed solution on a stretched grid of problem 3 with b = 1/20, n = 32 and different values of
1
= 0.1, = 0.01.
9328
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
An exponential stretched grid in both x- and y-directions is used. The interval (0,1) is divided into n-sub-intervals by n + 1 iQ
n ; nodes such that xi ¼ 1e 1eQ
i ¼ 0; . . . ; n. As i changes from 0 to n, xi gives the coordinates of nodes along the x-axis from 0 to
firstinterval ln b 1. If the stretching parameter Q is chosen such that the ratio b ¼ lengthof , then Q ¼ nn1 . lengthof lastinterval
The computed solution of problem 3 for = 0.1, = 0.01 on a stretched grid is shown in Fig. 12 where the solution has steep boundary layers near x = 0 and y = 0. The width of these boundary layers becomes thinner as e decreases and sharp variation in the solution is observed. It must be mentioned that the exact solutions are graphically indistinguishable from
Table 4 Convergence of MDC finite-element method for problem 3. e
0.10
0.05
0.01
Maximum discretization error
3.654 104
5.431 104
3.821 103
10 10
Residual error norm
10 10 10 10 10 10
GMRES with different preconditioners
1
AMG GMG ILU non
0
-1
-2
-3
-4
-5
-6
0
5
10
15
20
25
30
iterations Fig. 13. Convergence behavior of GMRES for different preconditioners for problem 3.
AMG Convergene for problem 3 with different stretching ratios B 1
10
B=0.01 B=0.1 B=0.05
0
Residual Error Norm
10
-1
10
-2
10
-3
10
-4
10
-5
10
-6
10
2
4
6
8
10
12
14
16
18
20
AMG Cycles Fig. 14. Convergence behavior of AMG for different stretching factors b for problem 3.
9329
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
the shown computed solutions. It is clear that although the solution seems smooth for e = 0.1, it becomes more steep as e decreases and hence obtaining accurate solutions becomes more difficult. This is clear in Table 4 which reports the maximum discretization error for different values of e on a grid with mesh size h = 1/32. The resulting algebraic linear system produced by applying the finite-element method with the modified diffusion coefficient of problem 3 is solved by the generalized minimum residual GMRES method [20] with and without preconditioners. The results in Fig. 13 show the effect of preconditioning on the convergence of the GMRES. For this problem, the best performance is attained by the algebraic multigrid (AMG) preconditioner; only 7 GMRES iterations were sufficient to reduce the residual error norm below 106 of its initial value. The convergence of the GMRES with geometric multigrid (GMG) and incomplete LU (ILU) preconditioners are slightly slower but are significantly better than GMRES without any preconditioner. It is important to mention that classical multigrid does not converge when used to solve the resulting system. However, AMG works efficiently as observed from Fig. 14 where the convergence behavior is shown for different stretching ratios b. 4.4. Problem 4 [2] Although the previous two problems are singularly perturbed, they seem to have quasi-1D character. Therefore, in this Section we consider a problem that has different variability in each direction and hence represents a real 2D convection diffusion problem. This problem is similar to problem 3 but the boundary conditions and right hand side f(x, y) are computed such that the exact solution is u(x, y) = sin px + cos 3py + exy.
Computed solution h=1/16
Computed solution h=1/64 3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0 1
1 0.5
1
0.8
0.6
0.5
0.4
0.5
0.2
0 0
0
0
Fig. 15. Computed solution of problem 4 on a uniform grid with different mesh sizes h = 1/16, h = 1/64.
AMG convergence h=1/16
0
10
h=1/32 h=1/64 h=1/128
Maximum Error Norm
0 1
-5
10
-10
10
-15
10
2
4
6
8
10
12
14
AMG V-cycles Fig. 16. Convergence behavior of AMG for different mesh sizes.
9330
S.A. Mohamed et al. / Applied Mathematics and Computation 219 (2013) 9317–9330
Table 5 Convergence of MDC finite-element method for problem 4. h
1/16
1/32
1/64
1/128
Maximum discretization error Error reduction factor
5.11 104 –
1.27 104 4.02
3.16 105 4.02
7.91 106 3.99
The computed solutions using MCD on uniform grids with different mesh sizes h = 1/16, h = 1/64 are plotted in Fig. 15. The algebraic systems produced by MCD were solved efficiently by AMG whose convergence rates are presented in Fig. 16. The discretization error versus mesh size are given in Table 5 which show that a reduction of the discretization error by a factor of 4 is attained when the mesh size is reduced by a factor of 2. 5. Conclusions In this paper, we investigate the idea that the discretization error introduced by a numerical scheme to solve a boundary value problem can be compensated. A one-dimensional convection–diffusion equation with known exact solution is considered. The problem is discretized by a classical linear finite-element technique (equivalent to central finite-difference) and the analytic solution of the discrete system was computed and compared to the exact one to get the condition for being identical. We choose this condition in the form of a change in the normalized diffusion coefficient and get an analytic formula for the modified diffusion coefficient as a function of the original coefficients and the mesh size. Numerical results are presented to demonstrate that implementation of the proposed technique produces the exact solution for this class of 1-D singularly-perturbed problems even on coarse grids with uniform or non-uniform mesh sizes. As an extension to 2-D and variable coefficients convection–diffusion problems, the obtained formula is computed locally within each element and each space direction to approximate the solution. To demonstrate this approach and investigate the convergence behavior of the multigrid method, several numerical experiments are presented and accurate solutions were obtained even on coarser grids. As a result, multigrid-based solvers retain its efficient convergence rates. References [1] M.M. Gupta, R.P. Manohar, J.W. Stephenson, A single cell high order scheme for the convection–diffusion equation with variable coefficients, Int. J. Numer. Methods Fluids 4 (1984) 641–651. [2] S. Karaa, J. Zhang, Convergence and performance of iterative methods for solving variable coefficient convection–diffusion equation with a fourth-order compact difference scheme, Comput. Math. Appl. 44 (2002) 457–479. [3] G.Q. Chen, Z. Gao, Z.F. Yang, A perturbational, h4 exponential finite difference scheme for convection diffusion equation, J. Comput. Phys. 104 (1993) 129–139. [4] Z.F. Tian, S.Q. Dai, High-order compact exponential finite difference methods for convection diffusion type problems, J. Comput. Phys. 220 (2007) 952– 974. [5] M.M. Gupta, J. Zhang, High accuracy multigrid solution of the 3D convection–diffusion equation, Appl. Math. Comput. 113 (2000) 249–274. [6] Y. Wang, J. Zhang, Fast and robust sixth-order multigrid computation for the three-dimensional convection-diffusion equation, J. Comput. Appl. Math. 234 (2010) 3496–3506. [7] L. Ge, J. Zhang, High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids, J. Comput. Phys. 171 (2001) 560–578. [8] Y. Ge, F. Cao, Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems, J. Comput. Phys. 230 (2011) 4051–4070. [9] J. Zhang, H. Sun, J.J. Zhao, High order compact scheme with multigrid local mesh refinement procedure for convection diffusion problems, Comput. Methods Appl. Mech. Eng. 191 (2002) 4661–4674. [10] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, Springer-Verlag, Berlin, Heidelberg, 2008. [11] A. Russo, Streamline-upwind Petrov/Galerkin method (SUPG) vs residual-free bubbles (RFB), Comput. Methods Appl. Mech. Eng. 195 (2006) 1608– 1620. [12] E.G. Carmo, G.B. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the streamline and approximate upwind/Petrov–Galerkin method, Comput. Methods Appl. Mech. Eng. 192 (2003) 3379–3396. [13] N. Madden, M. Stynes, Linear enhancements of the streamline diffusion method for convection–diffusion problems, Comput. Math. Appl. 32 (10) (1996) 29–42. [14] N. Kopteva, How accurate is the streamline-diffusion FEM inside characteristic (boundary and interior) layers?, Comput Methods Appl. Mech. Eng. 193 (2004) 4875–4889. [15] M. Farhloul, A.S. Mounim, A mixed-hybrid finite element method for convection–diffusion problems, Appl. Math. Comput. 171 (2005) 1037–1047. [16] X.-G. Li, C.K. Chan, S. Wang, The finite element method with weighted basis functions for singularly perturbed convection–diffusion problems, J. Comput. Phys. 195 (2004) 773–789. [17] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics, Oxford, New York, 2005. [18] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [19] M.M. Gupta, R.P. Manohar, J.W. Stephenson, High-order difference schemes for two-dimensional elliptic equations, Numer. Methods Partial Differ. Equ. 1 (1985) 71. [20] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, 2003.
