NASA Contr.lctor Report 177973 

lCASE REPORT NO. 85-26 

leASE 

NASA-CR-177973 

19860004486 

Multigrid Solutions to Quasi-Elliptic Schemes 

Achi Brandt 

Shlomo Ta'asan 

Contract No. NASl-17070 

May 1985 I."I,C:~L' 1\t:.~t.kl:C.- CEiJYCR 

U8R'\r.). n,c,s.\ 

, IMv,PTON. ViHGINI>' 

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING 

NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Univers~ties Space Research Association 

NI\SI\ 

NatJonal AeronautiCS and 

Space Administration 

Langley Research Cent. 

Hampton. Virginia 23665 111111111111111111111111111111111111111111111 

NF00720 
1 Report No NASA CR-177973 I 2 Government Accession No 3 Recipient's Catalog No 

ICASE Report No. 85-26 

4 Title and Subtitle 5 Report Dati 

Multigrid Solutions to Quasi-Elliptic Schemes May 1985 

6 Performing Organization Code 

7 Author(s} 8 Performing Organization Report No 

Achi Brandt and Shlomo Ta'asan 85-26 

10 Work Unot No 

9 Performing Organization Name and Address 

Institute for Computer Applications in Science 

and Engineering 11 Contract or Grant No 

Mail Stop 132C, NASA Langley Research Center NASl-17070 

Hampton, VA 23665 13 Type of Report and Period Covered 

12 Sponsorong Agency Name and Address Contractor Report 

National Aeronautics and Space Administration 

Washington, D.C. 20546 14 Sponsorong Agency Code 

505-31-83-01 

15 Supplementary Notes Appeared in Progress and Super- Langley Technical Monitor: J. C. South, Jr. computing in CFD, (E. M. Murman 

F1nal Report and S. Abarbanel, eds) BirkhauseI 

Rn<lt-nn Tnr 1 qR ') nn 235-256. 

16 Abstract -

Quasi-elliptic schemes ar1se from central differencing or finite element 

dl.scretization of elliptic systems with odd order derivatives on non-staggered 

grids. They are somewhat unstable and less accurate then corresponding 

staggered-grid schemes. When usual multigrid solvers are applied to them, the 

asymptotic algebra1c convergence is necessarily slow. Nevertheless, it is shown 

by mode analyses and numerical experiments that the usual FMG algorithm is very 

efficient in solving quasi-elliptl.c equations to the level of truncation 

errors. Also, a new type of multigrid algorithm 1S presented, mode analyzed and 

tested, for which even the asymptotic algebraic convergence is fast. The 

essence of that algorithm is applicable to other kinds of problems, including 

highly indefinite ones. 

17 Key Words (Suggested by Author(s)} 18 Dlstrobutlon Statement 

multigrid, 64 - Numerical Analysis 

FMG solutions 

quasi-elliptic schemes 

Unclassified - Unlimited 

19 Securoty Oasslf (of thiS report) 20 Securoty Classlf (of thiS page) 21 No of Pages 22 Price 

Unclassified Unclassified 22 A02 

N-305 For sale by the National Technical Information SerVice, Springfield, Virginia 22161 NASA Langley, 1985 
Multigrid Solutions to Quasi-Elliptic Schemes 

Achl BMndt and Shlomo Ta'asan 

Department of Applied MathematICS 

Welzmann Institute of Science, Rehovot, Israel 

and 

InstItute for Computer ApphcatIOns In Science and EngIneerIng 

NASA Langley Research Center Hampton V Irgmla. USA 

Abstract QuasI-ellIptIc schemes arIse from central differencIng or fimte element dIscretIzatIOn of elhptlc systems WIth odd order derIvatIves on non-staggered 

gnds They are somewhat unstable and less accurate then correspondmg staggered-grId schemes When usual multigrId solvers are apphed to them, the asymptotIC algebraIC convergence IS necessanly slow Nevertheless, It is shown by mode 

analyses and numencal experIments that the usual FMG algOrIthm IS vel') effiCIent m solvmg quasI-elhptIc equatIOns to the level of truncation errors Also a 

new type of multIgrId algOrIthm IS presented. mode analyzed and tested, for whIch 

even the asymptotic algebraIC convergence IS fast The essence of that algorIthm 

IS applicable to other kmds of problems, mcludIng hIghly mdefimte ones 

Tins research IS supported b} the AIr Force AeronautIcal Laboratones. Au Force Systems 

Command. Umted States Au Force, under Grant AFOSR 84-0070, and also by the NatIOnal 

AeronautIcs and Space AdnumstratlOn under ~ASA Contract No NASI-17070 wlule the 

second author was In resIdence at lCASE, NASA Langley Research Center, Hampton, VA 

23665 
1 Introduction 

2. Definitions and Examples 

3 Instability and Inaccuraty 

Content 

4. Multigrid Troubles and Their Implications 

5. Modified Mode Analysis 

6 FMG Solution to Truncation Level 

6 1 Multlgrid cycle 

6.2 Full MUltlgrld (FMG) 

6 3 Averaging 

6 4 MeasurIng convergence 

7 Algorithm for Fast Algebraic Convergence 

7 1 MultIple coarse grId corrections 

7.2 The modIfied algorithm 

7 3 ModIfied smoothmg analysis 

8. Numerical Experiments 

8.1 The skew Laplacian problem 

8 2 The Stokes and N aVler-Stokes problems 

8.3 Accuracy and stabibty 

8.4 Poor asymptotIC convergence 

8.5 FMG results 

8.6 AsymptotIC convergence wIth new algorithm 

- 2 -
1. Introduction 

QuasI-ellIptic schemes arIse. for example, when central dlfferencmg IS used to 

approXlmate odd-order derIvatIves m ellIptic systems of partIal dIfferential equatlOns, such as the Cauchy-RIemann, Stokes and Navier-Stokes systems Usual 

fimte element approximatlOns to such systems also lead to quasi-ellIptIc schemes. 

Such schemes are m some sense unstable certain highly-osClllatmg components are 

amplIfied in the discretized solutlOn much more than in the dIfferentIal solutlOn 

Instead of the quasI-ellIptic schemes, other discretIzatlOns of the same system 

can usually be constructed whIch are h-elhptic. hence fully stable, and whIch are 

also more accurate than the quasI-ellIptic schemes Sometimes. however, these 

fully elhptic schemes are Inconvement to use In case of ellIptic systems with oddorder derIvatives, for example, full ellIptIcIty is obtained by grId staggerIng. I e., 

by approximatmg dIfferent functlOns on dIfferent grIds (cf [3] and [8]) ThIS IS 

mconvement, espeClally near curved boundarIes. Also the InstabIlIty of qUaslellIptIc approxlmatlOns seldom really hurts, SInce the unstable components have 

very small amplItudes, whIch are still small even m the dIscrete solutlOn The 

inaccuracy IS modest The error In the quasI-elliptic solutlOn IS typIcally tWIce to 

four tImes larger than the error In an ellIptic solutlOn usmg the same grId size 

Thus, qUaSi-ellIptIc schemes are often preferred and are WIdely used 

The mstabIllty of quasi-ellIptIc schemes does seem to hurt when multlgrId 

solvers are apphed' The asymptotIC convergence turns out to be slow, and a 

SImple mode analYSIS traces thiS slowness to the unstable modes. One approach, 

perhaps the best, to deal WIth thIS difficulty IS SImply to ignore It· the algebraiC 

