Analysis of linear multigrid methods for elliptic differential equations with discontinuous and anisotropic coefficients

METHODS FOR ELLIPTIC DIFFERENTIAL

EQUATIONS WITH DISCONTINUOUS AND

ANISOTROPIC COEFFICIENTS

MOHAMMED KHALIL

TR diss

1744

U

^ >

ANALYSIS OF LINEAR MULTIGRID METHODS

FOR ELLIPTIC DIFFERENTIAL EQUATIONS

WITH DISCONTINUOUS AND ANISOTROPIC COEFFICIENTS

PROEFSCHRIFT

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus,

Prof. Drs. P.A. Schenck,

in het openbaar te verdedigen

ten overstaan van een commissie

door het College van Dekanen

daartoe aangewezen,

op 22 juni 1989 te 16.00 uur

door

MOHAMMED KHALIL

Docteur de Troisième Cycle

geboren te Bensümaae^Marokko)

TR diss

1744

Summary

Chapter I

GENERAL INTRODUCTION

Chapter II

CELL-CENTERED AND VERTEX-CENTERED DISCRETIZATION OF

TWO-DIMENSIONAL INTERFACE PROBLEMS

Chapter III

VERTEX-CENTERED AND CELL-CENTERED TWO-DIMENSIONAL

MULTIGRID METHODS FOR INTERFACE PROBLEMS

Chapter IV

A CELL-CENTERED MULTIGRID METHOD FOR THREE-DIMENSIONAL

ANISOTROPIC DIFFUSION AND INTERFACE PROBLEMS

Chapter V

LOCAL MODE SMOOTHING ANALYSIS FOR VARIOUS INCOMPLETE

FACTORIZATION ITERATIVE METHODS IN TWO DIMENSIONS

Chapter VI

LOCAL MODE SMOOTHING ANALYSIS FOR VARIOUS INCOMPLETE

FACTORIZATION ITERATWE METHODS IN THREE DIMENSIONS

Chapter VII

Curriculum Vitae

In this thesis linear multigrid methods are analysed and tested for application to the

solution of elliptic partial differential equations with discontinuous and anisotropic

coefficients, discretized with a finite—volume method.

Two types of multigrid methods are compared: the vertex—centered and the

cell—centered multigrid method.

In the presence of discontinuous coefficients the vertex-centered method requires rather

complicated (matrix—dependent) operators between the grids. To define these operators

and the coarse grid operators both in two and three dimensions, use is made of a stencil

notation, which results in a clear and simple description of these complicated operators,

and straightforward algorithms to compute them.

A number of two—dimensional test problems are solved to compare the vertex-centered

and the (not yet much used) cell-centered method. The latter method is found to be

equally efficiënt but simpler than the former. The greater simplicity makes extension to

three dimensions easier.

For the cell-centered multigrid method transfer operators are developed in three

dimensions. They are chosen such that the coarse grid matrices resulting from Galerkin

approximation have optimal sparsity.

In order to cope with anisotropic coefficients, the smoother must be carfully chosen.

Smoothers based on incomplete factorizations are studied and their smoothing behaviour

is analysed with local mode (Fourier) analysis. The presence of boundaries is accounted

for in an approximate way, in order to obtain realistic predictions of the behaviour of the

methods in practice.

A new robust and simple smoothing method is found in two dimensions. In three

dimensions, only complex and costly robust smoothing methods are found. Satisfactory

practical results are obtained with a not completely robust but simple smoothing

method, accelerating the multigrid method with a Lanczos—type method (namely

Orthomin).

CHAPTERI

GENERAL INTRODUCTION

This study started in 1984 at Société Nationale Elf Aquitaine (Production) at

Pau (France). It was thought that one can develop a fast solver based on a multigrid

method for three-dimensional reservoir simulation linear systems in a short time. Indeed

it did not take long to develop an efficiënt two—dimensional code with vertex—centered

coarsening, matrix—dependent transfer operators and incomplete line factorization

smoothing. The code could not be used immediately for testing on two—dimensional

reservoir simulation systems. The reason was that the grid dimensions were not of type

(2

P

x 2

q

), p, q integers. To eliminate this restriction, we spent some time to find an

optima! "padding" (introduction of additional artificial grid points).

The extension to the three—dimensional case was not carried out easily. The

obstacles were:

1. Implementation of matrix-dependent transfer operators.

2. The cost of computation of their 26 entries per coarse node.

3. The high computational cost of Galerkin coarse grid approximation.

4. The choice of smoothing and its implementation for matrices with twenty—seven

diagonals.

As we used the CRAY—XMP computer, the cost associated with points 2 and 3 above

was reduced by vectorization. For smoothing we used simple point Gauss—Seidel as

finest grid smoother and ILU(O) [5] as coarse grids smoother. In real tests the code

worked satisfactory for interface problems but not for problems with anisotropy, because

the smoothing method was not efficiënt. The code was more time consuming than

conjugate gradiënt type methods preconditioned by ILU(O), ILU(l) or ILU(2) [5]. The

time gained by acceleration of multigrid by Orthomin was quite insufficiënt to

compensate the time consumed by preliminary work. Another weakness was the large

storage required for the coarse grid matrices and transfer operators. Nevertheless the

code was useful for systems with multiple right—hand—sides. This work is described in

part in [12].

The challenge remained to reduce the computational cost and the storage

requirement of the method and to find a simple robust smoother for problems with

discontinuous and anisotropic coefficients.

Because of its greater simplicity, the recently developed cell—centered multigrid

method [14] promises greater efficiency. In the present work we compare cell—centered

and vertex—centered multigrid method in the two-dimensional case. A smoothing

analysis of two—dimensional anisotropic diffusion problems is given for various

incomplete factorizations. A robust and simple smoother is found for such problems. An

extension of the cell—centered multigrid method to the three—dimensional case is

presented. Local mode smoothing analysis of various three—dimensional incomplete

factorizations is performed for anisotropic diffusion problems. A robust and simple

smoother IPFM(a>) is found for the case of dominant direction (large diffusion coëfficiënt

in one dominant direction). Unfortunately robustness is not obtained with IPFM(w) in

the case of dominating plane (large diffusion coefficients in two coordinate directions). In

order to annihilate the short wavelengths for which IPFM(w) does not work, acceleration

of multigrid by methods of Lanczos type may be effective. The local mode analysis shows

that alternating incomplete plane factorization is a robust smoother. Because of its lower

cost and greater simplicity, we prefer IPFM(w).

In the analysis and comparisons we make use of special cases of the following linear

second order elliptic partial differential equation.

_ J _ ^ _ [ D

a

| ^ ] + a ( x )

¥

< x ) = /(x) (1)

z e n c R

d

, de {2,3}

with D

a

> 0, <r > 0. D

a

, o and ƒ may have strong discontinuities across interna!

interfaces.

Equation (1) occurs in reservoir simulation models [2]. For general presentations of

multigrid methods we refer to [6, 8, 9, 11]. The cell-centered multigrid method is

presented in [14]. Vertex-centered multigrid methods for equations of type (1) are

presented in [1, 3, 4, 7, 12, 13].

REFERENCES

[1] Alcouffe R.E., Brandt A., Dendy J.E., Jr., Painter, J.W.: The multigrid method

for diffusion equations with strongly discontinuous coefficients. SIAM J. Sci. Stat.

Comp. 2 (1981), 430-454.

[2] Aziz, K. and Settari, A.: Petroleum reservoir simulation. Applied Science

Publishers, London, 1979.

[3] Behie, A. and Forsyth, P.A.: Multigrid solution of the pressure equation in

reservoir simulation. In: Proceedings of the Sixth SPE Symposium on Reservoir

Simulation, 29-42. Society of Petroleum Engineers of AIME, 1982. Paper SPE

10492.

[4] Behie, A. and Forsyth, P.A.: Multigrid solution of three—dimensional problems

with discontinuous coefficients. Appl. Math. Comp. 13 (1983), 229-240.

[5] Behie, A. and Vinsome, P.K.W.: Block Iterative Methods for Fully Implicit

Reservoir Simulation. Soc. Pet. Eng. J. (Oct. 1982), 658-668.

[6] Brandt, A.: Multi—level adaptive solution to boundary value problems. Math. of

Comput. 31 (1977), 333-390.

[7] Dendy, J.E., Jr.: Two multigrid methods for three-dimensional problems with

discontinuous coefficients. SIAM J. Sci. Stat. Comp. 8 (1987), 673-685.

[8] Hackbusch, W.: Multigrid method and applications. Springer-Verlag, Berlin

1985.

[9] Hackbusch, W. and Trottenberg, U. (eds.): Multigrid methods. Proceedings,

Köln-Porz 1981. Lecture Notes in Mathematics 960, Springer-Verlag, Berlin

1982.

[10] Kettler, R.: Linear multigrid methods for numerical reservoir simulation.

