- Email: [email protected]
Downwind numbering: robust multigrid for convection-diffusion problems
MATHEMATICS ELSEVIER
Applied Numerical Mathematics 23 (1997) 177-192
Downwind numbering: robust multigrid for convection-diffusion problems Jtirgen Bey a, Gabriel Wittum b,* Mathematisches lnstitut, Universitiit Tiibingen, Auf der Morgenstelle 10, 72076 T~bingen, Germany b lnstitutfiir Computeranwendungen, Universitiit Stuttgart, Pfaffenwaldring 27, 70550 Stuttgart, Germany a
Received 8 December 1995; revised 11 September 1996
Abstract
In the present paper, we introduce and investigate a robust smoothing strategy for convection-diffusion problems in two and three space dimensions without any assumption on the grid structure. The main tool to obtain such a robust smoother for these problems is an ordering strategy for the grid points called "downwind numbering" which follows the flow direction and, combined with a Gauss-Seidel type smoother, yields robust multigrid convergence for adaptively refined grids, provided the convection field is cycle-free. The algorithms are of optimum complexity and the corresponding smoothers are shown to be robust in numerical tests.
1. Introduction
In recent years, the research interest in robust multigrid methods has grown significantly. This is due to the fact that modelling challenging application problems like problems from fluid mechanics typically results in singularly perturbed partial differential equations. Since early multigrid days, the basic problems with singularly perturbed equations have been known and several strategies to work around them have been developed. One of the most powerful of these is the robust smoothing strategy, introduced by Wesseling and Kettler in 1980 (cf. [20,26]). The basic idea is to choose a smoother which is a fast solver for the singular perturbation and thus allows the coarse-grid correction cycle to handle the rest. This strategy has been used in many practical applications and theoretically analyzed (see [6-8,28,29] and the references there). An overview can be found in [6-8]. Unfortunately, up to now, robust multigrid methods have usually required structured grids, which is in direct contradiction to the requirement for adaptive refinement caused by local phenomena which are typical for singularly perturbed problems, such as boundary and shear layers, transition regions * Corresponding author. E-mail: [email protected]. 0168-9274/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S 0 1 6 8 - 9 2 7 4 ( 9 6 ) 0 0 0 6 7 - 0
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
178
between phases, etc. Thus, we need robust multigrid methods which can be applied to locally refined grids. The topic of this paper is to introduce and investigate a robust smoothing strategy for convectiondiffusion problems in two and three space dimensions without any assumption on the grid structure. The main tool for obtaining such a robust smoother for these problems is an ordering strategy for the grid points which follows the flow direction and, combined with a Gauss-Seidel-type smoother, yields robust multigrid convergence for adaptively refined grids, provided the convection field is cycle-free. The idea of combining a Gauss-Seidel method with an ordering process following the flow direction to obtain a very fast solver for the convective terms is well known of old (see, e.g., [19,21]). Combinations of this idea with multigrid methods were given for special problems and structured grids only (see, e.g., [20]), where numbering is very easy. In the present paper, we give an ordering strategy called "downwind numbering" which is of optimum complexity. This downwind numbering proves to be robust not only for problems which are globally convection-dominated but also for problems which are convection-dominated in part of the region only. There, multigrid correction is shown to be really necessary. We base our investigations on the scalar equation -eAu+c.
Vu=f
in 12
(1.1)
with divc -- 0 and Dirichlet boundary conditions. Discretization grids are constructed by combining triangles and quadrilaterals in 2D and by tetrahedra in 3D. To obtain a stable refinement, we apply the strategy introduced by the first of the authors in [9]. On these grids, we discretize (1.1) by means of a finite volume method and a corresponding first-order upwinding of the convective terms. These strategies are described in Section 2. After reviewing the basic multigrid schemes, we describe two downwind numbering algorithms and some of their basic properties in Sections 3 and 4. Section 5 contains the numerical tests validating the robustness of the downwind numbered smoother.
2. Grid refinement and discretization Multigrid methods for the solution of partial differential equations are based on a successively refined hierarchy 7~, 7], T2,... of triangulations ("grids") 7~ of the solution domain 12. In the following, we consider 12 c IR3 and 7~ consisting of a finite number of tetrahedra, which cover 12 and do not intersect pairwise in their interior. Starting with an initial triangulation 7~, we require the sequence 7~ to satisfy the following conditions: Nestedness: Each tetrahedron T E 7~ is contained in one of 7k-1 which is called its father element. The comers of T are either comers or edge midpoints of its father. - Compatibility: The intersection of two adjacent elements is either empty or a common comer or a common edge or common side. - Stability: For a given tetrahedron T E 7~ we define the measure of degeneracy ~(T) by the quotient of the diameter of the biggest inscribed ball and the length of the longest edge of T. We then require 3(T) to be bounded away from zero uniformly for all T E 7~ and k > 0. In two dimensions, several well-known algorithms exist generating such a sequence of nested grids satisfying the requirements above. One of the best known was used by Bank in the p l t m g software package [3]. The strategy consists of so-called "red" refinements, subdividing triangles into four -
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
179
congruent ones by simply connecting the midpoints of their edges and a following green closure by bisection. Some additional rules ensure stability. Another common two-dimensional refinement strategy is the bisection of elements, used, e.g., by B~insch [3] and Rivara [24]. The refinement of tetrahedral grids is more complicated. In contrast to the 2D case, there is obviously no way of dividing any given tetrahedron into eight congruent ones. Nevertheless, it is possible to extend the 2D refinement methods to three dimensions. For this purpose, algorithms have been developed by several authors during the last years. Three-dimensional bisection methods have been presented by B~insch [3] and Maubach [22] whereas a 3D analogon to Banks refinement strategy was introduced by Zhang in [32] and one of the authors of the present paper [9]. The latter one is applied in our code and explained in the following. 2.1. Regular refinement strategy
First, the edges of each face triangle are connected in the same way as in the two-dimensional regular refinement (Fig. 1). Then we chop off four subtetrahedra at the corners, which are congruent to the original one. In the interior of the remaining octahedron, there are three parallelograms as shown in Fig. 2. Cutting the octahedron along two of these parallelograms, we obtain four more subtetrahedra.
Fig. 1. Connectededges of each face.
Fig. 2. Interior parallelograms.
Fig. 3. Interior diagonals.
180
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Each choice of two cut parallelograms can be characterized by one of three possible diagonals shown in Fig. 3. Note that the eight subtetrahedra are of equal volume, but the interior ones are not congruent to the father. In contrast to the bisection method, we do not have any compatibility problems here when refining neighbouring elements, because the subdivision of all faces is fully symmetric. To preserve stability, we have to choose the cutting diagonal carefully. The wrong choice may lead to degenerated elements. In [9, Section 2.1], however, one of the authors proposed a simple algorithm preserving the stability of the hierarchy of triangulations and generating at most three congruence classes for each initial tetrahedron no matter how many successive refinements are performed. The algorithm works as follows: The tetrahedron T with comers T = [x0, aT1, 3C2, if:3] is divided into the subtetrahedra Ti, 1 <, i <~ 8, given by
T1 :T2 := 7"3 := T4 := T5 := T6 := T7 := T8 :=
[xo, (xo + Xl)/2, (xo + x2)/2, (X0 -k- X3)/2], [(xo + a:l)/2, Xl, (Xl + x2)/2, (Xl q'- X3)/2], [(xo + x2)/2, (Xl + x2)/2, X2, (X2 q'- X3)/2], [(xo -4- x3)/2, (Xl q- 2:3)/2, (X2 -+-X3)/2, X3], [(x0 -+-Xl)/2, (oc0+ x2)/2, (x0 + x3)/2, (Xl + x3)/2], [(zo + :q)/2, (zo + z2)/2, (acl + x2)/2, (Xl+X3)/2], [(x0 -q- x2)/2, (Xl + X2)/2, (Xl q- X3)/2, (x2 + x3)/2], [(x0 + x2)/2, (xo q- x3)/2, (Xl + x3)/2, (x2 + x3)/2]-
In this description, the choice of the cutting diagonal corresponds to the order of the comers. The stability proof is based on the dissection of the unit cube into six tetrahedra which can be transformed into each other under a permutation of their coordinates. It turns out that applying the refinement algorithm to each of them yields the same triangulation as if we first divide the cube into eight subcubes which again are subdivided into six tetrahedra as indicated above. This process can be repeated and the stability follows by induction and the usual affine transformation argument. The details of the proof can be found in the appendix of [9].
2.2. The green closure in 3D As indicated above, there are no compatibility problems in the case of global refinements, but in order to solve our equations adaptively, i.e., on adaptively refined grids, we need a three-dimensional "green" closure to obtain compatible triangulations avoiding the so-called slave nodes. These do not represent degrees of freedom and can locally change sparsity. The other possibility for circumventing the closure would be a non-conforming approximation. As pointed out in [9, Section 3.4], it is sufficient to consider four additional refinement rules shown in Fig. 4, which are called "irregular" refinement in analogy to the 2D case. Likewise according to [9, Sections 1.2 and 2.2], we add the following two rules ensuring stability and compatibility of our locally refined triangulations: If three or more edges not lying on one face are to be refined, then the element is refined regularly. Elements resulting from irregular refinements must not be refined further. -
-
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
181
@ Fig. 4. Green closure in 3D.
The second rule forces the substitution of irregular refinements by regular ones, if at least one of the corresponding elements is marked for further subdivision. In the last section, we show some figures of tetrahedral meshes with locally refined regions.
2.3. Finite volume discretization and upwind stabilization We use a finite volume method (also called box method) for the discretisation of our convection diffusion equation, which can be interpreted as a conservation law. See, e.g., the paper of Hackbusch [17] for a detailed description in the 2D case, which is extended here to three dimensions. Following [17] for every triangulation T = 7~, we define a dual box mesh B consisting of boxes b~ for every node xi of 7- in such a way that xi E bi and the boxes bi E B cover ~ and have pairwise empty inner intersection. The three-dimensional box mesh is constructed by connecting the center of mass of every element to the center of mass of all its faces and these again to the midpoints of their edges. Thus, we obtain four portions of the element, each of the same volume and one for each vertex. The box bi belonging to xi is then defined to be the union of all portions containing xi in a tetrahedron which has xi as corner. According to Hackbusch [17] we are able to construct the following second-order discretization under usual finite element assumptions. Using piecewise linear ansatz functions on 7- represented in terms of the corresponding nodal basis {~i}, integration of our convection-diffusion equation over the box bi and Gauss' theorem lead us to
~_uj f (-cV~j +c.~pj)d~= / fdx. J
~bi
bi
This is the discrete linear system with respect to 7-, and the discrete solution u is given by u = E j Uj~)j. Unfortunately, finite volume as well as finite element discretizations are known to be unstable in the case of dominating convection, which is observable in strong oscillations in the discrete solution. These are explained by the fact that the stiffness matrix is nearly skew symmetric and not an M-matrix. In special cases, the linear system can even be singular. To avoid those stability problems, we use a first-order upwind scheme (see, e.g., [4]), resulting in an O(h)-discretization. In one dimension, this upwind method can be interpreted as the substitution of central by upstream differences or as introducing an artificial diffusion. In the following we describe the finite volume method in detail.
182
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Fig. 5. Box portion of an elements comer. "V_
X0
X2
Xl Fig. 6. Segment vector and center of mass. To calculate the integrals along the boundary 0bi we proceed as follows: Let T/ be the set of tetrahedra T E 7- with comer xi. Then bi is a subset of the union of the elements T C T/and we can split the integrals into
f(--eV~j ~bi
+c~)j)dcr = ~ f (-eV~j +c~j)da. TCTi TN~bi
On account of the Dirichlet boundary conditions, it is sufficient to consider interior points xi. In this case T N ~bi consists of three plane segments that are located in the interior of T, as can be seen from Fig. 5. Thus each integral over a boundary portion T N Obi can be split again into a sum of three integrals running over the three plane segments. For simplicity assume that the comers of T are numbered as in Fig. 6. To each pair (i, k) satisfying 0 ~< i ~ k ~< 3 we associate a vector Sik with the following properties: (i) The length of Sik is given by the area of the plane segment between bi and bk. (ii) Sik/[Sik[ coincides with the outer normal of bi at the plane segment between bi and bk. Denoting by sik the center of mass of the segment between bi and bk, we obtain the quadrature formula
f TN~bi
+
da
k~i
+
(2.1)
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
183
The vectors Sik and the integration points 8ik can be computed as follows: Let 0 ~< g, g' <~ 3 be the remaining comer subscripts not equal to i, k. Then we have = 17
+ xk/+
+
and Sik = + ~
2
xe
x
2
xe,
,
where x denotes the usual vector product in R 3 and the direction of Sik is chosen in such a way that Sik " (xk - xi) > O. In addition ~,
j=iorj=k,
holds. After computing these values the quadrature formula (2.1) is easily applied to evaluate the integrals over the boundary of a box bi. The numbering process described below is based on the directed graph generated by upwinding. The application of the same process to higher order upwind schemes is straightforward.
3. Adaptive multigrid methods for the convection-diffusion problem One of the main difficulties of applying multigrid to convection-diffusion equations is the singularly perturbed character of this problem. Discretizations mix up convection and diffusion operators specifically with respect to the grid scale. This deteriorates the coarse grid approximation within the multigrid process. Further convection-diffusion equations show local phenomena as boundary or interior layers which have to be resolved by adaptive gridding. Thus, we need a robust, adaptive multigrid method to treat this type of equations. An overview of robust and adaptive multigrid is given in [6-8]. As pointed out there, one typical approach is to use a robust smoother and a standard coarse-grid correction cycle. In this case, the smoother has to eliminate the influence of the singular perturbation part such that the coarse-grid correction is able to do its work on the elliptic part. That means that the smoother should be exact or very fast for the limit case. This idea was introduced quite early by Wesseling [26,27] and investigated and extended by a couple of authors in the meantime [20,28,29]. For theoretical analysis of robust multigrid, we refer to [28,30]. In the case of convection-diffusion problems discretized by upwinding, it is quite easy to obtain an iterative scheme which solves the case of pure convection exactly. We can just choose a Gauss-Seidellike method and order the unknowns in flow direction. This trivial fact is well known in literature. If the diffusion part does not vanish, however, such a scheme is not sufficient as solver for the whole equation. Nevertheless, it is an ideal candidate as robust smoother within a multigrid cycle, since it satisfies the main requirement for a robust smoother explained above. Thus Gauss-Seidel and ILU schemes have been used in literature as robust smoother (see [20,27]) for two-dimensional convection-diffusion problems discretized on structured grids, where ordering is quite easy. For adaptively refined grids, we have to be aware of the complexity problem as well as of the robustness as described in [6,7]. This means we have to choose a multigrid method allowing for
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
184 Table 1
Smoothing pattern Additive
Multiplicative
new points only
hierarchical basis, Yserentant, 1984 [31]
HBMG, Bank, Dupont and Yserentant, 1987 [2]
whole refined region
BPX, Bramble, Pasciak and Xu, 1989 [10]
local muhigrid, Rivara, 1983 [24], Bramble, Pasciak, Wang and Xu, 1991 [11], Bastian, 1992 [5]
all points
mgp, Greenbaum, 1986 [15]
standard multigrid, Fedorenko, 1961 [14]
robust smoothing on the one hand and having optimal complexity in case of strongly local refinement. Several classes of adaptive multigrid methods exist, additive and multiplicative ones with global or local smoothing. As described in [6,7], these can be classified by their basic additive or multiplicative structure and the pattern which is used for smoothing as shown in Table 1. Table 1 refers to multigrid methods for unstructured grids. Refining structured grids by adding structured local patches, it is much easier to derive optimal multigrid methods. This has been done quite early by Brandt [12], Hackbusch [16] and McCormick [23]. As pointed out in [6], only local multigrid allows robust smoothing and has optimal complexity for adaptively refined grids. Since local multigrid is twice as efficient as BPX, we base our numerical tests on this algorithm. The kernel of our robust smoother is the ordering algorithm which is described in the next chapter.
4. Downwind numbering
4.1. Ordering algorithms Let the convection-diffusion equation
-eAu + c. V u = f,
(4.1)
with a vector c and c > 0 be given. This equation is discretized on a polyhedral grid by some grid method, leading to the stiffness matrix A = (aKL)K,LEA/', ./kf denoting the set of nodes where the unknowns are localized. By Aconv we denote the part of A corresponding to the discretization of the convective term c. V, by G = {(K, L): K, L E .h/', aconv,gL ~ 0} we denote the matrix graph. The main aim of the ordering process is to bring the stiffness matrix for the convective terms in triangular form. This can be performed easily, provided the flow field is vortex-free. When we discretize the convective part by the aforementioned upwind scheme, each edge in the matrix graph is assigned a direction. If the graph generated this way is cycle-free, it defines a partial ordering of the unknowns. This partial ordering can be used to number the unknowns in streamwise direction. In the following, we discuss two algorithms which order the cycle-free part of G.
Algorithm 4.1. Downwind_numbering_l (1) Assign the downwind direction from the discretization of the convective term to each link in the stiffness matrix graph. Indifferent links are marked by 0. (2) Put n = number of unknowns. (3) Compute I(K), the number of inflow edges of node K for all K E A/'.
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
(4) Find all minimal elements K E iV" with respect to the number of inflow edges them on to a fifo ~ . (5) Loop: While ~" not empty repeat: Get E from $'. (5a) Put Index(E) := 1. Put E on to a fifo 79. Let i := 1. (5b) Loop: While 79 not empty and i < n: Get K from 79. Loop: For all neighbours L of K do: If (L downwind from K) and (Index(L) ~< Index(K)), then i :-- Index(L) := Index(K) + 1; Put L onto 79. (6) Sort the vertex list such that Index(L) < Index(K) ~ L < K. Output: Ordered vertex list. Remark
185
I(K)
and put
4.2.
- If the edge graph is cycle-free, loop (4b) terminates in O(n)-steps with 79 = 0. Loop (5) has complexity O(q. n) where q is the number of minimal elements in the edge graph, which is small. - In step (6) nodes with the same index may be numbered in arbitrary order, caring that the number of a node with higher index is increased at the same time. This is possible with O(n) complexity. Using quicksort for this purpose would spoil the optimal complexity. - If loop (5b) terminates with 79 ~ 0 and i ~> n, then the edge graph contains a cycle. The complexity of Algorithm 4.1 may be improved by the following algorithm. For the description, we assume that the matrix graph G does not contain cycles.
Downwind_Numbering_2 Downwind_Numbering_2(matrixgraph G) Algorithm
4.3.
begin (1) Assign the downwind direction from the discretization of the convective term to each link in the stiffness matrix graph. Indifferent links are marked by 0. (2) Compute I(K), the number of inflow edges of node K for all nodes of the matrix graph ~. (3) m := 0; /* global counter */ (4) For each K E A f do: #(K) := - 1 ; /* #(K) contains the number of node K. #(K) = - 1 means node not yet numbered */ (5) For each K E A/" do: if (I(K) : 0 and #(K) = - 1 ) Number(K); end;
Number(nodeK); begin #(K) := m; m : = r n + 1;
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
186
for all neighbours L of K do: if (L is downwind from K ) then begin I(L):=I(L)-I; if ( I ( i ) = O) Number(L); end end; This algorithm is of optimal complexity, since each edge and each node is processed only once. The correctness follows directly from the depth search applied here. In the case of cycles, one may use the algorithm described by Tarjan [25] to determine the strongly connected components containing cycles. Then apply the above algorithm to the cycle-free graph of the components. Using a block Gauss-Seidel smoother with blocks corresponding to the components, we obtain a robust smoother. This is reasonable for small cycles. The treatment of large cycles is discussed in [18].
4.2. Artificial cycles Cycles in the matrix graph G will typically occur as soon as the flow field c contains vortices. Even in the case of vortex-free flow fields, however, it is not always guaranteed, that on the discrete level the matrix graph produced by the finite volume upwinding is free of cycles as stated in the following remark. R e m a r k 4.4. Let c be vortex-free. Then inner element cycles may arise in the following situation: Consider the triangle ABC, which is shown in Fig. 7 with neighbouring elements and a section of the corresponding dual box mesh. The flux is assumed to be constant and coming from the right. The boundary between adjacent boxes consists of two straight line segments. The total flux from one box into another is given by the sum of the fluxes over both segments. Of course the total flux from box B into box C is positive, since the same is true for the two boundary segments in between. Therefore, the directed edge from C to B is marked as a downwind edge. A
Fig. 7. Examplesfor artificialcycles.
187
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
In contrast to that, consider the flux over the two segments between boxes A and B. The two partial fluxes have opposite signs, but cross the segments with the same angle by construction. Thus, the flux over the larger segment dominates, and the result is a positive flux from A to B. The same argument holds for the flux from C to A, which proves the occurrence of an inner element cycle. Of course, this counter example subsists on the partial fluxes with opposite signs. This case cannot occur if the two segments between boxes are parallel. In fact, that means that there is only one plane segment, as in the example between B and C. Therefore, if we use the center verticals for construction of the dual box mesh, instead of the connection to the center of mass, we can avoid inner element cycles at least in the case of constant convection. Unfortunately, this strategy is applicable only for non-obtuse triangles, since otherwise the center verticals do not meet inside. Furthermore, this box method usually just provides O(h) convergence, instead of the O(h 2) convergence of the center of mass strategy, which satisfies the so-called equilibrium condition (see [17]). Flow-fields c containing vortices are treated in [18].
5. Numerical experiments We show results for the convection--diffusion equation -eAu+c.
Vu=f
in 52
(5.1)
with Dirichlet boundary conditions in two and three space dimensions on adaptively refined grids. 5.1. Two-dimensional test
Experiment 5.1. We solve Eq. (5.1) on the unit square/2 :--= (0, 1) × (0, 1) with c ( x , y ) = (1, 1) "r
(5.2)
and the boundary conditions u(x, y) =,(01
for for(XE(xE{0'I}'[0.5,1],YEy=(0'I])U(xC0). [0,1], y = 1) U (x C [0,0.5), y----0),
(5.3)
Adaptive refinement is performed by a standard residual-based error estimator (see [13]), which resolves the boundary and shear layers introduced by the boundary conditions (5.3) by local refinement. As expected, local multigrid with a Gauss-Seidel smoother and Vl = l, v2 -- 2, V-cycle shows robustness for large e as shown in Table 2. Throughout the experiments we used standard bilinear grid transfers. E x p e r i m e n t 5.2. For the second experiment, we consider a problem which is convection-dominated only in part of the region. There, we really need multigrid to treat the diffusion-dominated parts. Using the smoother as solver, we will not get reasonable convergence. To test this effect, we choose c as follows: c -- (1 - sin[a[ [2{x + ¼} - 1] + 2 cos[a] [y
- ¼])4 ( cos[a],
sin[a]) T,
where a is the angle of attack. The boundary conditions are: u = 0 on
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
188
Table 2 C o n v e r g e n c e rate averaged over 10 steps Level
¢ = 10 - 2
e = 10 - 4
e -~ 10 - 6
1
8.47 x 10 -3
1.05 x 10 - 3
1.00 x 10 -3
2
3.10 x 10 - 2
4.92 x 10 - 3
4.88 x 10 - 3
3
5.10 x 10 - 2
8.58 x 10 - 2
9.91 x 10 - 2
4
8.10 x 10 - 2
1.53 X 10 - 2
1.43 x 10 - 2
5
1.36 x 10 - I
2.01 x 10 - 2
1.77 x 10 - 2
6
1.82 x 10 - j
2.90 x 10 - 2
2.23 x 10 - 2
7
1.96 X 10 - 1
5.38 X 10 - 2
3.60 x 10 - 2
0,15 l 0,1
0 1E+0
o-.-......~~^
I
I
!
1E+I
1E+2
1E+3
! 1E+4
l
I
I
!
1E+5
1E+6
1E+7
1E+8
1/e Fig. 8. R o b u s t n e s s o f a (1,1,V)-multiplicative multigrid with I L U - s m o o t h e r and d o w n w i n d n u m b e r i n g . T h e m e t h o d used 8 locally refined grids to discretize p r o b l e m (5.2)-(5.3) with over 10.000 u n k n o w n s on level 8. T h e c o n v e r g e n c e rate ~ ( 1 0 ) is averaged over 10 steps and refers to finest grid.
{(x,y): x = 0 ,
0~y~
1}U{(x,y): 0
y=l}U{(x,y):
x=l,
0~y
u {(x,y): 0
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
189
downwind numbering but without coarse grid correction as a solver, we end up with a convergence rate of 0.949 as well. This confirms the outlined concept of robust multigrid. 5.2. Three-dimensional results
All three dimensional tests were performed on a unit cube. The coarsest triangulation used consists of 48 tetrahedra obtained by dividing the cube into eight subcubes, which are further divided into 6 tetrahedra canonically. Since we use Dirichlet boundary conditions, there is exactly one unknown on the coarsest level To. This grid is uniformly refined twice leading to ~ , 72, and T3. Further refinement steps are performed adaptively by means of a residual based error indicator. Hereby we adjusted the tolerance value of the error indicator in such a way that the number of unknowns on level 4 respectively 5 appeared to be nearly the same for all problems under consideration. Each problem has been discretized by the finite volume upwind scheme described in Section 2.3. After numbering the unknowns of every level in downwind direction using algorithm D o w n w i n d _ N u m b e r i n g _ 2 , we applied a multigrid (1,1,V)-cycle with a Gauss-Seidel smoother and standard trilinear prolongations and restrictions to solve the arising linear system. Note that in all experiments the complete matrix graph for the convection terms on the adaptively refined grids of levels 4 and 5 contained artificial cycles as explained above. In these cases downwind numbering was applied to a reduced convection graph, containing only connections for which the corresponding matrix entry exceeded a given threshold. The latter one was choosen minimal, just big enough that the reduced graph was cycle-free. The starting vectors for all iterations have been obtained by nested iteration, i.e., by interpolation of the previous solution to the next finer grid. To illustrate the convergence behaviour, we define the average convergence rate n by
Ilrko 112 l/ko --
\ttr0LI2]
'
where IlrkH2 is the euclidean norm of the residual after k iteration steps, and k0 denotes the number of iteration steps necessary to reduce the norm of the initial residual r0 by a factor of at least 10 -14. Experiment 5.3. We solve Eq. (5.1) in f2 := (0, 1) 3 with c(x,y,z) :=c5.3:=(1-z
1-z
1) w in J),
(5.4)
and Dirichlet boundary conditions 1 u(z, y, z) =
0 <~ x , y ~< 0.5, z = 0 , '
0,
(5.5)
otherwise.
This is a strongly convection-dominated problem for small ~. The convergence results versus e on different levels are shown in Table 3. Apparently the average convergence rate decreases for small values of ~. For fixed ~ however, n increases in the case of adaptive refinements (levels 4 and 5). This effect is due to better resolution which makes the diffusive effects dominating locally. The behaviour of ~; on level 5 seems to reflect some influence of the reduced graph numbering mentioned above. Fig. 9 shows grid and solution of level 5 for the case ~ = l0 -6. The plot shows a cut of the unit cube along its main diagonal.
190
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Table 3 Convergence rate for different levels versus e Level
# unknowns
• (~ = 10-2)
e~(e = 10-4)
~ (e = 10 - 6 )
2
343
0.03
8.0 X 10-6
3.0 X 10-8
3
3375
0.09
5.5 X 10-5
4.0 X 10-8
4
"~ 16200
0.17
3.0 X 10-3
2.6 X 10-3
5
" 30500
0.22
1.3 X 10-2
0.11
Fig. 9. Grid and solution on level 5 (e = 10-6). Table 4 Convergence rate for different levels and values of e Level
# unknowns
~¢ (£ - 10-2)
~(c
~---
10-4)
/'~ (E = 10 -6 )
2
343
o. 17
o. 12
0.07
3
3375
0.25
0.17
0.13
4
~ 16200
0.25
0.21
0.15
5
,-~ 30500
0.25
0.22
0.19
E x p e r i m e n t 5.4. Finally, we c o n s i d e r a m i x e d p r o b l e m with both c o n v e c t i o n - a n d d i f f u s i o n - d o m i n a t e d parts o f the d o m a i n in a n a l o g y to E x p e r i m e n t 5.2. To this end, w e multiply the c o n v e c t i o n coefficient c5.3 f r o m E x p e r i m e n t 5.3 b y
q=
0,
x + y-
2z ~< 0,
1,
x+y-2zl>
(x + y - 2z) ~,
otherwise.
I,
(5.6)
q v a n i s h e s a b o v e the m a i n d i a g o n a l o f the cube, c o n s t a n t and equals 1 b e l o w the line f r o m ( 0 . 5 , 0 . 5 , 0) to (1, 1,0.5), and is interpolated quadratically in between.
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
191
Fig. 10. Grid and solution of Experiment 5.4 on level 5 (e = 10-6). Table 4 shows the convergence rates for this problem and different values of e. The robustness of this approach is fully confirmed by the results. Again we use a (1,1,V)-multigrid cycle as solver. In Fig. 10, grid and solution on level 4 are shown for e = 10 -6.
References [ 1] R.E. Bank, T. Dupont and H. Yserentant, The hierarchical basis multigrid method, Preprint SC-87-2 (1987). [2] R.E. Bank, PLTMG: A Software Package for Solving Elliptic Partial Differential Equations. Users Guide 6.0 (SIAM, Philadelphia, PA, 1990). [3] E. B~insch, Local mesh refinement in 2 and 3 dimensions, Impact Comput. Sci. Engrg. 3 (1991) 181-191. [4] E Bastian, Locally refined solution of unsymmetric and nonlinear problems, in: W. Hackbusch and G. Wittum, eds., Incomplete Decompositions--Theory, Algorithms, and Applications, Notes Numer. Fluid. Mech. 41 (Vieweg, Braunschweig, 1993). [5] E Bastian, Parallele adaptive Mehrgitterverfahren, Dissertation, UniversiO't Heidelberg (1994). [6] E Bastian, W. Hackbusch and G. Wittum, Additive and multiplicative multigrid--a comparison, ICA-Bericht N8/95, ICA, Universit~it Stuttgart (1995). [7] E Bastian and G. Wittum, Adaptivity and robustness, in: Adaptive Methods, Proceedings 9th GAMM Seminar, Notes Numer. Fluid Mech. 46 (Vieweg, Braunschweig, 1994). [8] P. Bastian and G. Wittum, On robust and adaptive multigrid methods, in: E Hemker and E Wesseling, eds., Multigrid Methods, Proceedings 4th European Multigrid Conference, Amsterdam (Birkh~iuser, Basel, 1994). [9] J. Bey, Tetrabedral grid refinement, Computing 55 (4) (1995) 355-378. [10] J.H. Bramble, J.E. Pasciak, J. Wang and J. Xu, Convergence estimates for multigrid algorithms without regularity assumptions, Math. Comp. 57 (1991) 23-45. [11] J.H. Bramble, J.E. Pasciak and J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (1990) 1-22. [12] A. Brandt, Multi-level adaptive solution to boundary-value problems, Math. Comp. 31 (1977). [13] K. Eriksson and C. Johnson, Error estimates and automatic time step control for non-linear parabolic problems, SIAM J. Numer. Anal 24 (1987) 12-23. [ 14] R.P. Fedorenko, Ein Relaxationsverfahren zur L0sung elliptiscber Differentialgleichungen, Zh. Vychisl. Mat. i Mat. Fiz. 1 (5) (1961) 1092-1096 (in Russian).
192
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
[15] A. Greenbaum, A multigrid method for multiprocessors, Appl. Math. Comput. 19 (1986) 75-88. [16] W. Hackbusch, Local defect correction method and domain decomposition techniques, Comput. Suppl. 5 (1984) 89-113. [17] W. Hackbusch, On first and second order box schemes, Computing (1989). [18] W. Hackbusch and T. Probst, Preprint, Kiel (1996). [19] C. Johnson, Flow-directed Gauss-Seidel iterative methods for stationary convection-diffusion problems, Preprint 1992-29, Department of Mathematics, Chalmers University of Technology, Grteborg (1992). [20] R. Kettler, Analysis and comparison of relaxation schemes in robust multigrid and preconditioned conjugate gradient methods, in: W. Hackbusch and U. Trottenberg, eds., Multigrid Methods, Lecture Notes in Mathematics 960 (Springer, Heidelberg, 1982). [21] R. Komhuber and G. Wittum, Discretization and iterative solution of convection-diffusion equations, in: W. Hackbusch and G. Wittum, eds., Incomplete Decompositions--Algorithms, Theory and Applications, Notes Numer. Fluid. Mech. 41 (Vieweg, Braunschweig, 1993). [22] J. Maubach, Iterative methods for non-linear partial differential equations, Proefschrift, Wiskunde en Informatica, Katholieke Universiteit Nijmegen (1991). [23] S. McCormick, Multilevel Adaptive Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1990). [24] M.C. Rivara, Algorithms for refining triangular grids for adaptive and multigrid techniques, Internat. J. Numer. Methods Engrg. 20 (1984) 745-756. [25] R. Tarjan, Depth-first search and linear graph algorithms, SlAM J. Comput. 1 (1972) 146-160. [26] P. Wesseling, A robust and efficient multigrid method, in: W. Hackbusch and U. Trottenberg, eds., Multigrid Methods, Lecture Notes in Mathematics 960 (Springer, Berlin, 1982). [27] P. Wesseling, Theoretical and practical aspects of a multigrid method, SIAMJ. Sci. Statist. Comput. 3 (1982) 387--407. [28] G. Wittum, On the robustness of ILU-smoothing, SIAMJ. Sci. Statist. Comput. 10 (1989) 699-717. [29] G. Wittum, Linear iterations as smoothers in multigrid methods, Impact Comput. Sci. Engrg. 1 (1989) 180-215. [30] G. Wittum, Shifted iterations, ICA-Bericht N4/95, ICA, Universit~t Stuttgart (1995). [31] H. Yserentant, U'ber die Aufspaltung von Finite-Element-R~iumen in Teilr~iume verschiedener Verfeinerungsstufen, Habilitationsschrift, RWTH Aachen (1984). [32] S. Zhang, Multi-level iterative techniques, Ph.D. thesis, Research Report No. 88020, Department of Mathematics, Pennstate University (1988).
Applied Numerical Mathematics 23 (1997) 177-192
Downwind numbering: robust multigrid for convection-diffusion problems Jtirgen Bey a, Gabriel Wittum b,* Mathematisches lnstitut, Universitiit Tiibingen, Auf der Morgenstelle 10, 72076 T~bingen, Germany b lnstitutfiir Computeranwendungen, Universitiit Stuttgart, Pfaffenwaldring 27, 70550 Stuttgart, Germany a
Received 8 December 1995; revised 11 September 1996
Abstract
In the present paper, we introduce and investigate a robust smoothing strategy for convection-diffusion problems in two and three space dimensions without any assumption on the grid structure. The main tool to obtain such a robust smoother for these problems is an ordering strategy for the grid points called "downwind numbering" which follows the flow direction and, combined with a Gauss-Seidel type smoother, yields robust multigrid convergence for adaptively refined grids, provided the convection field is cycle-free. The algorithms are of optimum complexity and the corresponding smoothers are shown to be robust in numerical tests.
1. Introduction
In recent years, the research interest in robust multigrid methods has grown significantly. This is due to the fact that modelling challenging application problems like problems from fluid mechanics typically results in singularly perturbed partial differential equations. Since early multigrid days, the basic problems with singularly perturbed equations have been known and several strategies to work around them have been developed. One of the most powerful of these is the robust smoothing strategy, introduced by Wesseling and Kettler in 1980 (cf. [20,26]). The basic idea is to choose a smoother which is a fast solver for the singular perturbation and thus allows the coarse-grid correction cycle to handle the rest. This strategy has been used in many practical applications and theoretically analyzed (see [6-8,28,29] and the references there). An overview can be found in [6-8]. Unfortunately, up to now, robust multigrid methods have usually required structured grids, which is in direct contradiction to the requirement for adaptive refinement caused by local phenomena which are typical for singularly perturbed problems, such as boundary and shear layers, transition regions * Corresponding author. E-mail: [email protected]. 0168-9274/97/$17.00 Copyright © 1997 Elsevier Science B.V. All rights reserved PII S 0 1 6 8 - 9 2 7 4 ( 9 6 ) 0 0 0 6 7 - 0
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
178
between phases, etc. Thus, we need robust multigrid methods which can be applied to locally refined grids. The topic of this paper is to introduce and investigate a robust smoothing strategy for convectiondiffusion problems in two and three space dimensions without any assumption on the grid structure. The main tool for obtaining such a robust smoother for these problems is an ordering strategy for the grid points which follows the flow direction and, combined with a Gauss-Seidel-type smoother, yields robust multigrid convergence for adaptively refined grids, provided the convection field is cycle-free. The idea of combining a Gauss-Seidel method with an ordering process following the flow direction to obtain a very fast solver for the convective terms is well known of old (see, e.g., [19,21]). Combinations of this idea with multigrid methods were given for special problems and structured grids only (see, e.g., [20]), where numbering is very easy. In the present paper, we give an ordering strategy called "downwind numbering" which is of optimum complexity. This downwind numbering proves to be robust not only for problems which are globally convection-dominated but also for problems which are convection-dominated in part of the region only. There, multigrid correction is shown to be really necessary. We base our investigations on the scalar equation -eAu+c.
Vu=f
in 12
(1.1)
with divc -- 0 and Dirichlet boundary conditions. Discretization grids are constructed by combining triangles and quadrilaterals in 2D and by tetrahedra in 3D. To obtain a stable refinement, we apply the strategy introduced by the first of the authors in [9]. On these grids, we discretize (1.1) by means of a finite volume method and a corresponding first-order upwinding of the convective terms. These strategies are described in Section 2. After reviewing the basic multigrid schemes, we describe two downwind numbering algorithms and some of their basic properties in Sections 3 and 4. Section 5 contains the numerical tests validating the robustness of the downwind numbered smoother.
2. Grid refinement and discretization Multigrid methods for the solution of partial differential equations are based on a successively refined hierarchy 7~, 7], T2,... of triangulations ("grids") 7~ of the solution domain 12. In the following, we consider 12 c IR3 and 7~ consisting of a finite number of tetrahedra, which cover 12 and do not intersect pairwise in their interior. Starting with an initial triangulation 7~, we require the sequence 7~ to satisfy the following conditions: Nestedness: Each tetrahedron T E 7~ is contained in one of 7k-1 which is called its father element. The comers of T are either comers or edge midpoints of its father. - Compatibility: The intersection of two adjacent elements is either empty or a common comer or a common edge or common side. - Stability: For a given tetrahedron T E 7~ we define the measure of degeneracy ~(T) by the quotient of the diameter of the biggest inscribed ball and the length of the longest edge of T. We then require 3(T) to be bounded away from zero uniformly for all T E 7~ and k > 0. In two dimensions, several well-known algorithms exist generating such a sequence of nested grids satisfying the requirements above. One of the best known was used by Bank in the p l t m g software package [3]. The strategy consists of so-called "red" refinements, subdividing triangles into four -
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
179
congruent ones by simply connecting the midpoints of their edges and a following green closure by bisection. Some additional rules ensure stability. Another common two-dimensional refinement strategy is the bisection of elements, used, e.g., by B~insch [3] and Rivara [24]. The refinement of tetrahedral grids is more complicated. In contrast to the 2D case, there is obviously no way of dividing any given tetrahedron into eight congruent ones. Nevertheless, it is possible to extend the 2D refinement methods to three dimensions. For this purpose, algorithms have been developed by several authors during the last years. Three-dimensional bisection methods have been presented by B~insch [3] and Maubach [22] whereas a 3D analogon to Banks refinement strategy was introduced by Zhang in [32] and one of the authors of the present paper [9]. The latter one is applied in our code and explained in the following. 2.1. Regular refinement strategy
First, the edges of each face triangle are connected in the same way as in the two-dimensional regular refinement (Fig. 1). Then we chop off four subtetrahedra at the corners, which are congruent to the original one. In the interior of the remaining octahedron, there are three parallelograms as shown in Fig. 2. Cutting the octahedron along two of these parallelograms, we obtain four more subtetrahedra.
Fig. 1. Connectededges of each face.
Fig. 2. Interior parallelograms.
Fig. 3. Interior diagonals.
180
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Each choice of two cut parallelograms can be characterized by one of three possible diagonals shown in Fig. 3. Note that the eight subtetrahedra are of equal volume, but the interior ones are not congruent to the father. In contrast to the bisection method, we do not have any compatibility problems here when refining neighbouring elements, because the subdivision of all faces is fully symmetric. To preserve stability, we have to choose the cutting diagonal carefully. The wrong choice may lead to degenerated elements. In [9, Section 2.1], however, one of the authors proposed a simple algorithm preserving the stability of the hierarchy of triangulations and generating at most three congruence classes for each initial tetrahedron no matter how many successive refinements are performed. The algorithm works as follows: The tetrahedron T with comers T = [x0, aT1, 3C2, if:3] is divided into the subtetrahedra Ti, 1 <, i <~ 8, given by
T1 :T2 := 7"3 := T4 := T5 := T6 := T7 := T8 :=
[xo, (xo + Xl)/2, (xo + x2)/2, (X0 -k- X3)/2], [(xo + a:l)/2, Xl, (Xl + x2)/2, (Xl q'- X3)/2], [(xo + x2)/2, (Xl + x2)/2, X2, (X2 q'- X3)/2], [(xo -4- x3)/2, (Xl q- 2:3)/2, (X2 -+-X3)/2, X3], [(x0 -+-Xl)/2, (oc0+ x2)/2, (x0 + x3)/2, (Xl + x3)/2], [(zo + :q)/2, (zo + z2)/2, (acl + x2)/2, (Xl+X3)/2], [(x0 -q- x2)/2, (Xl + X2)/2, (Xl q- X3)/2, (x2 + x3)/2], [(x0 + x2)/2, (xo q- x3)/2, (Xl + x3)/2, (x2 + x3)/2]-
In this description, the choice of the cutting diagonal corresponds to the order of the comers. The stability proof is based on the dissection of the unit cube into six tetrahedra which can be transformed into each other under a permutation of their coordinates. It turns out that applying the refinement algorithm to each of them yields the same triangulation as if we first divide the cube into eight subcubes which again are subdivided into six tetrahedra as indicated above. This process can be repeated and the stability follows by induction and the usual affine transformation argument. The details of the proof can be found in the appendix of [9].
2.2. The green closure in 3D As indicated above, there are no compatibility problems in the case of global refinements, but in order to solve our equations adaptively, i.e., on adaptively refined grids, we need a three-dimensional "green" closure to obtain compatible triangulations avoiding the so-called slave nodes. These do not represent degrees of freedom and can locally change sparsity. The other possibility for circumventing the closure would be a non-conforming approximation. As pointed out in [9, Section 3.4], it is sufficient to consider four additional refinement rules shown in Fig. 4, which are called "irregular" refinement in analogy to the 2D case. Likewise according to [9, Sections 1.2 and 2.2], we add the following two rules ensuring stability and compatibility of our locally refined triangulations: If three or more edges not lying on one face are to be refined, then the element is refined regularly. Elements resulting from irregular refinements must not be refined further. -
-
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
181
@ Fig. 4. Green closure in 3D.
The second rule forces the substitution of irregular refinements by regular ones, if at least one of the corresponding elements is marked for further subdivision. In the last section, we show some figures of tetrahedral meshes with locally refined regions.
2.3. Finite volume discretization and upwind stabilization We use a finite volume method (also called box method) for the discretisation of our convection diffusion equation, which can be interpreted as a conservation law. See, e.g., the paper of Hackbusch [17] for a detailed description in the 2D case, which is extended here to three dimensions. Following [17] for every triangulation T = 7~, we define a dual box mesh B consisting of boxes b~ for every node xi of 7- in such a way that xi E bi and the boxes bi E B cover ~ and have pairwise empty inner intersection. The three-dimensional box mesh is constructed by connecting the center of mass of every element to the center of mass of all its faces and these again to the midpoints of their edges. Thus, we obtain four portions of the element, each of the same volume and one for each vertex. The box bi belonging to xi is then defined to be the union of all portions containing xi in a tetrahedron which has xi as corner. According to Hackbusch [17] we are able to construct the following second-order discretization under usual finite element assumptions. Using piecewise linear ansatz functions on 7- represented in terms of the corresponding nodal basis {~i}, integration of our convection-diffusion equation over the box bi and Gauss' theorem lead us to
~_uj f (-cV~j +c.~pj)d~= / fdx. J
~bi
bi
This is the discrete linear system with respect to 7-, and the discrete solution u is given by u = E j Uj~)j. Unfortunately, finite volume as well as finite element discretizations are known to be unstable in the case of dominating convection, which is observable in strong oscillations in the discrete solution. These are explained by the fact that the stiffness matrix is nearly skew symmetric and not an M-matrix. In special cases, the linear system can even be singular. To avoid those stability problems, we use a first-order upwind scheme (see, e.g., [4]), resulting in an O(h)-discretization. In one dimension, this upwind method can be interpreted as the substitution of central by upstream differences or as introducing an artificial diffusion. In the following we describe the finite volume method in detail.
182
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Fig. 5. Box portion of an elements comer. "V_
X0
X2
Xl Fig. 6. Segment vector and center of mass. To calculate the integrals along the boundary 0bi we proceed as follows: Let T/ be the set of tetrahedra T E 7- with comer xi. Then bi is a subset of the union of the elements T C T/and we can split the integrals into
f(--eV~j ~bi
+c~)j)dcr = ~ f (-eV~j +c~j)da. TCTi TN~bi
On account of the Dirichlet boundary conditions, it is sufficient to consider interior points xi. In this case T N ~bi consists of three plane segments that are located in the interior of T, as can be seen from Fig. 5. Thus each integral over a boundary portion T N Obi can be split again into a sum of three integrals running over the three plane segments. For simplicity assume that the comers of T are numbered as in Fig. 6. To each pair (i, k) satisfying 0 ~< i ~ k ~< 3 we associate a vector Sik with the following properties: (i) The length of Sik is given by the area of the plane segment between bi and bk. (ii) Sik/[Sik[ coincides with the outer normal of bi at the plane segment between bi and bk. Denoting by sik the center of mass of the segment between bi and bk, we obtain the quadrature formula
f TN~bi
+
da
k~i
+
(2.1)
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
183
The vectors Sik and the integration points 8ik can be computed as follows: Let 0 ~< g, g' <~ 3 be the remaining comer subscripts not equal to i, k. Then we have = 17
+ xk/+
+
and Sik = + ~
2
xe
x
2
xe,
,
where x denotes the usual vector product in R 3 and the direction of Sik is chosen in such a way that Sik " (xk - xi) > O. In addition ~,
j=iorj=k,
holds. After computing these values the quadrature formula (2.1) is easily applied to evaluate the integrals over the boundary of a box bi. The numbering process described below is based on the directed graph generated by upwinding. The application of the same process to higher order upwind schemes is straightforward.
3. Adaptive multigrid methods for the convection-diffusion problem One of the main difficulties of applying multigrid to convection-diffusion equations is the singularly perturbed character of this problem. Discretizations mix up convection and diffusion operators specifically with respect to the grid scale. This deteriorates the coarse grid approximation within the multigrid process. Further convection-diffusion equations show local phenomena as boundary or interior layers which have to be resolved by adaptive gridding. Thus, we need a robust, adaptive multigrid method to treat this type of equations. An overview of robust and adaptive multigrid is given in [6-8]. As pointed out there, one typical approach is to use a robust smoother and a standard coarse-grid correction cycle. In this case, the smoother has to eliminate the influence of the singular perturbation part such that the coarse-grid correction is able to do its work on the elliptic part. That means that the smoother should be exact or very fast for the limit case. This idea was introduced quite early by Wesseling [26,27] and investigated and extended by a couple of authors in the meantime [20,28,29]. For theoretical analysis of robust multigrid, we refer to [28,30]. In the case of convection-diffusion problems discretized by upwinding, it is quite easy to obtain an iterative scheme which solves the case of pure convection exactly. We can just choose a Gauss-Seidellike method and order the unknowns in flow direction. This trivial fact is well known in literature. If the diffusion part does not vanish, however, such a scheme is not sufficient as solver for the whole equation. Nevertheless, it is an ideal candidate as robust smoother within a multigrid cycle, since it satisfies the main requirement for a robust smoother explained above. Thus Gauss-Seidel and ILU schemes have been used in literature as robust smoother (see [20,27]) for two-dimensional convection-diffusion problems discretized on structured grids, where ordering is quite easy. For adaptively refined grids, we have to be aware of the complexity problem as well as of the robustness as described in [6,7]. This means we have to choose a multigrid method allowing for
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
184 Table 1
Smoothing pattern Additive
Multiplicative
new points only
hierarchical basis, Yserentant, 1984 [31]
HBMG, Bank, Dupont and Yserentant, 1987 [2]
whole refined region
BPX, Bramble, Pasciak and Xu, 1989 [10]
local muhigrid, Rivara, 1983 [24], Bramble, Pasciak, Wang and Xu, 1991 [11], Bastian, 1992 [5]
all points
mgp, Greenbaum, 1986 [15]
standard multigrid, Fedorenko, 1961 [14]
robust smoothing on the one hand and having optimal complexity in case of strongly local refinement. Several classes of adaptive multigrid methods exist, additive and multiplicative ones with global or local smoothing. As described in [6,7], these can be classified by their basic additive or multiplicative structure and the pattern which is used for smoothing as shown in Table 1. Table 1 refers to multigrid methods for unstructured grids. Refining structured grids by adding structured local patches, it is much easier to derive optimal multigrid methods. This has been done quite early by Brandt [12], Hackbusch [16] and McCormick [23]. As pointed out in [6], only local multigrid allows robust smoothing and has optimal complexity for adaptively refined grids. Since local multigrid is twice as efficient as BPX, we base our numerical tests on this algorithm. The kernel of our robust smoother is the ordering algorithm which is described in the next chapter.
4. Downwind numbering
4.1. Ordering algorithms Let the convection-diffusion equation
-eAu + c. V u = f,
(4.1)
with a vector c and c > 0 be given. This equation is discretized on a polyhedral grid by some grid method, leading to the stiffness matrix A = (aKL)K,LEA/', ./kf denoting the set of nodes where the unknowns are localized. By Aconv we denote the part of A corresponding to the discretization of the convective term c. V, by G = {(K, L): K, L E .h/', aconv,gL ~ 0} we denote the matrix graph. The main aim of the ordering process is to bring the stiffness matrix for the convective terms in triangular form. This can be performed easily, provided the flow field is vortex-free. When we discretize the convective part by the aforementioned upwind scheme, each edge in the matrix graph is assigned a direction. If the graph generated this way is cycle-free, it defines a partial ordering of the unknowns. This partial ordering can be used to number the unknowns in streamwise direction. In the following, we discuss two algorithms which order the cycle-free part of G.
Algorithm 4.1. Downwind_numbering_l (1) Assign the downwind direction from the discretization of the convective term to each link in the stiffness matrix graph. Indifferent links are marked by 0. (2) Put n = number of unknowns. (3) Compute I(K), the number of inflow edges of node K for all K E A/'.
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
(4) Find all minimal elements K E iV" with respect to the number of inflow edges them on to a fifo ~ . (5) Loop: While ~" not empty repeat: Get E from $'. (5a) Put Index(E) := 1. Put E on to a fifo 79. Let i := 1. (5b) Loop: While 79 not empty and i < n: Get K from 79. Loop: For all neighbours L of K do: If (L downwind from K) and (Index(L) ~< Index(K)), then i :-- Index(L) := Index(K) + 1; Put L onto 79. (6) Sort the vertex list such that Index(L) < Index(K) ~ L < K. Output: Ordered vertex list. Remark
185
I(K)
and put
4.2.
- If the edge graph is cycle-free, loop (4b) terminates in O(n)-steps with 79 = 0. Loop (5) has complexity O(q. n) where q is the number of minimal elements in the edge graph, which is small. - In step (6) nodes with the same index may be numbered in arbitrary order, caring that the number of a node with higher index is increased at the same time. This is possible with O(n) complexity. Using quicksort for this purpose would spoil the optimal complexity. - If loop (5b) terminates with 79 ~ 0 and i ~> n, then the edge graph contains a cycle. The complexity of Algorithm 4.1 may be improved by the following algorithm. For the description, we assume that the matrix graph G does not contain cycles.
Downwind_Numbering_2 Downwind_Numbering_2(matrixgraph G) Algorithm
4.3.
begin (1) Assign the downwind direction from the discretization of the convective term to each link in the stiffness matrix graph. Indifferent links are marked by 0. (2) Compute I(K), the number of inflow edges of node K for all nodes of the matrix graph ~. (3) m := 0; /* global counter */ (4) For each K E A f do: #(K) := - 1 ; /* #(K) contains the number of node K. #(K) = - 1 means node not yet numbered */ (5) For each K E A/" do: if (I(K) : 0 and #(K) = - 1 ) Number(K); end;
Number(nodeK); begin #(K) := m; m : = r n + 1;
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
186
for all neighbours L of K do: if (L is downwind from K ) then begin I(L):=I(L)-I; if ( I ( i ) = O) Number(L); end end; This algorithm is of optimal complexity, since each edge and each node is processed only once. The correctness follows directly from the depth search applied here. In the case of cycles, one may use the algorithm described by Tarjan [25] to determine the strongly connected components containing cycles. Then apply the above algorithm to the cycle-free graph of the components. Using a block Gauss-Seidel smoother with blocks corresponding to the components, we obtain a robust smoother. This is reasonable for small cycles. The treatment of large cycles is discussed in [18].
4.2. Artificial cycles Cycles in the matrix graph G will typically occur as soon as the flow field c contains vortices. Even in the case of vortex-free flow fields, however, it is not always guaranteed, that on the discrete level the matrix graph produced by the finite volume upwinding is free of cycles as stated in the following remark. R e m a r k 4.4. Let c be vortex-free. Then inner element cycles may arise in the following situation: Consider the triangle ABC, which is shown in Fig. 7 with neighbouring elements and a section of the corresponding dual box mesh. The flux is assumed to be constant and coming from the right. The boundary between adjacent boxes consists of two straight line segments. The total flux from one box into another is given by the sum of the fluxes over both segments. Of course the total flux from box B into box C is positive, since the same is true for the two boundary segments in between. Therefore, the directed edge from C to B is marked as a downwind edge. A
Fig. 7. Examplesfor artificialcycles.
187
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
In contrast to that, consider the flux over the two segments between boxes A and B. The two partial fluxes have opposite signs, but cross the segments with the same angle by construction. Thus, the flux over the larger segment dominates, and the result is a positive flux from A to B. The same argument holds for the flux from C to A, which proves the occurrence of an inner element cycle. Of course, this counter example subsists on the partial fluxes with opposite signs. This case cannot occur if the two segments between boxes are parallel. In fact, that means that there is only one plane segment, as in the example between B and C. Therefore, if we use the center verticals for construction of the dual box mesh, instead of the connection to the center of mass, we can avoid inner element cycles at least in the case of constant convection. Unfortunately, this strategy is applicable only for non-obtuse triangles, since otherwise the center verticals do not meet inside. Furthermore, this box method usually just provides O(h) convergence, instead of the O(h 2) convergence of the center of mass strategy, which satisfies the so-called equilibrium condition (see [17]). Flow-fields c containing vortices are treated in [18].
5. Numerical experiments We show results for the convection--diffusion equation -eAu+c.
Vu=f
in 52
(5.1)
with Dirichlet boundary conditions in two and three space dimensions on adaptively refined grids. 5.1. Two-dimensional test
Experiment 5.1. We solve Eq. (5.1) on the unit square/2 :--= (0, 1) × (0, 1) with c ( x , y ) = (1, 1) "r
(5.2)
and the boundary conditions u(x, y) =,(01
for for(XE(xE{0'I}'[0.5,1],YEy=(0'I])U(xC0). [0,1], y = 1) U (x C [0,0.5), y----0),
(5.3)
Adaptive refinement is performed by a standard residual-based error estimator (see [13]), which resolves the boundary and shear layers introduced by the boundary conditions (5.3) by local refinement. As expected, local multigrid with a Gauss-Seidel smoother and Vl = l, v2 -- 2, V-cycle shows robustness for large e as shown in Table 2. Throughout the experiments we used standard bilinear grid transfers. E x p e r i m e n t 5.2. For the second experiment, we consider a problem which is convection-dominated only in part of the region. There, we really need multigrid to treat the diffusion-dominated parts. Using the smoother as solver, we will not get reasonable convergence. To test this effect, we choose c as follows: c -- (1 - sin[a[ [2{x + ¼} - 1] + 2 cos[a] [y
- ¼])4 ( cos[a],
sin[a]) T,
where a is the angle of attack. The boundary conditions are: u = 0 on
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
188
Table 2 C o n v e r g e n c e rate averaged over 10 steps Level
¢ = 10 - 2
e = 10 - 4
e -~ 10 - 6
1
8.47 x 10 -3
1.05 x 10 - 3
1.00 x 10 -3
2
3.10 x 10 - 2
4.92 x 10 - 3
4.88 x 10 - 3
3
5.10 x 10 - 2
8.58 x 10 - 2
9.91 x 10 - 2
4
8.10 x 10 - 2
1.53 X 10 - 2
1.43 x 10 - 2
5
1.36 x 10 - I
2.01 x 10 - 2
1.77 x 10 - 2
6
1.82 x 10 - j
2.90 x 10 - 2
2.23 x 10 - 2
7
1.96 X 10 - 1
5.38 X 10 - 2
3.60 x 10 - 2
0,15 l 0,1
0 1E+0
o-.-......~~^
I
I
!
1E+I
1E+2
1E+3
! 1E+4
l
I
I
!
1E+5
1E+6
1E+7
1E+8
1/e Fig. 8. R o b u s t n e s s o f a (1,1,V)-multiplicative multigrid with I L U - s m o o t h e r and d o w n w i n d n u m b e r i n g . T h e m e t h o d used 8 locally refined grids to discretize p r o b l e m (5.2)-(5.3) with over 10.000 u n k n o w n s on level 8. T h e c o n v e r g e n c e rate ~ ( 1 0 ) is averaged over 10 steps and refers to finest grid.
{(x,y): x = 0 ,
0~y~
1}U{(x,y): 0
y=l}U{(x,y):
x=l,
0~y
u {(x,y): 0
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
189
downwind numbering but without coarse grid correction as a solver, we end up with a convergence rate of 0.949 as well. This confirms the outlined concept of robust multigrid. 5.2. Three-dimensional results
All three dimensional tests were performed on a unit cube. The coarsest triangulation used consists of 48 tetrahedra obtained by dividing the cube into eight subcubes, which are further divided into 6 tetrahedra canonically. Since we use Dirichlet boundary conditions, there is exactly one unknown on the coarsest level To. This grid is uniformly refined twice leading to ~ , 72, and T3. Further refinement steps are performed adaptively by means of a residual based error indicator. Hereby we adjusted the tolerance value of the error indicator in such a way that the number of unknowns on level 4 respectively 5 appeared to be nearly the same for all problems under consideration. Each problem has been discretized by the finite volume upwind scheme described in Section 2.3. After numbering the unknowns of every level in downwind direction using algorithm D o w n w i n d _ N u m b e r i n g _ 2 , we applied a multigrid (1,1,V)-cycle with a Gauss-Seidel smoother and standard trilinear prolongations and restrictions to solve the arising linear system. Note that in all experiments the complete matrix graph for the convection terms on the adaptively refined grids of levels 4 and 5 contained artificial cycles as explained above. In these cases downwind numbering was applied to a reduced convection graph, containing only connections for which the corresponding matrix entry exceeded a given threshold. The latter one was choosen minimal, just big enough that the reduced graph was cycle-free. The starting vectors for all iterations have been obtained by nested iteration, i.e., by interpolation of the previous solution to the next finer grid. To illustrate the convergence behaviour, we define the average convergence rate n by
Ilrko 112 l/ko --
\ttr0LI2]
'
where IlrkH2 is the euclidean norm of the residual after k iteration steps, and k0 denotes the number of iteration steps necessary to reduce the norm of the initial residual r0 by a factor of at least 10 -14. Experiment 5.3. We solve Eq. (5.1) in f2 := (0, 1) 3 with c(x,y,z) :=c5.3:=(1-z
1-z
1) w in J),
(5.4)
and Dirichlet boundary conditions 1 u(z, y, z) =
0 <~ x , y ~< 0.5, z = 0 , '
0,
(5.5)
otherwise.
This is a strongly convection-dominated problem for small ~. The convergence results versus e on different levels are shown in Table 3. Apparently the average convergence rate decreases for small values of ~. For fixed ~ however, n increases in the case of adaptive refinements (levels 4 and 5). This effect is due to better resolution which makes the diffusive effects dominating locally. The behaviour of ~; on level 5 seems to reflect some influence of the reduced graph numbering mentioned above. Fig. 9 shows grid and solution of level 5 for the case ~ = l0 -6. The plot shows a cut of the unit cube along its main diagonal.
190
J. Bey, G. Wittum / Applied Numerical Mathematics 23 (1997) 177-192
Table 3 Convergence rate for different levels versus e Level
# unknowns
• (~ = 10-2)
e~(e = 10-4)
~ (e = 10 - 6 )
2
343
0.03
8.0 X 10-6
3.0 X 10-8
3
3375
0.09
5.5 X 10-5
4.0 X 10-8
4
"~ 16200
0.17
3.0 X 10-3
2.6 X 10-3
5
" 30500
0.22
1.3 X 10-2
0.11
Fig. 9. Grid and solution on level 5 (e = 10-6). Table 4 Convergence rate for different levels and values of e Level
# unknowns
~¢ (£ - 10-2)
~(c
~---
10-4)
/'~ (E = 10 -6 )
2
343
o. 17
o. 12
0.07
3
3375
0.25
0.17
0.13
4
~ 16200
0.25
0.21
0.15
5
,-~ 30500
0.25
0.22
0.19
E x p e r i m e n t 5.4. Finally, we c o n s i d e r a m i x e d p r o b l e m with both c o n v e c t i o n - a n d d i f f u s i o n - d o m i n a t e d parts o f the d o m a i n in a n a l o g y to E x p e r i m e n t 5.2. To this end, w e multiply the c o n v e c t i o n coefficient c5.3 f r o m E x p e r i m e n t 5.3 b y
q=
0,
x + y-
2z ~< 0,
1,
x+y-2zl>
(x + y - 2z) ~,
otherwise.
I,
(5.6)
q v a n i s h e s a b o v e the m a i n d i a g o n a l o f the cube, c o n s t a n t and equals 1 b e l o w the line f r o m ( 0 . 5 , 0 . 5 , 0) to (1, 1,0.5), and is interpolated quadratically in between.
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
191
Fig. 10. Grid and solution of Experiment 5.4 on level 5 (e = 10-6). Table 4 shows the convergence rates for this problem and different values of e. The robustness of this approach is fully confirmed by the results. Again we use a (1,1,V)-multigrid cycle as solver. In Fig. 10, grid and solution on level 4 are shown for e = 10 -6.
References [ 1] R.E. Bank, T. Dupont and H. Yserentant, The hierarchical basis multigrid method, Preprint SC-87-2 (1987). [2] R.E. Bank, PLTMG: A Software Package for Solving Elliptic Partial Differential Equations. Users Guide 6.0 (SIAM, Philadelphia, PA, 1990). [3] E. B~insch, Local mesh refinement in 2 and 3 dimensions, Impact Comput. Sci. Engrg. 3 (1991) 181-191. [4] E Bastian, Locally refined solution of unsymmetric and nonlinear problems, in: W. Hackbusch and G. Wittum, eds., Incomplete Decompositions--Theory, Algorithms, and Applications, Notes Numer. Fluid. Mech. 41 (Vieweg, Braunschweig, 1993). [5] E Bastian, Parallele adaptive Mehrgitterverfahren, Dissertation, UniversiO't Heidelberg (1994). [6] E Bastian, W. Hackbusch and G. Wittum, Additive and multiplicative multigrid--a comparison, ICA-Bericht N8/95, ICA, Universit~it Stuttgart (1995). [7] E Bastian and G. Wittum, Adaptivity and robustness, in: Adaptive Methods, Proceedings 9th GAMM Seminar, Notes Numer. Fluid Mech. 46 (Vieweg, Braunschweig, 1994). [8] P. Bastian and G. Wittum, On robust and adaptive multigrid methods, in: E Hemker and E Wesseling, eds., Multigrid Methods, Proceedings 4th European Multigrid Conference, Amsterdam (Birkh~iuser, Basel, 1994). [9] J. Bey, Tetrabedral grid refinement, Computing 55 (4) (1995) 355-378. [10] J.H. Bramble, J.E. Pasciak, J. Wang and J. Xu, Convergence estimates for multigrid algorithms without regularity assumptions, Math. Comp. 57 (1991) 23-45. [11] J.H. Bramble, J.E. Pasciak and J. Xu, Parallel multilevel preconditioners, Math. Comp. 55 (1990) 1-22. [12] A. Brandt, Multi-level adaptive solution to boundary-value problems, Math. Comp. 31 (1977). [13] K. Eriksson and C. Johnson, Error estimates and automatic time step control for non-linear parabolic problems, SIAM J. Numer. Anal 24 (1987) 12-23. [ 14] R.P. Fedorenko, Ein Relaxationsverfahren zur L0sung elliptiscber Differentialgleichungen, Zh. Vychisl. Mat. i Mat. Fiz. 1 (5) (1961) 1092-1096 (in Russian).
192
J. Bey, G. Wittum/ Applied Numerical Mathematics 23 (1997) 177-192
[15] A. Greenbaum, A multigrid method for multiprocessors, Appl. Math. Comput. 19 (1986) 75-88. [16] W. Hackbusch, Local defect correction method and domain decomposition techniques, Comput. Suppl. 5 (1984) 89-113. [17] W. Hackbusch, On first and second order box schemes, Computing (1989). [18] W. Hackbusch and T. Probst, Preprint, Kiel (1996). [19] C. Johnson, Flow-directed Gauss-Seidel iterative methods for stationary convection-diffusion problems, Preprint 1992-29, Department of Mathematics, Chalmers University of Technology, Grteborg (1992). [20] R. Kettler, Analysis and comparison of relaxation schemes in robust multigrid and preconditioned conjugate gradient methods, in: W. Hackbusch and U. Trottenberg, eds., Multigrid Methods, Lecture Notes in Mathematics 960 (Springer, Heidelberg, 1982). [21] R. Komhuber and G. Wittum, Discretization and iterative solution of convection-diffusion equations, in: W. Hackbusch and G. Wittum, eds., Incomplete Decompositions--Algorithms, Theory and Applications, Notes Numer. Fluid. Mech. 41 (Vieweg, Braunschweig, 1993). [22] J. Maubach, Iterative methods for non-linear partial differential equations, Proefschrift, Wiskunde en Informatica, Katholieke Universiteit Nijmegen (1991). [23] S. McCormick, Multilevel Adaptive Methods for Partial Differential Equations (SIAM, Philadelphia, PA, 1990). [24] M.C. Rivara, Algorithms for refining triangular grids for adaptive and multigrid techniques, Internat. J. Numer. Methods Engrg. 20 (1984) 745-756. [25] R. Tarjan, Depth-first search and linear graph algorithms, SlAM J. Comput. 1 (1972) 146-160. [26] P. Wesseling, A robust and efficient multigrid method, in: W. Hackbusch and U. Trottenberg, eds., Multigrid Methods, Lecture Notes in Mathematics 960 (Springer, Berlin, 1982). [27] P. Wesseling, Theoretical and practical aspects of a multigrid method, SIAMJ. Sci. Statist. Comput. 3 (1982) 387--407. [28] G. Wittum, On the robustness of ILU-smoothing, SIAMJ. Sci. Statist. Comput. 10 (1989) 699-717. [29] G. Wittum, Linear iterations as smoothers in multigrid methods, Impact Comput. Sci. Engrg. 1 (1989) 180-215. [30] G. Wittum, Shifted iterations, ICA-Bericht N4/95, ICA, Universit~t Stuttgart (1995). [31] H. Yserentant, U'ber die Aufspaltung von Finite-Element-R~iumen in Teilr~iume verschiedener Verfeinerungsstufen, Habilitationsschrift, RWTH Aachen (1984). [32] S. Zhang, Multi-level iterative techniques, Ph.D. thesis, Research Report No. 88020, Department of Mathematics, Pennstate University (1988).