slowness does not matter because It occurs m modes whose amphtudes in the 

algebraIC solutlOn are erroneous anyway, bearmg no relatlOn to theIr amplItudes 

in the true differential solutlOn One should only take care not to Initially admit 

large unstable amplItudes, and to average them out In case they must latter enter 

We show, by mode analyses and numerIcal experIments, that the usual FMG 

algorIthm IS very effective in solvmg quasI-elliptIC problems to truncatlOn level 

(i.e., to the pomt where algebraIC errors are dommated by dl8cretIzatlOn errors) 

SometImes the FMG solutlOn may even be better than the exact solutlOn of the 

discrete equations, because the unstable components of the latter are slow to enter 

Although this is the easIest approach for obtammg fast df.fferentf.al convergence (convergence to the dIfferential solutlOn through a sequence of grIds). another algorithm is presented below which does provIde fast algebra7.c convergence 

for quasi-elliptic schemes ThIS algOrIthm, based on multiple coarse-grId corrections, is interestmg in Its own right, since it IS the simplest example of a new 

kind of algorithms for solvmg problems WIth highly-oscIllatmg solutlOns, mcluding 

hIghly indefimte problems (see [1, §3.2]' [8] and a subsequent article) Smoothmg 

rate analysis, for one quasi-elliptic example, SUItably modified to account for the 

- s -
multIple coarse-grId correctIOns, shows that the new algorIthm should algebraIcally 

be as effiCIent as usual multigrIcl cycles are for fully ellIptIc schemes Numencal 

expenments exactly YIeld these expected convergence rates (see §8 6, tests of such 

algorIthms were also reported 10 [6] ) 

The sIgmficance of the present studIes goes beyond ellIptIc PDE systems 

many non-elhptic systems, such as all subsomc steady-state flow problems. have 

determmants with at least one ellIptIc factor Most discretizatIOns of such systems prOVide quasl-elhptIc apprOXImation to that factor, leadmg to troubles and 

reqUInng cures similar to those reported here 

Moreover. the techmques described 10 thiS artIcle Illustrate the followmg general mulhgnd approaches to general non-elhphc problems (I) Differential not 

algebraiC convergence IS sought and usually easily obtamed Modified methods 

for apnon analyzmg and aposterIorI measurmg such a convergence have been developed (11) WIth considerably more effort. fast algebra'Lc convergence can also be 

obtamed (1l1) The analYSIS of difference schemes. and the derIVatIOn of effiCIent 

smoothers, for any PDE system IS based on the factors of the h-prmclpal part of 

the operator determmant 

We thank Ruth Golubev for some of the calculatIOns reported III §8 

2. Definitions and Examples 

In the followmg Lh WIll represent a system of q real difference operators on 

q gnd functIOns where h, the meshsize of the grId, IS for sImphclty assumed 

to be umform and the same 10 all dIrectIOns That IS, Lh IS a q x q matrIX of 

real polynomials in T}, , Td , TI- I , • Ti 1, where Tt are the grId translatIOn 

operators, defined by 

T~l T:du(;r) = u(~ + !!.h), 

With ~ = (Xl, ,Xd), !!. = (VI, ,Vd) and d be 109 the dImenSIOn of the Euchdean 

space housmg the grid (In case of staggered grIds there may appear non-mtegral 

powers of TJ and Lh will most usually be a matrIX of polynomIals 10 T

J

I

/

2 and 

T- I

/

2 = 1 d) J ' J . 

Three common examples of dIfference operator are 

(I) The fi ve-pomt (com pact) LaplaCIan 

h 1 1 [ fj. = h2 (TO,I + TI,o + TO,-l -t- T_I,o - 4To,o) = h2 1 

1 

-4 

1 

(2 1) 

where TQ ,{3 T] Tg and the array on the left IS the usual pIctOrial deSCrIptIOn 

- 4 -
of the weights of the operator ThIs is the simplest approxunation to the twodImensIOnal Laplace operator ~ -= 02/oxi + 02/0x~ 

(il) The central non-staggered approXImatIOn to the Cauchy-Riemann operator 

h (a} LCR = . -ac 

2 

(22) 

where ac = .l(T. - T-I ) I 2h' , 

(IiI) The central non-staggered apprOXImatIOn to the Stokes operator in two 

dImenSIOns 

a} ) ac 

2 

o 

(2.3) 

For sImphcIty we WIll deal m thIS article only WIth constant-coefficient operators Lh In thIS case the symbol Lh(fl.) of Lh IS defined by 

(1fl.1 ~ 7r) 

for any q-vector A, where fl. = (81, ,8d), fl. ~ = 8I x} + ... + 8dXd and 1fl.1 = 

max(181I, .. , 18dl) Thus, Lh(fl.) IS a q x q matrIX of polynomials in e±I9, , J 

1, . ,d, obtained from Lh by replacing each T, with e,B,. 

Also for simplicity we WIll deal here only WIth homogeneous operators Lh, 

1 e, operators for which all terms in det Lh (the determinant of Lh) have the 

same power in h (ThIS means that Lh apprOXimates a homogeneous differentIal 

operator L, 1 e., det L IS a homogeneous polynomIal In 0/ OXI, •. ,0/ OXd All 

examples above are homogeneous) For homogeneous difference operators, the 

general notion of ellIptICIty measure on a given scale (cf [2, Sec. 3 1] or [3, Sec. 

2.1j) is not needed, and we can use the follOWIng slIDpler definition. 

Definition The homogeneous dIfference operator Lh is ellIptIC of order 2m 

Iff 

d 

I det Lh(fl.) I~ Ch-2m L 8;m for all 1fl.1~ 7r, (24) 

,=1 

where C IS positive and independent of fl.. 

EllipticIty of differential operators is defined In the same way (The parameter 

h IS arbItrary then, and the range of fl. IS unrestrIcted. It is thus more natural m the 

contInUOUS case to replace fl./ h by another phase variable, w = fl./ h say.) It 18 easy 

to see that both ~ and ~ h are second-order ellIptic. Generally, slIDplest central 

approxlIDations to second-order scalar (q = 1) elliptic operators are themselves 

- 5 -
elliptic. But not all central appro2ClmatlOns are For example, the "skew Laplacian" 

1 1 [1 ~x = 2h2 (TI,l + T1,-1 + T-1,1 + T_1,-I - 4To,o) = 2h2 ~ o 1] -4 0 

o 1 

(25) 

or the "long Laplacian" 

2h 1 ~ = 4h2 (T2.0 + TO,2 + T_ 2 ,o + TO,-2 - 4To,o) (26) 

both approximatmg ~, have the symbols 

~ 1 2 2 2 ~ x at) = h2 [(cos 01 - cos O2 ) + sm 01 + sm O2 ] 

~ 2h 1 2 2 ~ (fl) = h2 (sm 01 + sm O2 ) 

which clearly fall to satisfy (24) Indeed. ~ x C7r. 1T) = 0 and ~2h (1T. 0) = ~2h (0, 1T) 

= ~2h(1T, 1T) = 0 Whereas these examples seem somewhat artIficial (although the 

skew LaplaCian does naturally arise in various SItuations, e g., 10 seml-lmphclt 

Lagrange codes [4, § IV] and for some kmds of fimte elements [7]) non-ellIptiC 

operators are very common m apprOXImatIOns to elhptlc systems (q > 1) The 