Ph.D.Thesis, University of Technology, Dept. of Technical Math. and

Informaties, Delft, December 1987.

[11] McCormick, S.F. (ed.): Multigrid methods, theory, applications and

supercomputing. Lecture Notes in Pure and Applied Mathematics, Marcel Dekker

Ine. New York, 1988.

[12] Puiseux, P.: Methode Multigrille pour la résolution de modèles mathématiques de

1'industrie pétrolière. These de Doctorat de 1'Université de Pau, Juillet 1988.

[13] Scott, T.: Multigrid method for oil reservoir simulation in two and three

dimensions. J. Comp. Phys. 59 (1985), 290-307.

[14] Wesseling, P.: Cell-centered multigrid for interface problems. J. Comp. Phys. 79

(1988), 85-91.

CHAPTER II*

CELL-CENTERED AND VERTEX-CENTERED DISCRETIZATION OF TWO-DIMENSIONAL

INTERFACE PROBLEMS

Page

l.INTRODUCTION 2.2 2. COMPUTATIONAL GRID 2.3

3. CELL-CENTERED DISCRETIZATION, CASE 1 2.4

3.1. Interior finite volumes 2.5 3.2. Boundary finite volumes 2.7 4. CELL-CENTERED DISCRETIZATION, CASE 2 2.9

4.1. Interior finite volumes 2.9 4.2. Boundary finite volumes 2.10 5. VERTEX-CENTERED DISCRETIZATION, CASE 3 2.11

5.1. Interior finite volumes 2.12 5.2. Boundary finite volumes 2.12 5.3. Corner finite volumes 2.13 6. VERTEX-CENTERED DISCRETIZATION, CASE 4 2.15

6.1. Interior finite volumes 2.15 6.2. Boundary finite volumes 2.15 6.3. Corner finite volumes 2.16

REFERENCES 2.17

1. INTRODUCTION

We consider the following elliptic boundary value problem in two dimensions:

a=l a a a=l

D (> 0), o (> 0) and f are given functions allowed to be discontinuous across internal interfaces,

b - f ^ + <tf = g, b > 0 , c > 0 , b + c # 0 , on dfl. (l.l.b) ön

Here b and c are functions. At each face of the d-dimensional domain f) we have one of the following cases:

Dirichlet: b = 0, c = 1; Neumann: b = 1, c = 0; Robbins: b = 1, c # 0.

Because of its importance in various fields of physics and engineering, fast numerical methods for equation (1.1) are of interest. Such methods have been proposed in [1,2,3,4,5]. The suitability of these methods has to be investigated by application to test problems. Numerical experiments have been described in the publications just mentioned.

In order to test a method developed in [7] and to compare it with some of the previous methods we have developed a computer program to generate discretizations of (l.l.a,b) using various methods. Cell-centered and vertex-centered finite volume discretization of the two-dimensional problem above is described. The material in this chapter is elemen-tary, but is presented in detail to facilitate reproduction of our numerical tests.

2. COMPUTATIONAL GRID

Let d = 2. The rectangle fl is divided uniformly into cells. The length of the cells in x -direction is h = L /N , with N

a a a a a Figure 2.1 gives a two- dimensional example.

the x -direction is h = L /N , with N the number of cells in the x -direction

a a a Q a a

Figure 2.1. Subdivision into cells.

The numerical unknowns are assigned to nodes. These nodes are positioned either at the centers of the cells (cell-centered grid) or at the vertices of the cells (vertex-centered grid). See Figure 2.2.

< 1 * 1 >

: 1 ~c 1 :

< l v I :

Cell-centered Vertex-centered

Figure 2.2. Cell-centered and vertex-centered grid.

For discretization we use the finite volume method. With cell-centered discretization the cells are the finite volumes. With vertex-centered discretization the finite volumes are centered on the vertices. Their boundaries are indicated by dashed Unes in Figure 2.2; at the boundaries their size differs from the size in the interior.

In the cell-centered case, D may be discontinuous across lines connecting cell centers, i.e. interfaces consist of line segments connecting cell centers (case 1); or interfaces coincide with cell boundaries (case-2)..In the vertex-centered case, interfaces consist of line segments connecting vertices (case 3), or along line segments midway between vertices ( — in Figure 2.2) (case 4).

The functions D , o and f are assumed to be given in certain points in the grid, de-pending on which case we have, according to Table 2.1. Vertices and cell centers are indicated by x and •, respectively, in Figure 2.2.

Cell-centered Vertex-centered Case 1 2 3 4

Points where D , o and f are given a

Vertices Cell centers

Cell centers Vertices

Table 2.1. Points where D . o and f are to be given.

3. CELL-CENTERED DISCRETIZATION, CASE 1

The cell whose center is at (i-l/2)h (j-l/2)h2 is called fl.. (1 < i < N j , 1 < j < Nj)

and is depicted in Figure 3.1. This cell is also called a finite volume. The approximation of 4> in the center of this cell is called <j>..\ when this is convenient this is also regarded as an average value of <f> in fl...

Integrating (1.1) over fl.. gives, using Green's theorem,

2

f E D r ^ n dr + f o^df) = f

Jan.,

ir,

ad

*

a

Jn.. Jn..

ij a~l a ij ij

fdfi, (3.1.1)

2 D A 1 I IP * -1

1*

C F B

Figure 3.1. Finite volume (ï..; P has coordinates ((i-l/2)h ,,(j-l/2)h?).

3.1. Interior finite volumes

a. The line integral. We have

■f E D ^ n d r = [ B D , ^ d x , - fC D, p- dx- + fD D , # ^ dx

Jan., „ t i adxn a J A 2ax 1 J

B

i dx 2 J

c

2 ax

A

J DDI

^ d x ax, a

2-The integral along AB is approximated by:

I A S ^ M E

0

2 <B> a ^2 d* , -T( D2( A ) + D2( B )^ i j - V l) / h2 = (3.1.2) = s . W . . .(<)>.. - <j>. . . ) , 1 i , j - 1 i j i . j - 1 with s. = h . / h - , and: W2. = -\(D. .. + D , . , .), ij 2 2 , i j 2 , i - l , j

with D_ .. the value of the function D . in the point (ih. ,jh«).

Note that in this case the diffusivity at the interface is an arithmetic average of diffusion coefficients in the adjoining cells.

The other faces of f}., are treated similarly.

ij

b. The surf ace integrals.

These integrals are approximated as follows:

J n . . I 2 i j i j (3.1.3)

and

f fdfl = h.h,f.., (3.1.4)

where o., (f..) is an approximation of o (f) in the point ((i-l/2)h ,(j-l/2)h_), given by

o.. = - j (o(ih,,ih2) + o ( ( i - l ) hr( j - l ) h2) + oOhj.ü-Dhj) + a ( ( i - l ) h j , j h2) )

and similarly for f... >J The resulting stencil is:

-s.W.1 2 ï-

U

-s,W2. 1 ij sum + h,h„<7.. 1 2 i j -S1WÜ - 1 -s-W.1. 2 ij (3.1.5) ,1 1 where s_ = h . / h . , W'. = -^ (D. . . . + D . ..) and 2 2 1 ij 2 v l,i,j-l 1,IJ sum = s. ( w2 . . + W2.) + s, ( w ! . . + W.'.l. 1 v i , j - l iy 2v ï - l , j ij-*

3.2. Boundary finite volumes

The sides of fl are labeled r . , . . . , r . as in Figure 3.2.

Figure 3.2. Nomenclature for the sides of n.

When n.. has a face on the boundary, the line integral in equation (3.1) has to be treated differently.

a. Dirichlet.

Suppose AB (Figure 3.1) is on T.. We have <j>v = g

The line integral is approximated by

-B

"2