dIscrete Cauchy-RIemann (2.2) and Stokes (23) operators well represent thlS SItuation· They are the simplest (non-staggered) central approxlmatlans to ellIptIC 

operators, but det L~R = ~2h and det L~ = ~h~2h, hence they do not satISfy 

(24), theIr symbol vanIshmg wherever ~2h does Note that takmg the determinant commutes with passing to the symbol, hence elhptlcity of Lh IS eqUIvalent to 

ellipticity of det Lh, whIch in turn IS equivalent to ellIptICIty of all factors of det 

Lh. 

Fmite element discretlzatlans of the same elhptIc systems. wIth umform nonstaggered partitions, gIve rlSe to simIlarly non-ellIptIC difference operators ThIS IS 

not usually recognized because fimte element dlscretlzatIOns are seldom Founeranalyzed as umform-grid operators 

In all the above examples, even when Lh fails to satisfy (24), It stIll satisfies 

the weaker condItIOn 

d 

/detih(fl) I~ Ch- 2m Lsm2rn (O,), for all/fl/~ 1T, (2 7) 

,=1 

where C is positIve and independent of ~ The term quas1.-eU,pt,c was mtroduced 

in [8J to describe such operators 

Perhaps all reasonable apprOJomatlons to homogeneous ellIptlc equatIOns satISfy (2.7), but for the purpose of including some addItlonal, not-so-reasonable approximatIOns, we can extend the class of operators, and admit any homogeneous 

- 6 -
operators Lh for whlch det ih((l.} vanishes only at a fimte number of points. Thls 

class mcludes for example ~2hx = (T2,2 + T2,-2 + T-2,2 + T-2,-2 - 4To,o)/(8h2), 

whlch satisfies neither (24) nor (2.7), but for which the methods descrlbed below 

are still applicable 

More generally when mhomogeneous operators are also admltted, our methods w1l1 extend to any operator Lh wlth 0(1) "measure of quasi-elliptlcity", defined 

by 

Eh,Q(Lh) = min IdetLh((l.) I /ldetLh((l.') I, Q~lftl~I~'1 

(2.8) 

for some reasonable 0' > 0 Eh.7r is the usual measure of ellipticity Eh described 

in [3J The methods here will in prmClple work for any positive 0', although they 

will gradually deteriorate with the decrease of a: for which Eh,cr(Lh) is still 0(1). 

For clarity, we dlscuss below only homogeneous operators. and the strict qUaslelhptiClty (2 7) is assumed. 

3. Instability and Inaccuracy 

QuasI-elliptIC operators do meet some general stabIlIty reqUlrements even if 

they do not satISfy (24) For example, the skew LaplaCIan (25) is a positive type 

operator, hence satIsfymg the maxImum prmciple. The associated matrIX has a 

dommant diagonal Nevertheless, m a certam sense such operators are not qUlte 

stable. Namely, smce det ih(!D = 0 for some fl ~ 0, in an mfimte space, or under 

periodic boundary conditions, there eXISts a highly-oscillatmg function v h (~) = 

.4 exp( l.fl ~/ h) which satisfies the homogeneous equatIon Lhvh (~) = O. Hence 

the solution, unlIke the correspondmg dIfferential solution, IS not unique (upto 

an additIve constant); it contams an undetermmed hlghly-oscillatmg component. 

Slmllarly, m any bounded domam with any boundary condItions, functIOns wh C~J 

close to vh(~) (e.g. w h = ,,01 v h + ,,02, "oJ being smooth) exist which satisfy the 

boundary condItions and for which Lhwh is everywhere very small. Such wh 

therefore forms an unstable mode. A small change in the equation can introduce 

a large change proportIOnal to w h . This is a kmd of numerical instability, since a 

corresponding large change in the differential solutIOn cannot occur 

This numencal mstability need not hurt much' If the differential system 15 

LU = F and the dlscrete system is LhUh = Fh, all one has to do IS to define 

Fh = ]h F, say, through an averagmg operator ]h which lIquidates the unstable 

modes, I.e jh(O) = 1 and the ratlo jh(fl)/ih(!l.) IS uniformly bounded for all 

1!l.1 ~ E > 0 For example, one can take ]h = sh I'h, where I'h IS any F averaging 

sUltable for the fully elliptic case and Sh is hke the solutIOn averaging Sh described 

below. Even thIS is unnecessary in the usual, smooth case (m the same way that 

the above rule for ]h is frequently neglected for fully ellIptic Lh), because the 

- 7 -
unstabl€' modes. even when un~uly magnified by the dIscretIzation, are usually 

still small 

Generally, the mam dIsadvantage of quasi-elhptic operators IS a certam loss of 

accuracy compared to correspondmg truely elhptic operators, whIch IS sImply due 

to the larger differencmg steps taken in certain terms of the quasI-ellIptic scheme 

In some cases thIS IS particularly obvious, since the grId is locally decoupled into 

severa1 subgrids whIch are not connected to each other by the quasI-elliptic operator For example, the skew-Laplacian (2.5) introduces no coupling between red 

and black pomts (m the usual sense of checkerboard colonng. one color being associated WIth grldpomts where (Xl + X2)/h IS odd, the other WIth even) On each 

subgrid the dIScretizatIon looks like the compact LaplaCIan (2 1) on a rotated grId 

with meshslze hI = V2h SimIlarly, m case of (22). the grId IS decoupled Into 

4 staggered subgnds WIth meshsizes 2h, on each of whIch the operator has good 

elhptlcity (being m fact eqUIvalent to the staggered-grId apprOXImatIOn described 

In [3, §17.2] or [5. §52]) Thus, smce the approXImatIOn IS O(h2). the error in 

case of (25) IS on the average tWIce larger, and m case of (2 2) four hmes larger 

than the errors ill correspondIng fully ellIptIC apprOXImations (assumIng other dIScretIzatIOn errors, related for example t.o the representatIOn of nght-hand sides or 

boundary condItIOns, behave SImIlarly) In these cases, ill other words, each of 

the subgrids can produce the resu1tIng accuracy by itself, other subgrids only add 

work 

When derIvatIves are calculated from the solutIOn, however, the approximating difference quotients may show much greater loss of accuracy, because they 

involve differences between values belongmg to different subgnds The error In 

i-order derivatIves will generally be O(h-l ) tImes the errors m the functIOn Itself 

This excessive error can be aVOIded by takmg dIfferences only from one subgrid at 

a time, or, more generally, by using only dIfference operator Dh such that b h (fl) 

vanishes wherever Lh(fl) does In case of (23), for example: derIvatIves of the 

thIrd unknown functIOn (the pressure) should be approximated by long dIfferences 

such as a;, a; ak. etc. 

The instability described above can also be removed, and the inaccuracy in 

derIvatives proportIOnally reduced, by averaging the solut,on. that IS, by replacing 

the computed solution u h by Shuh, where Sh IS an averaging operator which 

removes all the unstable components In other words, §h(O) = 1 and outside 

a neIghborhood of fl = 0 the ratio §h(fl)/Lh(fl.) should be uniformly bounded 

(wherever defined) For the quasI-elliptic Lh satIsfyIng (2.7) there always eXIsts 

such an averagmg operator of the form 

(3.1) 

WIth mtegral mJ ~ m. On the other hand, the averagmg may further reduce the 

accuracy of the solutIOn. WIth the averaging (31) the lost accuracy IS O(h2) One 

- 8 -
can make that loss O( h2~} by ta.J<mg for example 

Sh = rrd (1 - (1 _ 1T _ 1T-1)II)m, ]=1 2 4] 4] . (32) 

Another slight difficulty tyPlcal to quasi-elliptic approxunatlOns lS the need to 

define extra boundary relatlOns This can satisfactorily be done by extrapolation 

(cf, e.g., §8 2) 

In summary. although quasI-elliptlc dlscretIzations are in prmclple mfenor 

to fully elliptic ones (obtamable for systems by grid staggermg), they can be 

used Smce many programmers conslder grid staggermg a senous complicatlOn, 

especially near general boundanes, quasI-elliptlc schemes and theIr fast solutlOn 

become lIDportant 

4. Multigrid Troubles and Their Implications 

Usual multlgnd solvers Yleld poor asymptotIc convergence rates when applied 

to quasI-elliptIc schemes (see [4J and §8 4 below). The reason IS simple. Slow to 

converge are the unstable modes, such as v h or w h above. They cannot slgmficantly converge by coarse-grld correctlOns, smce they are hlgh-frequency modes. 

essentially mVlslble on coarser levels NeIther can they slgmficantly converge by 

any type of relaxatlOn, smce an error like wh shows a very small resldual function 

Lhwh (compared wlth residuals shown by other modes wlth comparable amphtude) and the correctlOns mtroduced by any relaxation scheme are proportlOnal 

to the size of the reslduals (d. [3. Sec 1 1]) In particlllar, the amplificatlOn 

factor p(!l) of the error mode exp( 7.!l ~/ h) per relaxatlOn sweep must be 1 when 

ih(!l) = 0, and smce the latter equality holds for some l!li = 7r, the smoothing 

factor Ji. = m~/2~1~1~7r Ip(!l} I cannot be smaller than 1. 

The poor asymptotic rates are not a real trouble, though The modes slow to 

converge are exactly those unstable modes for whlch algebralc convergence lS not 

really desired, theIr amplitudes m the algebraic solution being unrelated to their 

amplItudes in the dIfferential solution The only concern lS that these amplitudes 

will rem am sUltably small 

This situatIon IS typ'&cal to all proble~ whzch are not fully elllphc, mcludmg most problems in fluid dynamIcs. Slow asymptotlc convergence of sUltable 

muitigrid cycles occur exactly in those components where not much convergence 

IS needed anyway Whenever thIS sltuatlOn arISes, It is m a sense an absurd to 

try and fix the algorithm (although we show in Sec. 7 below how to do It), smce 

one would then often end up mvesting most of hls human and computer resources 

to obtam lIDprovements whlch are meaningless in terms of solvmg the ongmal 

dt.fferentlal equations. 

- 9 -
Thus, the real ObjectIve of.multIgnd solvers should not be a fast algebraIC 

convergence (convergence of the computed solutlOn u h to the exact dIscrete solution 

U h ), but fast d,&iferenhal convergence (convergence of u h to the true dIfferential 

solution U), usmg any sequence of meshsizes h and measured dIrectly m terms 

of the decrease m II uh - U II as functlOn of the overall computatlOnal work (cf 

13, §13]). ThIS modIfied obJectlye allows for SImpler algonthms, but also calls for 

some modificatlOns in our approach for analyzmg algorIthms, for apnorl predIctmg 

and apostenon measurmg theIr performance The next two sectlOns wIll Illustrate 

these modificatlOns for the case of quasI-elhptic schemes 

5. Modified Mode AnalysIs 

It was shown In Sec 4 that m case of quasI-elhptIC systems Ji ~ 1, but that 

thIS bad smoothmg factor IS not relevant to our real objectIve To analyze a gIven 

relaxatlOn scheme, assume first that It 18 as effiCIent as needed for the dt.iferenttal 

convergence ofthe hIghest frequency modes (whIch should latter be checked by the 

2-level FMG mode analysis mentlOned below) The questlOn then IS what effiCIency 

one should expect from the multigrid cycle (employmg the gIven scheme on all 

levels) in reducmg all other modes As m the conventlOnal smoothmg analysis, our 

SImplifying assumptlOn here WIll be that relaxation on each level should effiCIently 

treat all modes m only one segment of modes, and that the umon of these segments 

should cover all relevant modes Instead of assignmg to grId h the conventlOnal 

segment 7r /2 ~ I~I ~ 7r, however, we can asSIgn to It any segment of the form 

0:/2 ~II ~ II~ 0, WIth any norm II ~ II. That would automatIcally assign to gnd 

h/2 the segment 0:/4 ~II ~ I!~ 0/2, and so on It means that we allow some of the 

hIghest frequency components on any mtermedlate level not to converge effiCIently 

by relaxatlOn on that level, as long as those components effiCiently converge by 

the next-finer-level relaxatlOn ThIS only leaves the hIghest frequency modes on 

the finest grId unaccounted for. whIch IS exactl)- the segment where we do not 

seek simple algebraic convergence Thus. the modIfied defimtlOn of the smoothmg 

factor relevant for our purpose here IS 

(5 1) 

where Jl(~) IS the amphfication factor of exp(1.~ ~/h) per relaxation sweep. and 

the mIDlmum can be taken over all 0 > 0 and over all possible chOIces of the norm 

/I . /I (For a generalizatlOn of thIS definition to cases of semI coarsening. cf 13, 

§12]). 

In case of the skew Laplacian (25), for example. the leXIcographIcally order 

Gauss-Seidel relaxation YIelds the amplification factor 

(5.2) 

- 10 -
so that f.I.( 7r, 7r) = 1 and the conventional smoothmg factor is 1. But choosing 

0: = 7r and II ~ II = max( 181 + 82 j, 181 - 82 1) easily shows that the modIfied factor 

(5 I) YIelds jl S; 5. The same can be shown for the long LaplaCIan (2.6), by takmg 

0: = 7r/2 and II ~ 11= I~: 

In case of systems (q > 1). the q x q amplificatIon matrlX f.I. of the mode 

A exp( tIl- ~I h) depends both or; ~ and on the q-vector A The modIfied smoothmg 

factor jllS then defined by 

ji = mm max Ii f.l.A II / II A II Q/2~1I~II:SQ 

where 0: IS allowed to depend on both ~!!f!.l and A/IAI WIth these defimtIOns and 