D . f*- dx. = s . w r . .(#.. - 0_)

i ^ 2 9x0 1 I i,j-l ij E (3.2.1)

and similarly for other parts of 3fl. The resulting stencil is:

-s-W.1 . . 2 i - l , j

-s.W?.

1 ij sum + h.h-cr.. 1 2 ij 0 -s W.. 2 ij _ (3.2.2)

b. Neumann or Robbins.

Again, let AB (Figure 3.1) lie on r . . We have

(-§£

_an

_T + c

_^

«0

_J

_E

F

" Sp­

_Ö

_E

-The following approximation is made:

(3.2.3)

h - *? - 1

h

2

ak

Hence

E

° \ V * "

C<é)

E-

(3.2.4)

^E=^

h

2

g

E

+

^ / (

1 +

i

h

2

C

E ^

(3.2.5)

The line integral is approximated as follows:

f

D

-M

d x

„

h w

2 J± „

J

A D

2 3 x

2 d X

l -

h

l

W

i , j - 1 3 x

2 E

-• » -

h

l

W

r j - l t e E -

c

E ^

By substitution of (3.2.5) one obtains

f D„ -^- dx, ~

2 ^

d x

l " 2 7 1 ^

W

?,j-l

(c

E*ij -

8

E>-The resulting stencil is:

(3.2.6) (3.2.7) -s2w; - . , j -s.W2. 1 i j 2 h

l

C

E 2

s™ + h

1

V

u + r r S

^ ;

w

ij-i

-s,W.. 2 i j (3.2.8)

where sum = s„ I W. , .+ W..J +s,W...

2 i - l , J ij 1 ij

4. CELL-CENTERED DISCRETIZATION, CASE 2

Consider the finite volume Cl.. given in Figure 3.1. Now the diffusion coefficients D may be discontinuous along 80.. and are approximated by constant values D . . i n

a ij a,ij

Ct... Again, we have to approximate the integrals in (3.1).

4.1. Interior finite volumes a. The line integral.

The integral along AB is approximated by:

.B

I

A

°2 £2 d x

l =

h

,

D

2,i^ij - ^ I ^ = 2s,D

2

..(*.. - <»

E

). (4.1.1)

In order to obtain <j>r we use continuity of <j> and D . -r-^-:

2,ijvyij ^E-1 2 , i , j - r * E *i,j-l

Hence ^_ = ( D . .4.. + D . . . , f . , ) / ( D , .. + D , . . , ) , E 2,ijrij 2 , i , j - lTi , j - r 'v 2,ij 2 , i , j - l " so that (4.1.1) becomes: .B f D_ ^-dx. ~ s . W2. ,(tf.. - f . ,) (4.1.2) JA 2 a x , 1 1 i , j - rri j M , J - I with s. = h / h? and W?. . = 2D- . . D , . . , / ( D . .. + D . . . .). (4.1.3) i,j-l 2,ij 2 , i , j - rv 2,ij 2 , i , j - l '

Note that in this case the diffusivity at the cell face is an harmonie average of the diffusion coefficients in the two adjoining cells.

b. The surface integrals.

The surface integrals are approximated by (3.1.3) and (3.1.4), but now o.. and f.. are defined by:

a .j- a ( ( i - l / 2 ) h1, ( j - l / 2 ) h2) , (4.1.4)

and similarly for f ... 'J

The resulting stencil for the interior points is given by (3.1.5), but with,

W. . . = 2 D , . . D . . . . / ( D . .. + D . . , .)

i - I , j l.ij 1 , i - l,j l,ij l . i - l , / (4.1.5)

and W. . . given by (4.1.3).

4.2. Boundary finite volumes

At the boundaries the line integrals are approximated as follows. Let AB (Figure 3.1) He on the boundary.

a. Dirichlet.

.B

h°

, TT^-dx, = ; 2 s . D , ..((/>.. - g_).

2 3 x . 1 1 2,ij ri j ° E (4.2.1)

The resulting stencil is:

-s,W. . . 2 ï - 1 ,j

-s.W?.

1 ij sum + h. h-a.. + 2s, D-1 2 ij D-1 2,ij -s,W .. 2 l,ij (4.2.2)

b. Neumann or Robbins.

In the same way as in section 3 we obtain:

•B »* 2 h 1 A

f

D 5

f

JA 2 ax DT ;;(cr^;; " 8p)-2 1 _ 2 + h2c£ "2,ijv"Erij BE' (4.2.3)

The resulting stencil is:

-s,W. . . 2 i - I , j -s.W2. I i j 2 hlCE sum + h, h.cr.. + •=—r D , .. 1 2 ij 2 + h-CF 2,ij -s,W.. _{2 ij} (4.2.4)

where sum is defined after equation (3.2.8).

5. VERTEX-CENTERED DISCRETIZATION, CASE 3

As in case 2, D , o and f are approximated by constants in cells tl.., and their values

a ' i j

in n.. are indicated by a subscript ij, as before. But now the unknowns are located in the vertices of the cells, and é.. is the value at the vertex with coordinates (ih. ,jh_). From

TJJ 1" 2

the context it will be clear whether (i,j) refers to a vertex with coordinates (ih ,jh_) or to a cell center with coordinates ((i-l/2)h.,(j-l/2)h.).

5.1. Interior finite volumes

Proceeding as before, the following stencil is obtained:

-*2*Ui

-s.W2. 1 ij sum + h,h„<7.. 1 2 ij - s . W2. . 1 i , j - l -s.W.. 2 ij (5.1.1) with W.1.

= 4

(D,

• , • , + D, . .

.), W?.

= 4

(D.

. . . . + D- . . .)

ij 2 v l,i+l,j+l l,i+l,j" ij 2 v 2 , i + l , j + l 2 , i , j + l ' <7.. = — (o. , . , + o. , . + o. . , + CT..) i j 4 ï + l , j + 1 ï + l , j i , j + l i j

and sum defined after equation (3.1.5).

5.2. Boundary finite volume

Let the finite volume fl.. = ABCD lie at T (Figure 5.2.1). Note that its size is half that of interior finite volumes.

In the Dirichlet case we simply have 0 _ = gp, and #._ is not considered as an

unknown, but transferred to the r i g h t - h a n d side, in order to obtain a symmetrie matrix. In the case of Neumann or Robbins the following approximations are m a d e :

ïl "A

- -z-2- d x . =; - h . W . . ~^-2 dx. 1 1 ij dn E

=-

h

,

w

rj

(

%-

c E

V'

d x . =s - -=■ s - D . . . .(<f>. , . - $..), 2 2 2 l , i + l , jv ri + l , j ij (5.2.1) (5.2.2)

and similarly for A D ; for DC the same approximation is used as for interior finite volumes, of course. F u r t h e r m o r e ,

f . (T0dn=:4 h.h.tr..^..

J n . . 2 1 2 r f i j (5.2.3)

with <T- = "^ (p. . , + o. , . , ) , and similarly for f. The following stencil is obtained:

wiui Pjj 2 i,j+l i + ] , j + r " 2S2Dl , i - l j -s.W2. 1 ij 2 1 sum + h , W . . c _ + -r h . h . c r . . 1 ij E 2 1 2 ij -=• s - D . .. 2 2 l , i j (5.2.4) 1 2 where sum = - r s . ( D . . . . + D , ..) + s,W... 2 2V l , i - l j l , i j ' 1 ij

5.3. Corner finite volumes

Let the finite volume Ci.. = ABCD lie at the corner at the intersection of V, and V.

ij 1 4

( F i g u r e 5.3.1). In the Dirichlet case we simply have <j> = g_., and this value is trans­ ferred to the r i g h t - h a n d side as before.

D .

A Tj B

Figure 5.3.1. Finite volume in corner.

In the case of Neumann or Robbins we have, in general,

t - ! ^ + Cl ^ A = gl,A> ( - ^+ c4 ^ A -g. 4,A'

(g, = g | r , g4 = g | r , c , = c | r , c4 = c | r ). (5.3.1)

! " 4 1 4 The line integrals are approximated as follows:

Tt

f D. -i$- dx, ~ 4 ",D-,(C)

3

t

( A )

- 4

h

.

D

-, ■ , ■ .(C. A^-- - g. A). (

5

-

3

-

2

)

J A 2 d x . 1 2 1 2 3 x . 2 1 2,i+l,j+r 1,A ij 1,A

f D, -j^- dx, ~ 4 h-D,(C)

3

ff

A)

~ 4

h

->

D

, • ■ • i(c, A-- - S

A

A)» (

5

-

3

-

4

)

J A 1 3x. 2 2 2 r 3x. 2 2 l , i + l , j + r 4,A ij °4,A

J B

D

l ^

d x

2 * T « 2

D

l . i

+

l J

+

I < V W

- f D, - ^ - dx. = 4 s,D, . . . M.. - <j>. . .). JD 2 d x . 1 2 1 2 , i + l , j + l i j i,J+l (5.3.5) (5.3.6)

The surface integral is approximated by

IJ

r 4 1 2 i + l , j + lri j

(5.3.7)

" 2 Sl ° 2 , i + l , j + l

0,sum+-rh,h,<7. , . , + —h,D- . . . , c , j + ^ h - D . . . . , c . . , - - r s . D . . . 4 1 2 i+l,j+l 2 1 2,i+l,j+l I,A 2 2 l,i+l,j+l 4,A 2 2 ],i+l,j+l

(5.3.8) w h e r e s u m - - i ( s1D2 + s2D,).+ 1 J + 1.

6. VERTEX-CENTERED DISCRETIZATION, CASE 4

As in case 1, D , o and f are approximated by constants in cells whose centers are the vertices; indices (ij) refer to the vertex with coordinates ( i h . J h . ) . As in case 3, the finite volumes are indicated with dashed lines in Figure 2.2.

6.1. Interior finite volumes

Proceeding as before, the same stencil as in section 5.1 is obtained, but with

W.1. = 2D, ..D, . . ./(D. .. + D. . . .), (6.1.1)

ij l,ij l , i + l , j l,ij l,i+ l , j "

W?. = 2 D . . . D , . . , / ( D , .. + D , . . . ) , (6.1.2) ij 2,ij 2,i,j+l'v 2,ij 2,i,j+r'

a.. = a... (6.1.3)

ij i j

6.2. Boundary finite volumes

Let the finite volume fi.. = ABCD He at r . , as in Figure 5.2.1. In the Dirichlet case,

ij 1

the unknown <f>- = gF is eliminated. In the Neumann or Robbins case, the following

approximations are made: .B

f

D

M

dx

„ _

h D

_M

T h e following stencil is obtained: - S - W . , . 2 ï - l , j - s . W2. 1 O sum + -=■ h . h - a . . + h.CT-D, .. 2 1 2 ij I E 2 , i j -S-W.. 2 ij (6.2.2) /ith sum = s , (w.1 , . + W.'.l + s.W2.. 2 v i - I , j ij-* 1 ij

6.3. Corner finite volumes

Let the finite volume fl.. = ABCD be situated as in F i g u r e 5.3.1. In the Dirichlet

ij

case, the unknown <t>. = g . is eliminated. In the Neumann or Robbins case, we have (5.3.1). The line integral along AB is approximated by:

f D,£tdx.^±h.D...-?4^- = ±h.D...(c. J..-ë. .). (6.3.1) J A 2 dx_ 1 2 1 2 , i j d x - 2 1 2,ij 1,A ij 1,A

Similarly, . D

D,

P-

dx,=

4

^D,

..

- ^ -

= 4

h,D. ..(c,

.

#..

-

g. . ).

A 1 3 x . 2 2 2 l , i j d x . 2 2 l , i j 4,A*ij ° 4 , A ' (6.3.2) F u r t h e r m o r e , r C

T

D

JB d<t> i , „ i , . _ d x . = -^ s , W . . ( 0 . . - 4>. , . ) , 1 dx 2 2 2 i jv ri j ï + l , j (6.3.4) fC 2 d x - 1 2 1 i jk ri j i,J+l (6.3.5)

The following stencil is obtained:

1 «;2

- -z s,W..

2 1 ij

0 sum + -7 h . h - a . . + -r ( h . D - ..c. . + h - D , . . c . . ) - -=• S-W.. 4 1 2 ij 2 1 2,ij 1,A 2 1 .ij 4,A' 2 2 ij

(6.3.6) with sum = -r (s.W + S-W J ...

2 v 1 2 ■* ij

REFERENCES

[1] Alcouffe, R.E., Brandt, A., Dendy, J.E. and Painter, J.W.: The multigrid method for the diffusion equation with strongly discontinuous coefficients. SIAM, J. Sci. Stat. Comp. 2 (1981), 430-454.

[2] Behie, A. and Forsyth, P.A.: Multigrid solution of the pressure equation in reservoir simulation. Proceeding of the sixth SPE Symposium on Reservoir Simulation, pp. 29-42. SPE of AIME, 1982. Paper SPE 10492.

[3] Dendy, J.E., Jr.: Black box multigrid. J. Comput. Phys. 48 (1982), 366-386. [4] Kettler, R.: Analysis and comparison of relaxation schemes in robust multigrid and

conjugate gradients methods. In: W. Hackbusch and U. Trottenberg (eds.): Multigrid methods. Proceedings, Köln-Porz, 1981. Lecture Notes in Mathematics 460: 502-534. Springer-Verlag, Berlin, 1982.

[5] Kettler, R.: Linear multigrid methods for numerical reservoir simulation. Ph.D. Thesis, University of Technology, Delft, 1987.

[6] Khalil, M. and Wesseling, P.: Cell-centered and vertex-centered discretization of interface problems. Report 88-18, Dept. of Techn. Math. and Informaties, Univer­ sity of Technology, Delft, 1988.

[7] Wesseling, P.: Cell-centered multigrid for interface problems. J. Comp. Phys. 79 (1988), 85-91.

CHAPTER III

VERTEX-CENTERED AND CELL-CENTERED TWO-DIMENSIONAL MULTIGRID METHODS

FOR INTERFACE PROBLEMS

Page

1. INTRODUCTION 3.2 2. CELL-CENTERED AND VERTEX-CENTERED COARSENING 3.2

2.1. Vertex-centered coarsening 3.3 2.2. Cell-centered coarsening 3.4 3. EQUATION AND STENCIL NOTATION 3.5

3.1. Equation 3.5 3.2. Stencil notation 3.5 4. MULTIGRID ALGORITHM 3.7

5. TRANSFER OPERATORS 3.8 5.1. Transfer operators for MGVC 3.8 5.2. Transfer operators for MGCC 3.13

5.2.1. Restriction type D 3.13 5.2.2. Restriction type N 3.13 6. COARSE GRID OPERATORS 3.17 7. NUMERICAL EXPERIMENTS 3.23 8. STORAGE AND PRELIMINARY WORK 3.33

9. CONCLUDING REMARKS 3.34

REFERENCES 3.34

l.INTRODUCTION

Many methods have been proposed to solve large systems that arise from discretizing elliptic partial differential equations. A class of iterative methods which currently plays an important role in engineering practice, such as reservoir simulation, is formed by SOR and iterative methods based on incomplete factorization accelerated by Lanczos type methods. A new class of iterative methods is that of multignd methods (MG), which use a sequence of increasingly fine (or coarse) grids [3,6].

MG methods have two main components: coarse grid approximation and smoothing. In problems with interfaces, vertex-centered coarsening together with Standard bilinear interpolation leads to inaccurate coarse grid approximation, resulting in deterioration of the rate of convergence. In [1,2,4,5,6,7,8,9], vertex-centered MG methods with matrix-dependent transfer operators are presented.

Such methods have a good rate of convergence but need much storage and preliminary work. A new method has been developed in [11,12,13], based on cell-centered coarsening and simple interpolatory transfer operators.

In this chapter cell-centered and vertex-centered multigrid methods for solving interface problems are compared. These methods differ in the location of the nodes in the grids and in the transfer operators. Numerical results of several test problems with strong discontinuities in equation coefficients are presented. We compare also storage and work required by each method.

2. CELL-CENTERED AND VERTEX-CENTERED COARSENING

There are basically two ways to construct coarse grids in multigrid methods. Suppose that we have a two-dimensional domain subdivided in cells of equal size. Suppose that we have n cells, a = 1,2 in each direction. The unknowns may be located at vertices of the

a

X X X X

X X X X

X X X X

X X X X

Figure 2.1.a. Vertex-centered grid. Figure 2.1 b. Cell-centered grid.

In the vertex-centered case the computational grid G is defined by:

G = | ( x , , x - ) : x = i h , i =0,l,...,n , h = L /n , a = l , 2 V 1 1 2 a a a a at at a at J

In the cell-centered case G is defined by:

G = l ( x . , x . , ) : x = (i - 4 ) h , i = l,...,n , h = L /n , a = 1,2V

1 1 2 a v a 2J a a a a a a )

To each case corresponds a coarsening which characterizes the corresponding multigrid method:

- Vertex-centered multigrid (MGVC) associated with vertex-centered coarsening. - Cell-centered multigrid (MGCC) associated with cell-centered coarsening.

In this section we need to consider only two grids, called G and G; • indicates a coarse grid quantity. Let n be even.

2.1. Vertex-centered coarsening

The coarse grid G is selected as follows:

G = ( x = ( x . , x . ) : x = i h , i = 0,l,...,n , h = 2h , n = - y - V V. 1 2 a a a a a a a a 2 )

This means that G is obtained from G by deleting every other vertex in each direction (cf. Figure 2.2.a).

r i

^

-3

« H

_-a-

-i)

®-

_-e-

_-«

Figure 2.2.a. V'ertex-centered coarse grid.

Figure 2.2.b. Cell-centered coarse grid. 2i-l,2j 2 i - l , 2 j - l 2i,2j 2i,2j-l Fine

J

t _ ,

Coarse

Figure 2.3: Cell numbering in two dimensions.

2.2. Cell-centered coarsening

The coarse grid G is selected as follows:

. _ n _ . G = -I x = (x. ,x_): x = [i - -T) h , i = l,...,n , n = —z-, h = 2h \.

L 1 2 ' a K a 2y a a a a 2 a a)

G is the set of coarse grid cell-centers where each coarse cell is the union of four fine cells (cf. Figure 2.2.b and 2.3).

3. EQUATION AND STENCIL NOTATION

3.1. Equation

We consider the following two-dimensional diffusion equation with discontinuous and/or anisotropic coefficients:

2 - . . , 2 . ^ . J + (70 = f, x E n = , dx v a öx J n-\ a=\ ct at U _ I £ -ir ( D # ) + ^ = f- x s n = n ( 0 , L ) , (3.1.a) D > 0, CT> 0,

and boundary condition

b - f ^ - + c 0 = g, b > 0, c > 0, b + c * 0, on 3D (3.1.b) on

where n is the outward normal to the boundary.

D , f, g, CT, b and c are given functions. D and CT may have strong discontinuities across

Ce ct

internat interfaces. This equation is similar to the pressure equation in oil reservoir simu-lation. We use finite volume discretizations as described in [10,11,13]. The resulting system on the finest grid nas a matrix A with five non-zero diagonals.

3.2. Stencil notation

We will make use of the following stencil notation. Let the grids be G ,G ,...,G , with G the finest grid and G the coarsest grid. Let the space of grid functions on G be denoted by $ :

*k: Gk -* K.

k k k k k

With A : * -+ $ , the stencil representation of A 0 is defined by:

( AkA = £ Ak( i , j ) 0k. , i £ Gk, (3.2)

• ™d J

j e E

The set of values:

S

k

= / j e Z

d

: A

k

(i,j)#0, f o r s o m e i e G

k

\ ,

(3.3)

k k k

is called the structure of A . The stencil [A ] of A is defined in the usual way. For example, if A has the familiar 5-point structure, we have:

kl i = Ak(i

-v

Ak(i,e2) Ak(i,0) Ak(i,-e2) 1 Ak(i,ej)

-where e. = (1,0) and e . = (0,1). k k-1 k k k-1 With R : $ —► * a restriction operator the stencil representation of R 0 is defined as:

( R

k

*

k _ 1

) . - z R

k

a,j) ^

k

: ' .

j 6 Zd

(3.4)

Define the inner product on $ by:

i e Zd

(3.5)

If i $. G the corresponding grid function values are defined to be zero. k * k k-1 With this inner product, one finds for the adjoint R ' : $ —► $

( R

k

' * A - £ R

k

(j,i-2j) rf

k

k k+1

From the preceding equation follows the stencil notation for P <f> for prolongation _k «k+1 , k operators P : * —» * : (3.6) „k.k+K , k + l ( P T ); = T. PK'*(j.i-2j) *j . j e

z

d (3.7)

4. MULTIGRID ALGORITHM

Let the discrete fine grid problem to be solved be denoted by

A«S = f. (4.1)

We present a simple well-structured FORTRAN description, using only one goto, of the fundamental MG algorithm (correction variant, V- or W-cycle), see [12]. Some notations:

M : the number of grids used;

k : grid index (k = 1, finest grid, k = M, coarsest grid); S : smoothing subroutine;

f : current right-hand side (= restriction of residuals on coarse grids); $ : correction (k > 2) or solution (k = 1);

ic(k): counter of multigrid cycles, ic(l) = maxit, number of MG iterations desired; 7 : type of cycle (T = 1: V-cycle, 7 = 2: W-cycle).

The FORTRAN description:

Subroutine MG(0,f) k = 1 ic(l) = maxit 10 if (ic(l).gt«0) then if (ic(k>eq-0»or-k-eq«M) then if (k-eq-M) then call S(*M,fM) ic(M) = ic(M) - 1 endif k = k - 1 ,k ,k „k,k+l call S ( ^ \ fk) ic(k) = ic(k) - 1

else call S ( / , fk) fk + l = Rk + l( fk _ Ak^kj k = k + 1 0k = O ic(k) = 7 endif goto 10 endif end. 5. TRANSFER OPERATORS

5.1. Transfer operators for MGVC

In vertex-centered multigrid for interface problems, matrix-dependent transfer operators have to be used. Several definitions of such operators has been given [1,2,4,5,6, 7,8,9]. Here we restrict ourselves to the 'collecting' matrix dependent operators proposed in [1,2]. They can be conveniently described by the stencil notation introduced in section 3.

P<t> is represented as (cf. (3.7)):

(P« = £ P*(j,i-2j)*., i € G , (5.1)

J E Zd '

where we drop the superscript k for convenience and indicate coarse grid quantities by an overbar. P is defined by the stencil [P ]..

Let e. = (1,0), e. = (0,1) and e, = (1,1) (= e. + e-). The vertex i e G coincides with the vertex 2i € G (Figure 5.1).

We define:

In other words, we have interpolation between nearest neighbours only. Using (5.1) we see immediately that

P*(i,0)= 1. (5.3)

2i+e„

!f__ ,

2i+e3

2i+e,

Figure 5.1. Indices of grid points.

For the grid points 2i + e , a = 1,2, P is defined as follows. Let [A]. be given by:

[ A ] j

-A(i,-ej+e2) A(i,e2) A(i,e3)

A(i,-ej) A(i,0) A(i,ej)

A(i,-e3) A(i,-e2) A(i,ej-e2)

(5.4)

Let C a matrix which is A or a slight modification of A as defined below. By summing the columns of [C]. we obtain the operator H defined by:

H(i,(a,0))- E C(i,(a,/?)) a = -1,0,1.

By summing the rows of [A]. we obtain the operator V defined by:

V(i,(0,/J))= E C(i,(a,/3)) a=-l

For grid point k = 2i + e. we define P<j> by:

P= -1,0,1.

(5.5)

(5.6)

which gives

( P «k = - C E H(k,(Q,0)XPk , ]/H(k,0). K a*0 K + ( a'U )

Using (5.1) and (5.2) we have

W )k " " C S H(k.(o,0)) E P*(U+(a,0)-20 * , ] /H(k,0) K a#0 * ' = E (" E H(k,(a.0))PV,k+(a,0)-2O) *,/H(k,0) £ a*0 E PV,k-2*) % . (5.8) e c Hence pV,k-2£) = - Y. H(k,(a,0))PV,k+(a,0)-2/)/H(k,0). a*0

For k = 2i + e . , l can take only the values i or i + e.. For l - ï, we obtain:

H(2i+e.,-e.) P*(i,e.) = - — (5.9) 1 H(2i+e,,0) and f or l = i + e . : H(2i+e.,e.) P*(i+e,,-e,)= — — . (5.10) 1 • H(2i+er0) We deduce: H(2i-e e ) P ( i , - e , ) = - (5.11) 1 H(2i-ep0)

Using the same method we deduce for the grid point 2i + e»;

V(2i+e -e )

P*(i,e,) = - * — * - (5.12)

2 V(2i+e2,0)

*

V(2i-e2,e2)

p (i.-e,) = ' ■ (5-13) ₂

V(2i-e2.0)

For the points k = 2i + e , , we define P<£ such that:

( C ( P « )k- 0 . (5.14) We obtain:

W V = " C E C(k,j)(P« . ] /C(k,0). (5.15)

k j # 0 k+j Using (5.1) we obtain: W V " " E ( E C(k,j)P*(U+J-2i)) ü,/C(k,0). (5.16) k £ j*0 l This gives: P*(t,k-2l) = - Y. C(k,j)P*(/,k+j-2Z)/C(k,0). (5.17) }*0

For k = 2i + e-, £ can take the values i, i + e . , i + e- and i + e. only. For t = i we obtain: Y, C(2i+evJ)P*(i,e.+j) P*(i,e.) = - - ^ . (5.18) 3 C(2i+e3,0) For t = i + e . : E C(2i+e3,j)P*(i+e1,e2-e1+j) P*(i+e,,e--e.)= - - ^ . 1 2 1' C(2i+e3,0) We deduce by using i := i + e . :

£ C(2i-e1+e2,j)P*(i,e2-e,+j)

P * ( i , e , - e . ) - - -1*5 : • (5.19)

Similarly we obtain: r » * / -Y C(2i+e]-e2,j)P*(i,e1-e2+j) PT(i,e. - e , ) = - - ^ (5.20) 1 2 C(2i+ e ]-e2,0) Y C(2i-e3,j)P*(i,-e3+j) P*(i,-e.) = - - ^ . (5.21) 3 C(2i-e3,0)

Let for e = (a.,a_), ||e|| = | a . | + | a - | . Then (5.18-5.21) can be summarized, for e such that ||e|| = 2, as follows:

Y C(2i+e,j)P*(i,j+e)

P*(i,e) = - j :' I J+ el ' * ' . (5.22)

C(2i+e,0)

The elements of P occurring in the right-hand side are defined by (5.2), (5.3), (5.9), (5.11), (5.12), (5.13).

The restriction operator will be defined by:

R = P*. (5.23)

We will choose C in two ways:

1) C = A (5.24)

as proposed in [1,2,4,7,8,9].

The transfer operators following from this choice will be called type A.

2) C(i,j) = A(i,j) for j # 0 (5.25a)

C(i,0)

£ A(i,j) if I £ A ( i , j ) / £ A(i,j)| < 10" j*0 ' j j*0 '

A(i,0) otherwise

(5.25b)

where p is an integer to be chosen.

The transfer operators related to this choice of C will be called type B.

5.2. Transfer operators for MGCC

5.2.1. Restriction type D

The restriction type D is the scaled adjoint of linear interpolation in triangles with zero on the boundaries ([11,13]). This is appropriate for Dirichlet boundary conditions. But we will also use this restriction for other boundary conditions. The operator R has the following stencil: [ R ]; = l 16 w. w. 1 0 0 w. 1 w.+2 2 0 0 2 e.+2 i e. i 0 0 e. e. ï 1 0 - 1 - 2 (5.26) -1

where w. = 0 in cells along the west and north boundaries, and e. = 0 along the east and south boundaries: otherwise, w. = 1 and e. = 1. The indices of the rows and columns of

ï ï

[R]. give j , , j . , respectively, for R(i,j).

5.2.2. Restriction type N

In the interior of the domain R type N is defined as type D: the (scaled) adjoint of linear interpolation in triangles. Interpolation takes place in triangles such as HK.F and

HEF (Figure 5.2). H E X X X X K F

Figure 5.2. Interpolation for ad joint of restriction, x .• center of fine cells. , ; center of coarse cell.

A boundary modification is obtained by the following considerations. Suppose we have a Neumann boundary condition only and that o = 0 in (3.1.a), so that the matrix A is singular. A base of Ker (A) is the grid function e, which equals 1 everywhere. Let e also

T

be a base of Ker (A ). The fine grid equation to be solved

A<j> = f (5.27)

has a solution cj> only if

( f , 0 = O (5.28)

which is assumed to hold. Consider the two-grid method, indicating coarse grid quantities by an overbar. The coarse grid equation to be solved is

A <j> = R ( f - K<j>) (5.29)

with <j> an approximation to <j>. Our prolongation (see below) has the property

Pe = e (5.30)

so that

Hence A is singular, and (5.29) has a solution only if it is consistent. This is ensured by choosing R such that

T

-R e (5.32)

As a consequence we have

—T T T T

-A e = P -A R1 e = 0 (5.33)

so that the consistency condition for (5.29) is

(Rf - R A f e) = 0 (5.34)

which is satisfied because (Rf J ) = (f,R e) = (f,e) = 0 and (RA^.ë) = (tf,A e) = 0. To (5.32) we add the condition

R e = e (5.35)

which does not play a role here, but is required for nonlinear problems.

It is found that a suitable boundary modification of the stencil of R satisfying (5.32) and (5.35) is obtained by the following rule: elements to be deleted because they refer to function values outside the domain are added to the nearest neighbour in the "north-west/ south-east direction". This leads to the following stencils at the west boundary, north-west corner and south-north-west corner, respectively:

[R]; = 1 6 0 1 0 0" 0 4 2 0 0 3 3 1 0 0 1 1 , [Rij = ₁₆ 1 " 0 0 0 0 0 4 3 0 0 3 3 1 0 0 1 1

rRi

s

= -ür

0 1 0 0' 0 4 2 0 0 4 4 1 0 0 0 0 (5.36)

In general, R has the following stencil: [R]; = V3,i Vl,i 0 0 4,i w . .+2 2 , i 2 O e3 , i+ 2 O 0 ZA ■ 4,1 (5.37)

The values of e and w , a = 1,2,3,4, depend on the position of point i relative to the boundaries, in the way just discussed.

In multigrid methods P and R must satisfy the following requirement. Let mp- l ,

mD- l be the maximum degree of polynomials that are interpolated exactly by sP or tR

respectively for some real value of s or t. Then we must have ([3,6]):

m_ + m > 2m, (5.38)

with 2m the order of the partial differential equation to be solved. The above restriction has mD = 2. It is therefore sufficiënt here to take an interpolation P with m_ = 1. We

take for restriction R type D or type N and define P by:

[P'L =

1 l

(5.39)

Hence

6. COARSE GRID OPERATORS

In this section we need to consider only two grids. Coarse grid quantities are denoted by an overbar. The coarse grid operator A is defined by:

A = RAP. (6.1)

In the following we give a method for the computation of A in the case of MGVC and MGCC.

Using (3.2) and (3.7) we obtain:

(AP^)j = E A(i,k)(P«.+ k = £ A(i,k) E P*(j,i+k-2j)0 (6.2)

k k j J

and using (3.4)

( A « . = {R(AP^)}. = £ R(i,m)(AP«2.+ m

m

= I R(i,m) E A(2i+m,k) E P*(j,2i+m+k-2j) *.. m

Let j .= i + j , then

m k j J

(A*)= = E E R(i.m) E A(2i+m,k)P*(i+j>m+k-2j) * (6.3)

1 j m k 1+J

He nee:

A(i,j) = E E R(i,m)A(2i+m,k)P*(i+j,m+k-2j). (6.4) m k

In order to determine the storage required, we have to compute the set:

S _ = | j G E : A(i,j) # 0 for at least one i e G V (6.5)

Algorithm structure: A for m e S „ do K for k € S . do — A — for p € S * do P — begin j = ( m + k - p ) / 2 if j € Zd then S _ = — A end. = S U{j} A

Given S _ , the following algorithm compute A: A Algorithm RAP: A(i,j) = 0 for j g S— do A for m e SD do for k e S while (m+k-2j 6 Sp*) do

for i € G while (2i+m € G and i+j e G) do

A(iJ) = A(i,j) + R(i,m)A(2i+m,k)P*(i+j,m+k-2j).

The inner loop does not allow vectorization because of the while clause. We therefore re-place the last two lines by

G j = ■[i e G: 2i+m e G\ (6.6)

G2= l i e G: i+j € G } (6.7)

G3 = G j n G2 (6.8)

for i € G , do

A(i,j) = A(i,j) + R(i,m)A(2i+m,k)P*(i+j,m+k-2j).

The inner loop vectorizes along grid lines.

This algorithm is completely general. It covers cell-centered and vertex-centered multi-grid and matrix-dependent P and R. Notice that in cell-centered multimulti-grid, the

algo-rithm can be simplified because P (i,j) = 1, j 6 Sp*.

To illustrate the computation of G , we give an example in two dimensions. Let G and G given by: G = 4i = ( ' i . ' n ^ 0 5 i. < nxc, 0 < L < n y c | G = 4 i = (i. , i , ) : 0 < i < 2nxc, 0 < i . < 2nyc V i 6 G , is equivalent to: m = ( mrm2) e SR, j = ( jpJ2) € S _ m. m. max ( - J , , - —^~,0j < i. < min (nxc - —y-,nxc-j. ,nxcj

m.

max

(-i,.- -^.o)

< i . < min ^nyc

c»

,nyc-j_,nyc

)•

The subroutine to compute A = RAP in the matrix-dependent case in two dimensions is given in [10].

In the cell-centered case it is simple to give A explicitly. Let A be given by:

i j = "ij 7ij 0 fij "J a.. ij 0 e.. ij

'ü

(6.9)

where (temporarily) i,j e 71. It is found that if R is given by (5.37) and P by (5.39), then A = RAP is given by:

Ji j fi j

's..

i j o.. 0 e . . i j

h

(6.10) where

•1 + (6.11.4)

V C

e

l(

a+

-ï

+

*0

2

i,2j-2

+ ( e

2'

T )

2i

+

l,2j-2

+ 2 ( 0 [ + /

'

)

2i-1.2j-l

+ + ( ( e3 +2) a)2.2. _1 (6.11.1) fi... ( e10 J « ) ) )2 i j 2 j.2+ ( e2 ( a +^+ £ )^ 2 i+l , 2 j - 2+ t < V2 )^ 2 i , 2 j - l + + (e4(«+/9))2.+ 1 > 2 j_1 (6.H.2)

V

2 ( T H ? )

2 i - l , 2 j - l

+

(

w 1

( ^

+

* ) )

2

i - 2 , 2 j

+

C(w

2

+2)-,D

2i

_

1>2j

+

+ (W3Q>2i-2,2j+l ( 6 1 1 3 ) «ij = ( e , ( ^ f ) ) 2 U 2 j_2 + ( e2r ? )2 i + l ! 2 j.2 + W««)2Ï_U2]_X + + ((e3+2)(7+5+r)+f)]2 i 2j.1 + U4(l+l)) 2 i + u j_ + (w 1( /5 + £) ]2 i.2 j 2j+ ( ( w2+ 2 ) ( a + / ? + « + 0 )2 i.] ) 2j ' + 2 ( a+ 1 +« )2 i > 2. + ( w3«2 i_2 i 2 j + 1 + [ w4( a+/ 3 ) ]2. _] 2 j + 1 V( e2f )2 i+l , 2 j - 2 + ^( e3+ 2^ 2 i , 2 j - l + (e4(5+e+(T))2 i + 1 2 j_ j + + 2 ( ^ )2 i 2 j (6.11.5) ï i j j - (wj(»ï+f))2i-2,2j+ t( w2+ 2^ 2 i - l , 2 j+ Cw3(-7+5+r7 +r))2 i_2 2 j +i + + K ^ 2 i - l , 2 j+l ( 6 1 ,-6 ) f ij " ^( w2+ 2 )^ 2 i - l , 2 j + 2 ("+ f )2i,2j + ( w3£ )2i-2,2j+l + + ( w4( 5+ £ +f ) )2 i_1 2 j + r (6.11.7)

As R * P*, symmetry is not conserved in general. Note also that by setting eQ = e,

w = w, a = 1,2,3,4 we obtain A = RAP for R given by (5.26). a

Let us consider the following special case: [A].. = 1 'J r.. 0 ij 8.. e.. i j 'J i j a.. 0 ij (6.12)

with a.. = f. . -y.. = £. . (symmetry) and (a+i+S+e+t;).. = 0 in interior cells. Then we have [A]; f.. 0 i j 'S.. 't.. u n IJ o... 0 i j (6.13) with: «y = C ( e , ( a+ 1 +5 ) ) 2 U j.2 + ( e2e )2 i ; 2. _2 + ^ f ) ^ ^ + + 2 ( a2 i - l , 2 j - l+ Q2 i , 2 j - l Ö/ 1 6' \ = DCwjCa^+S)) 2. _2 2. + ( w202. _2 j 2 j + ( w4f )2. .2 i 2 j + ï y - [ ( e4( * «+f ) )2 i + 1 > 2 j. , + ( e2« )2 i + I ) 2 j.1 + < e37 )2 i + l f 2 j_ , + + 2 ( €2il 2j + e2 i , 2 j - i a /1 6' f i j = C ( w4( f i «+f ) )2. _1 > 2 j + 1 +( w3 l)2. _1 > 2 j + 1 +( w2a )2. _ ,i 2 j + 1 + 2^ 2 i - l , 2 j + «-2i,2j>3/16 -(6.13.1) (6.13.2) (6.13.3) (6.13.4)

If in the interior where w = e = 1 we have (a + 7 + 5 + e + r ) . . = 0, the formulas

a a * ij

above simplify to:

V

( a

2 i - l , 2 j - l

+ a

2 i , 2 j - l > /

8 ( 6 I 4

-

1 )

V ^ i - l ^ j - l ^ i - u /

8 ( 6 1 4 2 ) "£ij = (£2i,2j + e2 i , 2 j - l ) /8 <6 1 4-3)

V

(f

2i-l,2j

+

W

8 ( 6 1 4

-

4 ) Furthermore, we have ( a + 1 + S + € + <")•• = 0. ' v ' * ij

At boundaries and corners, a.., 1... «.., and C.. are obtained in the same way, by

substi-i j ' ' substi-i j ' substi-i j ' si j

tution of the corresponding values of e and w in (6.11.1), (6.11.3), (6.11.5) and (6.11.7) respectively.

It is found that if A is given by (6.12) and if

(a + T + 5 + e + f) „ = 0 (6.15)

in all (interior, boundary and corner) cells, then A is symmetrie with both restriction type D and type N. When (6.15) is fulfilled only for interior points, only restriction type D conserves symmetry.

7. NUMERICAL EXPERIMENTS

In this section we use the vertex-centered multigrid method with matrix-dependent transfer operators and the cell-centered multigrid method to solve three special cases of (3.1). In all tests W-cycles (i = 2) are used. We use one presmoothing and one postsmooth-ing by the symmetrie (forward-backward) point Gauss-Seidel method, or the alternatpostsmooth-ing zebra method (line-Gauss-Seidel, ordering the lines in a zebra pattern, using lines in both directions) when anisotropy is involved. The coarsest grids consist of lxl or 3x3 cells in the cell-centered case and 2x2 or 3x3 vertices in the vertex-centered case. The solution method on the coarsest grid is a direct method using QR decomposition. When on the coarsest grid the system is singular, we take the solution (j> = 0 on a lxl grid; on other grids, we replace one equation i of the system by 4>. = 0. The results of the tests are collected in the tables given later. The quantities listed in the tables are the reduction factors K, defined by

!lr

W

H

2

| r

( ,

"

, )

l l

2

v> 1 (7.1)

where r is the residual at the iteration v (r = f - A<j> ) and the average reduction factor K defined by:

- r

llrM

hy/u

K

= { (0) J ■

{12)

Of course, /c and n depend not only on the multigrid method used, but also on v, f and <t> . However, it was found that in many cases « tends rapidly to the spectral radius. In order to facilitate comparison with other work [4,5] we take v such that K < 10 .In cases where the methods converge well, the value of K thus obtained is not far from the spectral radius.

The discretization used is that described in Chapter II. Two positions of interfaces will be considered:

- case 1: interfaces are midway between two mesh lines, - case 2: interfaces are mesh lines.

Note that in case 1 (resp. case 2) diffusivities across interfaces are the harmonie (resp. arithmetic) average of the diffusion coefficients in the adjacent cells in the vertex-centered case. In the cell-vertex-centered case, arithmetic average is used in case 1 and harmonie average in case 2.

Test problem 1

[11,13].-The domain is the unit square discretized into cells of size h, h = l/n. [11,13].-The diffusion coefficients D in (3.1.a) are given in the following subdomains (Figure 7.1):

n

2 = {(

X

l'

X

2>

n

3

= {(x

1

,x

2

)

0 < x < w, 0 < x , < 1 ]•

ai < x. < w+a, 0 < x . < i l , and

u+a < x. < 1, 0 < x_ < 1 V.

x. = w and x. = w+a are interfaces (lines of discontinuities) with a the width of the band fi„. D is defined by:

2 a D = D = i ' 2 for (x x.) e n /3= 1 and 3 1 0 "1 0 f o r ( xrx2) e n2. a = 1,2 D=2

°1

D=10-10 " 2 D=2 fi3

The right-hand side is f(x.,x_) = x.x«, and o = 0. We take a Dirichlet boundary 2 2

condition (b = 0, c = 1 in 3.1 .b) with g(x.,x-) = x. + x_. The results of the computations in case 2 are given in Tables 7.1.a,b. Results for case 1 are nearly the same as for case 2, and are not reported here. By inspecting Table 7.1 .b we see that restriction type D is better than restriction type N, which is expected because the boundary condition taken here is of Dirichlet type, but both types are applicable. The comparison of Tables 7.1.a,b leads to the conclusion that MGCC is as good as MGVC for problem 1.

h 1/8 1/14 1/18 1/26 1/34 1/50 w 0.5 0.5 0.5 0.5 0.5 0.5 a 1/8 1/7 1/9 3/26 2/17 3 / 2 5 M G V C V 5 6 6 6 6 6 K 0.079 0.105 0.103 0.092 0.102 0.101 K 0.049 0.094 0.092 0.082 0.092 0.092

Table 7.La. Results for problem 1 using MGVC.

h 1/8 1/12 1/16 1/24 1/32 1/48 w 0.6250 0.5833 0.5625 0.5417 0.5313 0.6875 a 1/8 1/6 1/8 1/8 1/8 1/8 MGCCN V 6 8 8 9 9 9 K 0.160 0.205 0.227 0.225 0.229 0.233 K 0.088 0.164 0.175 0.184 0:188 0.192 M G C C D u 5 5 6 5 5 5 K 0.075 0.071 0.135 0.070 0.066 0.065 K 0.061 0.055 0.080 0.048 0.046 0.045

Table 7.Lb. Results for problem I using MGCC.

For test problems 2 and 3 we define the following subdivision of the domain. Let Q = [0,L] x [0,L], L given. Let fl be subdivided in four subdomains fi. i = 1,2,3,4 (Figure 7.2). "3 " l " 4 " 2 0 OJ

Figure 7.2. Subdivision of fi.

Here u designates the location of the interfaces in x or x . direction.

Test problem 2 [4]:

For this problem D. and D_ in (3.1.a) are taken equals to D. f and D are discontinu-ous across the internal boundaries of fl., i = 1,2,3,4 and are defined as in Table 7.2.

ï D f " l 1 0 " 2 1 03 1 n3 1 0

°4

103 1

Table 7.2. D and f for lest problem 2.

o = taken -r=r. The boundary conditions are defined as follows

d<j>

dn for x. = 0 or x . = 0

-£- + ^TFT <\> = 0 for x. = L or x . = L. I. Sn 2D 1 2

Table 7.3 gives results for problem 2 when L = 1 in case 1 for MGVC and MGCC using Dirichlet type restriction (MGCCD) or Neumann type restriction (MGCCN). The smoother used is the symmetrie point Gauss-Seidel method. We see that MGVC is two times faster than MGCCN. MGCCN is slightly better than MGCCD, which is to be expected here.

Table 7.4 gives results for problem 2 in case 2. In this case MGCCN is as fast as MGVC and requires about 25% fewer iterations than MGCCD. For the three methods, the reduc-tion factors are independent of h.

Tables 7.5 and 7.6 give results for problem 2 in case 1 and case 2 respectively for various values of L (cf. [4]). This is a somewhat artificial test problem, since it does not approxi-mate a differential equation, because h does not tend to zero. Table 7.5 shows that for case 1 MGVC is two times faster than MGCCN. MGCCN and MGCCD are about equally fast. In case 2 (Table 7.6) MGVC performs as well as for case 1, but MGCCN and MGCCD did not converge. This non-convergence is may be due to inefficiency of the smoothing on the coarse grids (the Galerkin approximation coarse grid matrices are found to be not M-matrices). It was found that this non-convergence does not occur when o is not discontinuous and large enough.

h 1/8 1/12 1/16 1/24 1/32 1/64 01 0.6875 0.6250 0.5938 0.5625 0.5469 0.5234 MGVC V 7 7 7 7 7 7 K 0.113 0.110 0.105 0.102 0.102 0.102 K 0.114 0.111 0.105 0.104 0.104 0.104 MGCCD V 19 18 19 24 18 20 K. 0.689 0.449 0.453 0.524 0.418 0.437 K 0.480 0.460 0.471 0.558 0.457 0.489 MGCCN V 15 14 16 16 17 18 K, 0.377 0.345 0.401 0.380 0.403 0.400 K 0.382 0.358 0.421 0.410 0.440 0.452

Table 7.3. Results for problem 2, case 1, h = l/n.

h 1/8 1/12 1/16 1/24 1/32 1/64 w 0.6250 0.5833 0.5625 0.5417 0.5313 0.5156 MGVC V 8 8 8 8 8 9 K. 0.142 0.138 0.138 0.138 0.138 0.139 K 0.157 0.161 0.167 0.177 0.178 0.193 MGCCD V 13 11 13 13 13 13 K 0.316 0.257 0.299 0.303 0.288 0.301 K 0.327 0.275 0.319 0.329 0.318 0.342 MGCCN V 7 8 8 8 8 8 K 0.134 0.152 0.130 0.129 0.125 0.122 K 0.126 0.142 0.139 0.143 0.141 0.140

L 8 12 16 24 32 64 u 5.5 7.5 9.5 13.5 17.5 33.5 MGVC V 7 7 7 7 7 7 (C 0.113 0.111 0.106 0.104 0.103 0.103 K, 0.114 0.112 0.107 0.104 0.103 0.103 MGCCD V 16 13 13 12 13 14 K 0.320 0.320 0.259 0.268 0.256 0.306 K, 0.400 0.342 0.315 0.303 0.312 0.353 MGCCN V 13 12 14 13 14 15 K. 0.329 0.284 0.334 0.319 0.336 0.332 K 0.334 0.291 0.351 0.314 0.368 0.376

Table 7.5. Results for problem 2, case 1. h = 1.

L 8 12 16 24 32 64 0) 5 7 9 13 17 33 MGVC V 8 8 8 8 8 8 K. 0.121 0.121 0.123 0.124 0.123 0.121 K 0.139 0.146 0.150 0.157 0.161 0.169

Table 7.6. Results for problem 2, case 2. h = 1. MCCC did nol converge.

Test problem 3 [5]:

This is a two-dimensional version.of the three-dimensional problem (4.4) reported in [5]. For this problem L = 1 and D., i = 1,2 are given in Table 7.7.

Dl D2

°1

1 102 " 2 I Q "2 102 " 3 1 I Q "2 "4 I Q "2 ■

o"

2

Table 7.7. Diffusion coefficients for problem 3.

Problem 3a will refer to the case where unit sources of opposite sign are present at the "south-west" and "north-east" corner cells, and o = 0. Problem 3b will refer to the case -4 where a unit source is present at the "south-west" corner cell only, and o = 10 Problems 3a and 3b have homogeneous Neumann boundary conditions.

Tables 7.8 and 7.9 give results for problem 3a for case 1 and case 2. The rate of conver-gence is independent of h. MGVC is about 50% faster than MGCCN. MGCCN involves roughly 25% fewer iterations than MGCCD. Many different values of w were tried, with similar results.

Tables 7.10 and 7.11 give results for problem 3b in case 1 and case 2. We have the same conclusions as for problem 3a, except that when h goes to zero we need use of transfer operators type B instead of type A for MGVC, since type A gave divergence.

h 1/8 1/12 1/16 1/24 1/32 1/64 (i) 0.6875 0.6250 0.5938 0.5625 0.5469 0.5234 MGVC V 4 4 4 4 4 4 K 0.011 0.010 0.023 0.029 0.033 0.044 K. 0.015 0.011 0.018 0.020 0.021 0.028 MGCCD V 5 7 10 11 11 10 K 0.065 0.200 0.259 0.317 0.308 0.281 K 0.055 0.127 0.237 0.276 0.263 0.247 MGCCN V 6 6 6 7 '1 7 K 0.098 0.111 0.112 0.136 0.147 0.145 K 0.097 0.073 0.093 0.117 0.118 0.115

Table 7.8. Results for problem 3a, case 1.

h 1/8 1/12 1/16 1/24 1/32 1/64 w 0.6250 0.5833 0.5625 0.5417 0.5313 0.5156 MGVC V 4 4 4 4 4 4 K 0.011 0.014 0.025 0.029 0.033 0.046 K 0.016 0.013 0.023 0.024 0.025 0.037 MGCCD V 8 9 11 11 11 10 K. 0.180 0.237 0.303 0.325 0.312 0.292 K 0.164 0.208 0.267 0.281 0.266 0.246 MGCCN V 5 6 6 7 7 6 K 0.050 0.041 0.142 0.123 0.137 0.154 K 0.042 0.071 0.095 0.112 0.111 0.100

h 1/8 1/12 1/16 1/24 1/32 1/64 w 0.6875 0.6250 0.7188 0.6458 0.6719 0.6484 MGVC V 3 3 3 3 4* 4* K 0.003 0.009 0.011 0.015 0.025 0.048 K 0.005 0.007 0.005 0.007 0.017 0.031 MGCCD V 1 8 8 9 9 10 K 0.145 0.175 0.196 0.228 0.273 0.259 K 0.122 0.145 0.160 0.195 0.200 0.222 MGCCN V 6 6 6 6 6 7 K 0.097 0.106 0.121 0.139 0.153 0.134 K 0.085 0.071 0.091 0.089 0.100 0.107

Table 7.10. Results for problem 3b, case 1; *: operator type B used with p = 3.

h 1/8 1/12 1/16 1/24 1/32 1/64 0) 0.6250 0.5833 0.6875 0.6250 0.6563 0.6406 MGVC V 3 3 3 4 4* 4* K 0.007 0.009 0.008 0.019 0.026 0.049 K, 0.009 0.009 0.006 0.012 0.023 0.040 MGCCD V 6 7 7 8 8 8 K 0.088 0.165 0.183 0.222 0.246 0.267 K 0.090 0.132 0.122 0.148 0.167 0.161 MGCCN V 4 6 7 7 6 6 K 0.085 0.055 0.096 0.130 0.131 0.148 K 0.029 0.074 0.104 0.112 0.099 0.095

Table 7.11. Results for problem 3b, case 2; * .■ operator type B used with p = 3.

8. STORAGE AND PRELIMINARY WORK

We present in this section the storage and preliminary work requirements of the multigrid method in the cell-centered and the vertex-centered case. The preliminary work consists here of the computation of the coarse grid matrices and, in the vertex-centered case, of the transfer operators. We consider only the coarse grids, as on the fine grid the storage required is nearly the same in the two cases.

In the cell-centered case, we store 7 reals per node for the matrices, while in the vertex-centered case we need 9 reals per coarse node for the matrices and 8 reals for the transfer operators.

Assuming that the cell-numbers of a rectangular domain are n and n . in the x - and x.-direction respectively, Table 8.1 gives the storage and the preliminary work (operations per coarse grid point) required in the cell-centered and the vertex-centered case. Vertex-centered Cell-centered STORAGE (words) Coarse matrices -3(n1 + l)(n2+l) 7 ~ Tnl 'n2 Transfer operators ~ - j ( n] + l)(n2+l) 0 PRELIMINARY WORK Coarse matrices (a) 169+ 338* 70+ 70* (b) 96+ 120* 63 + 21* Transfer operators 24+ 16* 0+ 0*

Table 8.1. Storage and preliminary work per coarse grid node. (a) Using algorithm RAP of section 6;

(b) Using explicit expressions (6.11.1-7 ) for the cell-centered case and explicit expressions in [2] for the vertex-centered case.