sUItable dIstrIbuted Gauss Seidel (DGS) relaxatIOn schemes (see e g [3, §18.6j) 

thIS agam YIelds jl S; .5, for both the Cauchy-RIemann (2.2) and the Stokes (2 3) 

operators In all these cases. stIll better factors are obtamed by four-color ordenng, 

for which defimtIOns (5 1) and (52) should further be extended (cf (32) m [3]) 

As for two-level analyses (cf [3. §4 1] or [5, §4.6]). they always couple lowest 

WIth hIghest frequency modes In non-elhptlc cases some hIghest-frequency modes 

are not expected to converge fast What the analYSIS should then tell us IS how 

effiCIent IS the entIre multlgnd algorithm m reducing the algebraiC errors below 

the truncatIOn errors This can be done by a two-level FMG mode analys,s. whIch 

Founer analyzes the N-FMG algorithm described below (usually for N = 1) by 

assummg exact solutIOn of the coarse gnd equatIOns (both for obtainmg the first 

approximatIOn and m each of the N cycles) and by comparmg for each mode the 

final algebraic error WIth the truncation error (see [3, §7.4]) 

6. FMC Solution to Truncation Level 

Smce the multIgnd cyclmg lS mefficient m reducing unstable mode errors, 

the multignd solver should take care not to start with an IDltial solutIOn whIch 

contams large amplitudes of such errors The overallmItial error in unstable modes 

should better be smaller than the overall truncatIOn error. ThIS lS easily obtamed 

by takmg a first approximatIOn from a coarser grid, employmg interpolation of 

suitable order The usual "Full multIgrid" (FMG, also called "nested Iteration") 

algOrIthm can therefore be used. WIth shght modIfications The usual algonthm 

and ItS modificatIOns are brIefly described in the following For a flowchart, and 

a detaIled dISCUSSIon of FMG algorithms and the order of the first mterpolatIOn, 

see Secs 1.6 and 7 in [3] For SImplICIty we describe here the CorrectIOn Scheme 

(CS) version of the algorithm, so the problems are assumed lmear; It should be 

converted to Full ApprOXImatIOn Scheme (FAS) to treat nonlinear problems [3, 

§8] 

- 11 -
6.1 Multigrid cycle A sequence of grids 18 gIven with meshsizes hk (k = 1,2,3, 

), where hk+l = h k /2 On tne hk grId the dIscrete equatlOns have the form 

(61) 

where Lk approXImates Lk+I GIven u~, an apprOXImate solutlOn to (6.1), the 

multi grId cycle MG for producmg an Improved approxImation u~ 

(6.2) 

IS recursIvely defined as follows 

If k = 1 solve (6 1) by any dIrect or IteratIve method, yieldmg the final result 

uf OtherwIse do (A) through (D) 

(A) Perform VI relaxatlOn sweeps on (6 1), resultmg mane\\' apprOXImatlOn 

(B) StartIng wIth u~-l = 0 make..., succeSSIve cycles 

where 1:- 1 IS a transfer ("reductlOn") of reSIduals from grId hk to grid hk- 1 We 

have used the "full weightmg" 

It-l = ~To.o .... l (To, I + TI,o + TO,-l + T_I,o) 

+ 116 (TI,1 + TI,-l -+- T-I,l + T_I,_J) 

(63) 

(C) Calculate uk = uk + IL1U~-I, where 1;_1 IS a SUItable mterpolation 

("prolongatlOn") from grId hk - 1 to grId hk For problems conSIdered here, bIhnear 

mterpolation is used 

(D) Perform V2 relaxatlOn sweeps on (6 1), startmg WIth uk and YIeldmg the 

final result uf 

The cycle WIth "I = 1 IS called V cycle or V (VI, V2), and the one WIth "I = 2 

IS called W cycle or W (VI, V2 ) 

6.2 Full Multigrid (FMG) The N-FMG IS an algorithm for calculatmg an approximate solution 

u~ = FA1G(k, Fk, N) (6.4) 

to equatlOn (6.1), defined recursIvely as the followmg two succeSSIve steps 

(a) Calculatmg a first approximatlOn u~' If k = 1, put u~ = 0 OtherwIse put 

(6.5) 

- 12 -
where D~_l is an mterpolatlon of solutIOns from gnd hk - 1 to grid hk' and I~-I 18 

some transfer ("averagmg") from gnd k to gnd k - 1. usually a full weIghtIng of 

the type (63). The 1OterpolatIOn DZ- 1 should usually be of order hIgher than that 

of the correction 1Oterpolation 1;_1 mentIoned above [3, §7.1J. In our experIments 

bicubic mterpolatIOn was used 

(b) Improve the first apprpxlmation by N succeSSIve MG cycles 

(J = 1. .. N) 

as defined 10 Sec 6.1. 

6.3 Averag10g The algonthm above IS the conventIonal one, and for equatIOm: 

wIth constant coefficients It reqUIres no modificatIOns In case of quasI-ellIptIC 

equatIons with vanable coeffiCIents. and 10 partIcular in case of nonlmear equatIons. It IS not enough to prevent unstable InItIal errors. because such errors can 

also later be mtroduced due to Interaction between modes It IS then better to 

explICItly reduce the unstable modes by averagmg. such as (3 1) or (32) It may 

also then be Important to replace I;_1uk-1 In step (C) above by I:_ 1Sh./C-luk- 1 

In fact, experIments with non-staggered Navier-Stokes equatIOns (cf Sec 82) 

gave slowly diverging MG cycles unless thIS averagIng was used 

6.4 Measurmg convergence In varIOUS SItuatIOns where algebraIC convergence is 

not attempted, as In the present algorithm and double dIscretIzatIOn [3, §10 2] 

and other algonthms, the questIOn IS raIsed hoV\ to measure convergence, how to 

know, In partIcular, that a solutIOn to the truncatIOn level (1 e , WIth algebraIC 

errors dOmlnated by dIscretIzatIOn errors) has been obtaIned 

The answer is that solutIOn to the truncatIOn level IS not really the Important 

mformatIOn when dtfJerentf.al convergence IS our objective (as It should most often 

be). because (1) SolVIng to truncatIOn level tells us nothing about the truncation error Itself We may for Instance be dOIng good Job In solVIng the algebraIC 

system due to haVIng chosen an easy-to-solve but badly-approxImatIng discretizatIOn (11) A smaller differentIal error may often be obtaIned faster by SWItchIng to 

a finer gnd before the equatIOns on the present grId have been solved to truncation 

level 

The Important InformatIon IS the dIfferential convergence itself as functIOn 

of computatIOnal work ThIS very informatIOn can dIrectly be obtamed from the 

N-FMG algOrIthm Indeed, the sequence of apprOXimatIOns u~, (k = 1,2, ) is a 

sequence convergIng to the dIfferential solutIOn, hence the decrease m the sequence 

of dIfferences 8k =Ii l;-lu7v - u~l II exactly exhIbit the speed of dIfferentIal 

convergence, where the norm II II used to measure 8k can be chosen to exactly 

represent the sense 10 whIch convergence IS sought. One only has to check that 

- is -
the smallness of 61e IS not governed by lack of change from u~ to u~ It IS enough 

for thIS purpose to check that the sUltable residual norm T~ = II Fie - L ku'Jv II IS 

considerably smaller than T~ One can usually also verIf~ that the algebraiC errors 

are below truncatIOn level, e g , by confirmmg that T~ /r~ is considerably smaller 

than 81e /8/c-1 

7 Algorithm for Fast Algebraic Convergence 

Although fast asymptotIC algebraiC convergence IS not needed for fast dIfferentIal convergence. It can stIll be produced by a more mvolved multIgnd algorIthm 

ThIS algOrIthm (also deSCrIbed m [6]) may be interestmg m Its own rIght, smce It IS 

the slIDplest example of a new kmd of algOrIthms (first mentIOned m [1, §3.2]' and 

more fully m [81) for solvmg problems WIth hlghly-osClllator~ solutIOns, mciudmg 

highly mdefiDlte problems 

7 1 Multiple coarse gnd corrections Let fll, fl2. . flf. be all the components for 

which Lh vaDlshes, or. more generally, the centers of all neighborhoods m whIch ~h I L (fl) IS small Usually fl = 0 Then (by [3, §IIj, for example) there eXIsts a 

relaxatIOn scheme WIth fast convergence for all Founer components except those 

close to some flJ The error after few such relaxation sweeps must therefore have 

the form 

l 

Vh(.~) = L ~h(~) exp(zfl] ~/h). (7.I) 

J=I 

where VJ

h are smooth functIOns. Whereas classical multIgnd seeks to apprmomate 

Vh on a coarser grId and the algonthm of Sec 6 approxlIDates V1

h, the new 

algOrIthm will separately apprOXImate each of the ~ ~h, by successively employmg 

l different coarse-grId correctIOns 

Generally, denotmg by H the coarser-grId mesh SIze (H = 2h), the equatIOns 

for VJH, the coarse-grId approximatlOn to VJ

h , should have the form 

(72) 

where Lf(fl) ~ L;(flJ +fl) for small fl, and if!.J(flk) ~ bJ/c (=0 except for bJ] = 1) 

The boundary condItions may couple V, Hand VkH on any pIece of boundary along 

whIch exp(z(!lJ _!lk) ~/h) IS a smooth functIOn There are vanous ways, variational 

ones and more dIrect ones to derIve Lf, I f!.J and the boundary conditIOns There 

also eXIst various ways for solving (7.2). In highly indefiDlte problems the latter 

leads to creating more components on gnd 4h, etc., so that on mcreasmgly coarser 

grids the representation tends to a Founer representatIon 

- 14 -
Ht>rt> we gIve only the very sImple example of solving for the skew Laplacian 

(2.5). (For a more genera] caSe. see [6J) In thIS case t = 2. fll = (0,0) and 

fl2 = (TI, 71"), and one may slIDply take L{f = Lr to be any H-approxlIDatlon to the 

Laplace operator In some sItuations, where the same mechamsm creates both the 

fine gnd and the coarse gnd equatIons, these L~ may agam be skew-Laplacians. 

As transfer operators one can use 

. 2 1] 4 2 

2 1 

1 [ 1 and 1:;2 = - -2 , 16 

1 

-2 

4 

-2 

(7.3) 

Considermg the case thdt the fine-grId boundary condItIons are Dinchlet condItions IdentIcally satIsfied by any fine-gnd approXImatIOn. the coarse-grid boundary 

condItIOns for both V1

H and VZ

H are the homogeneous Dmchlet conditIons For 

solvmg the coarse-grId equatIons (72) the MG cycle of Sec 6 1 can be used, even 

m the case that L~ are themselves quast-elltphc. because. for the purpose of acceleratmg the fine-gnd algebraic conyergence. equatIOns (7 2) need to be solved each 

tlIDe only to theIr truncatIOn level (1 e , only to the level of the error '\;~ H - V]h) 

In case of SImIlar equatIOns but WIth non-constant coeffiCIent, averagmg as m Sec 

6 3 should better be used 

7.2 The modIfied algorithm GIven an apprOXImate solutIon u~ to (6.1), the modIfied multignd cycle MMG for producmg the Improved solution u~ 

u~ -- Ml\lG(k,u~,Fk) (74) 

IS defined non-recursIvely as follows 

If k = 1 solve (6 1) by any duect or iteratIve method. YIeldmg the final u~ 

OtherwIse, perform v relaxatIOn sweeps on (6.1), resultmg m a new approximation 

uk,o, and then. for J = 1,2. . ,l, calculate 

V k- 1,] -- AfG(k - 1, 0, 1~~,~ (Fk - Lkuk,]-l)) 

uk,] = uk,] -1 + exp( ttl] ~! h) 1~_1 v k - 1,] 

with ut = uk,l being the final result 1~_1 agam denotes hnear mterpolatIOn MG 

is the cycle defined m Sec 6.1, WIth a choice of ,/, VI. Vz 

With thIS MMG cycle replacmg the MG cycle, the modIfied FMG algOrIthm 

IS defined m the same way as FMG m Sec 6.2 

7.3 ModIfied smoothmg analYSIS. The smoothmg factor for the above MMG cycle. 

I.e., the ideal factor of convergence one can expect from such a cycle per relaxation 

sweep on the finest grid is defined by 

jl = max IIL(fl) I· 'lr/2:::;1~-~'1 for one], 1~1:::;'lr 

(7 5) 

- 15 -
where \#lUn I is the spectral radms of the amphficatlOn matrIX (or the absolute value 

of the amplIficatlOn factor. If q -= 1) Note that for I ~ 2d , the domam of ~ over 

whIch the maxImum is taken may be empty. In such a situatlOn convergence can 

m prmciple be obtamed wIthout any relaxation on the finest grId ThIS does not 

mean that the algorithm IS more effiCIent than a conventlOnal multlgrid. because 

It employs at least I tImes as many relaxatlOn sweeps on each coarser grId . A more preCIse two-level analysts can of course be made here m the conventlOnal way [3, §4 1] 

For the skew-LaplacIan and the algOrIthm descnbed above, the leXIcographIC 

Gauss-Seidel amplIficatlOn factor (52) attams ItS maxImum (7.5) at (±7l" /2.0) and 

at (±7l" /2. 7l"), Yleldmg jl = 447 

8 Numerical Experiments 

8 1 The skew Laplacian problem Our mam expenmental studIes were conducted 

WIth the skew Laplacian scheme (2.5) m the rectangle {O :::; Xl :::; 2,0:::; X2 :::; 3} 

WIth DIrIchlet boundary conditlOns These conditlOns and the rIght-hand SIde 

of the dIfferentIal equatlOn !:l.U = F were chosen so that the solutIOn U of the 

dIfferentIal equatlOns IS known, to allow direct measurements of dIscretIzatIOn 

errors The sequence of grIds have meshslzes hk = 21- k (k = 1. 2, .), each 

posItIOned so that the boundarIes of n comcide with grId lInes On every level 

Lk IS the skew LaplaCIan, and the relaxatIOn IS leXIcographIC Gauss-SeIdel The 

algOrIthms were those deSCrIbed m Secs 6 and 7. 

Table 1 shows the maxImal dIfferentIal error (maxImal dIfferences between 

computed and dIfferentIal solutlOns) on varIOUS grIds In additlOn, columns headed 

by a or ae show maxImal error m first denvatives, approxImated eIther at gnd 

mldpomt by short difference quotients (the a columns), or at grIdpomts by a~ (the 

ae columns). The upper part of the table gIves these errors for the exact dIscrete 

solutIOn, the lower part - for the solution obtamed by a I-FMG algonthm WIth 

V(2,1) cycles For grId 5 an addItIOnal result (5a) is sometimes given It shows 

errors measured after the solutlOn IS ~veraged by (T1

I

/

2 + TI

- 1/

2

)/2 (cf. (3 1)) 

The table compares skew-LaplacIan with usual (compact) LaplacIan (usmg the 

same meshsIze and the same relaxation), and a case of smooth solution WIth a 

hIghly-oscillatory case The latter IS shown m order to emphasIZe how bad quaslellIptIC schemes can be. In practice such hIghly oscillatory components have very 

small amplitudes If their amplitudes are bIgger than O(h2) (here h~ = .001), 

then second-order apprOXImatIOns cannot be obtamed by any discretizatIOn In 

the hIghly oscillatmg case It was of course necessary to use the full weighting (63) 

for 1;-1 m (6.5), thIS was started WIth k = 7. In the smooth case, however, 

mJectIOn of F was used, in order to obtain a clearer pIcture, clean of F-averaging 

- 16 -
errors 

I 

i 

TABLE u = sm{3x + 2y) u = x{2 - x)y(3 - y) cos 7r(:Z:+Y) 

h h!> I 

1 I 

grid I 

Exact! 3 I 1703 

14 0417 

5 ,.0104 

5a i 0084 

1-FMG! 3 ! 1709 

I I 4 1.0418 I I , I 

: 5 ! 0105 

I ' 

1 5a I 0085 I ! 

~x h ! ~ ~x 1 I I I 

a acl a ! a ae ! 

410 2711 .0517 108 i 2.25 393 3371 0152 

i I I 

115 092 i 0129 033; 225 422 393 i 0038 

031 .027 : .0032 0091608 19494 979 1.0009 

.027 0311 306 979 38 11 I 

436 .276! 0606 151 I 225 393 337 '.0198 I ! I I I ~ 109 .092 I 0169 048, 2 25 4.22 3 93 0050 

I I 

030 027 I 0045 014 I 21 5 752 64 7 i 0014 I I 

.027 .031 ! 2 13 64 7 6 411 

~ 

a 

030 

.008 

002 

.039 

.012 

003 

82 The Stokes and Navier-Stokes problems We have also conducted expenments 

WIth the Stokes operator (2.3). described m detaIl m [3, §18.6] (WIth slIght Improvements, to be described m the new edItion) The unknown grId functIOns of thIS 

operator are Uh: V h and ph - the dIscrete hOrIzontal velOCIty, vertIcal velOCIty 

and pressure, respectively 

In the dIfferentIal problem only velOCities are normally gIven on the boundary. 

In the non-staggered dIscretIzatIOn (2 3) some boundary condItions for ph should 

be introduced (whIch IS a dIsadvantage typIcal to many quasI-elliptic operators) 

For clarIty of expOSItIOn we here aVOld thIS ISsue by showmg results for penodlc 

boundary cond"hons (adJustmg undetermmed additive constants before measurmg 

errors). 

The exact treatment of boundary conditions is Important only m measurmg 

asymptott.c convergence rates It does not much affect results of 1-FMG Therefore 

we will show such results also for the Dirichlet boundary conditIOns In these 

experIments ph at each boundary pomt IS taken equal to the nearest mtenor 

value of ph, and It changes whenever the latter does ThIS does not correspond 

to Neumann boundary condItions, but to couplIng the four subgrids mto which 

the ph grId decouples. A partIal relaxation sweep near DirIchlet boundarIes is 

performed before each full relaxatIOn sweep 

The relaxatIOn employed IS dt.stnbuted Gauss-Bet.del (DGS), a speCIal case of 

a scheme for relaxmg general PDE systems, explained in [3. §3.7J Bnefly, It IS 

equivalent to wrIting Uh = cp~ - alcp~, Vh = cp~ - a~cp~ and ph = _~hcp~, 

- 17 -
and relaxIng by usual Gauss-Seidel thE' resultmg equatIOns m CP~. The changes 

10 the latter Imply changes m 1Jh, Vh and ph, whIch define the actual changes 

performed by thE' DGS relaxatIOn The relaxatIOn ordermg IS 4-colored, relaxmg 

the four mentIoned subgnds one at a tIme 

The domam for this problem IS the square {o S; x] S; 271"} The meshsizes are 

hk = 2-kr. The nght-hand sI~e and the boundary condItions are chosen so as to 

gIve the preSCrIbed solutIOn V = V = P = sm(cos(xl + 2X2»), a perIodic solutIOn 

whIch mcludes many Founer modes The dIscrete rIght-hand sIdes were calculated 

by Fk - 1 = 1;-1 Fk, using (63), starting at k = 8 

Some expenments were conducted wIth averagmg (cf Sec 63). In the 

present case thIS means averagmg of ph only, smce Uh and Vh vamsh m the 

unstable modes When used, thIS P-averagmg employed (3.1), wIth m} = m2 = 2. 

performed on ph m any solutIOn or correctIon Just before mterpolatmg it to a finer 

gnd 

Also mentIOned below are experIments wIth non-stagered mcompresslble 

Navter-Stokes (INS) equat'tons. with procedures sImIlar to those for Stokes For 

detaIls see [3. §19] the modIficatIOn from staggered to non-staggered formulatIOn 

and processmg are the same as for the Stokes equatIOns Results are gIven for the 

Dmchlet problem (U and V given on the boundary, P on the boundary treated as 

above), for the case U = V = P = 1+.2sm(cos(xl +2X2))' We have experimented 

with small and large Reynolds numbers, Re In the latter case anisotropic artificial 

VISCOSity was used m relaxation, its magmtude bemg 14 times the VIscosity mtroduced by upstream differencmg Central differencmg Without artifiCial viscoslty 

was used for the fine-to-coarse reSIdual calculations. allowmg O(h2) solutions to be 

obtamed The large Re PDE problem 18 not elhptic (more precisely, it has small 

ellIptICIty measure), so Its detaIled diSCUSSIOn is beyond the scope here Indeed. 

the present example IS not fully tYPIcal for large Re, because it has no boundary 

layers and no gndhne-streamllne alIgnments 

Table 2 summarIZes four numerical experiments Three with the Stokes (Re = 

0) problem (exact solutIOns for the perIOdIC ("Per") boundary condItions, I-FMG 

solutions WIth W(2, I) cycles for the same problem, and SImilar I-FMG solutlOns 

for the Dmchlet ("Dir.") problem), and one experIment for "mfinite" Re. I e., 

WIth VISCOSIty completely dommated by artificial viscosity The latter experIment 

uses 2-FMG algorithm with W(2,O) cycles, because double discretizatIon (dIfferent artIficial VISCOSIties at different stages) IS mvolved (cf. [3, §1O.2]). For each 

experIment and each grid k, the three numbers shown m the first column are 

max(1I uk - U II, II vk - V II), II pk - p 11 and I~ jik - p II, where (uk,vk,pk) is 

the solutlOn obtamed for that grId, jik = iII;=l (T;/2 + T)-1/2)pk and II illS the 

discrete Ll norm per umt area. The three numbers m the next column (headed 

by "a") are max]=1.2 max(1I a;uk - a)u II, II a;vk - a] v 11), II a;pk - a]p II, and 

II a;jik - a]p II, where a) = ajax] and a; = (T)1/2 - T]-1/2)jhk In the next col-

- 18 -
TABLE ~ on-staggered N aVler-Stokes Staggered I, 

2 i N av -Stokes I I No P-a~raglDg: P-averagmg 

gndl 

I Rf = 0 15 00084 

! Per I I 00500 

I I I 00542 

I I I 

! Exact 16 I 00024 

I I ' 

: Sol i 00123 I 

I I a I 

0030 0108' .00084 .0030 0108 I 00079 .0055 1 

.0394 03921.00500 0394 0392 I 00150.00941 

.0395 .03951.00542 0395 0395 II 00513 0212 I 

.0007 0026; 00024 0007 0026 I 00018 .0013 ! , I I 

0098 0098: 00123 0098 .0098 I 00037 0025 i 

I Re = 0 i 5: 00090 0036 0113 i 00097 0031 0113 I 00080 0052 ! 

I Per I ! .00661 2086 0555\ 00978 .3452 .0682 I 00163 .0215 ! 

! I I 00562 .0447 0445

1

: 00977 0670 .0691 i .00540 0243 

! ! 00135 0100 0108! 00135 0100 0108 i 00127 0054 : 

I I I 

W(2,1) 6 '.00024 .0008 0027' .00025 .0008 .0027 00018 0013 

11-FMG .00136 0536 0146 00216 1346 .0181 00036 .0036 'I I I 

I I 00136 0119 0117 I 00209 0180 0189 00129 0057 I 

I Re = 0 : 51 00104 0041 0111 1 00097 0035 0109 I 00076 0055 

: Dlf : ! 01285 3715 .0851 I 01480 3946 0763 : 00198 .0176 I I I I I I 

j ; i .00712 0665 .05301.00971 0719.0649 .00544 0246 i 

II ! ii, I W(2,1) 16 i 00027 0011 .0028: .00026.0008.0027 I 00017 0013 ! 

I 1-FMG I 1.00337 1586 04511 00371 .1776 0264 00047 0038 I 

I 00191 .0348 .0271 .00223 .0252 .0191 00132 .0059 I 

IRe = 00 15! 00272 0433 .0253 I .00215 .0250 0097 I 00168 0180 I 

I Du i .01515 .4536 .2382\.00637 .2098 .0256 .00242 0832 I 

: .00957 .1945 .15941.00273.0168 .0147 ,.00106 0076 

I I I I I W(2,0) ,6! 00138 0547 .03571 00088.0154 .0046 ,.00039.0064 

2-FMG I .01517 .8227 .4913 I .00142 .0830 .0085 00066 .0436 I 

.01051 4062.3413 .00074.0048 .0033 .00038 .0016 I 

- 19 -
umn (headed by ae ). sImilar numbers are gIven, WIth thc long dIfference quotIent 

a; = a~ = (TJ - T

J

-] ) / (2hk ) replaung a; The next 3 columns show sImIlar sets 

of results for the case that P-averagmg IS used The remammg 2 columns gIve for 

companson results obtamed on a staggered grId wIth the same meshsIze, wIthout 

P-averagmg (Usmg pk for approximatlOn, especIally of denvatlves. shll may pay 

if 112 = 0) 

8 3 Accuracy and StabIhty Tables 1 and 2 dearly show that the exact qUasIellIptic solutIon~ (~). and the non-staggered Stokes, the latter mamly m terms 

of P} are several tImes less accurate than the correspondmg fully elhptic ones 

(~h and staggered Stokes, respectIvely), but they are stIll 0(h2) Errors m the 

hIghly oscillatmg case, exhibItmg mstablhty, could of course all be reduced to 0(1) 

(or O(h;-I) in denvabves} by enough F-averagmg (see §3) Averagmg the solutIon (row 5a, or the pk results). or taking SUitable long dIfference quotIents. cure 

the worst behaVIOr too but also some\\ hat furt her reduce the smooth-component 

accuracy, whIch nevertheless remams 0(h2) 

8 4 Poor asymptotic algebraIC convergence Denote by ). the asYmptotIC convergence factor per multIgnd cycle, 1 e . ). = (rt!rm)I/(l-m) for suffiCIently large f., 

m and f. - m, where rl IS any error (or reSIdual) norm measured at any fixed stage 

of the l-th cycle As expected (see §4 1). the usual cycles MG(k, ) yielded poor 

). for quasi-ellIptIc schemes 

In case of the Skew LaplaCIan and \' (2 1) cycles. our experIments exhIbIted 

). = .845 and ). = 96 for levels k = 4 and k = 5. respectIvely The convergence 

rate log 1/), is clearly 0(h2), as the rate of a SImple Gauss-SeIdel solver for the 

compact LaplaCIan ~h Indeed, on each subgrid (red or black) the relaxatIOn does 

look like Gauss-SeIdel for ~ hI. and the coarst' grId correctlOns are no help m case 

the black reSIduals cancel the red ones m the transfer to grId k - 1 For comparIson 

V(2, 1) cycles for the compact LaplaCian t::.,h with lexicographIC Gauss-Seidel yield 

>. ~ 12 on all grids 

SImIlarly. for the periodic Stokes problem and W(2.1) cycles, >. = 80 and 

>. = 945 were obtamed on levels 4 and 5, respectIvely, exhibItmg agam O(h2) 

rate. The rates were almost identically the same whether P averagmg was used 

or not. For companson, for staggered-grId Stokes discretIzatIons the red-black 

DGS relaxatlOn gIves>' = 30 and>' = 20 for the W(I,O) and the W(2,0) cycles. 

respectively These same excellent rates are obtamed both for the penodic and 

the DIrichlet boundary condItlOns (prOVIded some local relaxatlOn near boundaries 

IS added in the latter case) The same results are obtained for the NavIer-Stokes 

problem with small Re For large Re, dIvergence occur unless P-averaging is used 

(cf §63). 

- 20 -
85 FMG results DespIte the bad asymptotIc convergence, Tables 1 and 2 clearly 

shoV'. that results obtamed for the quasl-elhptlC cases by short FMG algorIthms 

are very good In smooth case~ they YIeld dIfferential errors practIcally as small 

as m the exact dIscrete solutIOns Moreover, m case of the unstable mode, the 

FMG results are VIsIbly much better than the exact solutIOn (preCisely because the 

bad behavIOr I~ slow to enter) In case of non-Imear equatIOns (Table 2, Re = (0) 

proper averagmg (Sec 63) IS t!vldentl) necessary for good FMG results 

86 As~'mptotlc comergence wIth new algOrIthm The _\LVG(5 ) cycle of §7 2 

ha..<; been employed to solve the ske\\ LaplaCIan problem wIth II = 3 relaxatIOn 

sweeps per cycle and WIth V(2, 1) used as the MG(4. ) mner cycle For many 

cycle~ the convergence factor per cycle wa.." steadIly between 07 and 08, or a 

convergence factor of 425 per fine-gnd relaxatIOn close to the value 447 expected 

by the smoothmg mode analYSIS (§ 7 3 J 

[1] 

[2] 

[31 , J 

[41 

[5] 

[6] 

[7J 

'81 l J 

References 

A Brandt. Stages m developmg multlgnd solutIOns Numeneal Methods for 

Engmeenng (E AbsI, R Glowmskl P LascaID.. H Veysseyre, eds ), Dunod: 

Pans, 1980. pp 23-44 

A Brandt Multlgnd solvers for non-elhptIC and smgular-perturbatIOn steadystate problems The Welzmann Institute of SCIence 1981 

A Brandt ~lultlgnd techlllques 1984 GUIde WIth applIcatIOns to flUId dynamICS Monograph AvaIlable as GMD-Studle No 85 from GMD-F1T. 

Postfach 1240, D-5205, St Augustm 1. W Germany 

A Brandt. J E Dendy. Jr and H Ruppel The multIgnd method for semlImplICit hydrodynamICS codes J Comp Phys 34 (1980). pp 348-370 

A Brandt and K Dmar, MUltigrId solutIOns to elhptlc flow problems Numer-

'Leal Methods for PDE (S V Parter. ed) AcademIC Press. 1979, pp 53-147 

J E Dendy, Jr , MultIgnd seml-lmpiIcit hydrodynamICS reVISIted Large Scale 

Sc'Lent'Lfie Computahon (S.Y Parter, ed ), AcademIC Press 1984 

V Glrault Theory of finIte difference methods on Irregular networks, SIAM 

J Num Anal 11 (1974), pp 260-282 

S Ta'asan Multlgnd Methods for HIghly OSCIllatory Problems, Ph D TheSIS 

The Weizmann InstItute of SCIence, Rehovot, Israel, 1984 

- 21 -
End of Document