Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Delaney, Allen Daniel
Master of Science - MSc
University of British Columbia
A number of experimental implementations of the multigrid algorithm for the solution of systems of partial differential equations have been produced. One program is applicable to simple nonlinear scalar equations, the others to linear equations, scalar and systems, which may be mildly stiff. All use nested grids and residual extrapolation techniques to compute solution and error estimates very economically. One version implements list based adaptive grids to further decrease both computation and storage needed for comparable problems. Each experiment was demonstrated using a set of problems with known solutions and the program performance or nonperformance discussed. Several techniques were examined to ensure that the system of difference equations representing a given problem would be convergent. The use of artificial viscosity was found to be practical in the general case, though for linear problems the use of one-sided differencing may be superior.
ON THE IMPLEMENTATION OF MULTIGRTD METHODS FOR THE NUMERICAL SOLUTION OF P A R T I A L D I F F E R E N T I A L EQUATIONS By A L L E N DANIEL DELANEY B.Sc, M c G i l l University, 1964 Ph.D., M c G i l l University, 1969 A THESIS SUBMITTED I N P A R T I A L F U L F I L L M E N T OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES The Department of Computer Science We accept t h i s t h e s i s as conforming to the required standard THE U N I V E R S I T Y OF B R I T I S H COLUMBIA November 1984 (?) A l l e n Daniel Delaney, 1984 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree a t the U n i v e r s i t y o f B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and study. I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head of my department o r by h i s o r her r e p r e s e n t a t i v e s . I t i s understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be allowed without my w r i t t e n p e r m i s s i o n . _ . . . Computer Science Department of The U n i v e r s i t y of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date I^Un^^^ /7y/fr/ DE-6 (.3/81) Multigrid Methods ( i i ) ABSTRACT A number of experimental implementations of the multigrid algorithm for the solution of systems of p a r t i a l d i f ferent ia l equations have been pro-duced . One program is applicable to simple non-linear scalar equations, the others to l inear equations, scalar and systems, which may be mildly s t i f f . A l l use nested grids and residual extrapolation techniques to compute solution and error estimates very economically. One version implements l i s t based adaptive grids to further decrease both computation and storage needed for comparable problems. Each experiment was demon-strated using a set of problems with known solu-tions and the program performance or non-performance discussed. Several techniques were examined to ensure that the system of difference equations representing a given problem would be convergent. The use of a r t i f i c i a l viscosity was found to be pract ica l in the general case, though for l inear problems the use of one-sided differencing may be superior. Multigrid Methods ( i i i ) TABLE OF CONTENTS Abstract i i Table of Contents i i i L i s t of Tables v L i s t of Figures v i Acknowledgement v i i 1. Introduction 1 1.1 The Method 1 1.2 Other implementations 3 1.3 This Implementation 4 2. Method and Implementation 8 2.1 The Algorithm 8 2.2 Grid Implementations 13 2.2.1 Minimum storage 13 2.2.2 Minimum trouble 14 2.2.3 Adaptive grids 15 3. Numerical experiments and t h e i r r e s u l t s 22 3.1 Simple nonlinear s c a l a r equations 22 3.2 Linear s c a l a r equations with 25 poss i b l e mild s t i f f n e s s 3.3 Linear systems with mild s t i f f n e s s 33 3.4 Linear systems with adaptive meshes 36 Multigrid Methods (iv) 4. Discussion 44 4.1 What has been done 44 4.2 Take home message 45 4.3 Omissions from t h i s work 46 5. References 48 6. Appendices 50 A. User i n s t r u c t i o n s 50 B. Sample source code 52 C. A sample parameter f i l e 54 D. A sample output f i l e 55 E. Source f o r the c e n t r a l routines 57 solve, relax, and taufunc Multigrid Methods (v) LIST OF TABLES I. Single nonlinear equations 24 I I . One-sided differences with FAS 28 for l i n e a r s i n g l e equations I I I . One-sided differences with 31 r e s i d u a l corrections i n relaxations only IV. One-sided differences with 32 r e s i d u a l corrections i n relaxations and coarse g r i d corrections V. A r t i f i c i a l v i s c o s i t y with FAS 34 i n both relaxations and coarse g r i d c o r r e c t i o n VI. A r t i f i c i a l v i s c o s i t y with FAS 36 for systems of equations VII. The use of the adaptive g r i d 37 implementation with a very low tolerance VIII. Actual use of adaptive grid s , 38 with a high er r o r tolerance Multigrid Methods (vi) LIST OF FIGURES 1. The M u l t i g r i d Algorithm 9 2. Finest g r i d examples 16 3. Coarsest g r i d examples 17 4. Composite of many grids 18 5. Adaptive i n t e r p o l a t i o n 20 6. Sample problems (1) to (5) 23 7. Sample problems (6) to (12 ) 26 8. Sample problems (13) to (16) 35 9. Level 8 g r i d , example (15) 40 10. Level 16 g r i d , example (15) 41 11. Level 32 g r i d , example (15) 42 12. Level 64 g r i d , example (15) 43 M u l t i g r i d Methods ( v i i ) ACKNOWLEDGEMENTS My supervisor, Dr. U r i Ascher, deserves many thanks f o r h i s patience, since I was so i n f r e -quent i n my v i s i t s to discuss t h i s p r o j e c t . Also, I much appreciated the f l e x i b i l i t y allowed me by my employer, Dr. Ryk Ward, both i n time and i n the use of computer f a c i l i t i e s . The use of the same system I normally work on was quite a boon. Last, but not l e a s t , I must thank T r i -c i a f o r not l e t t i n g me throw i n the towel when l i f e got h e c t i c and everything was taking too long. Multigrid Methods (1) 1. Introduction 1.1. The method The numerical s o l u t i o n of boundary value problems f o r p a r t i a l d i f f e r e n t i a l equations has t r a d i t i o n a l l y been performed by the d i s c r e t i z a t i o n of the domain upon which the problem i s defined and the use of i t e r a t i v e s o l u t i o n methods f or the r e s u l t i n g large sparse system of l i n e a r equations. M u l t i g r i d methods f o r the computation of solutions of such sys-tems have been discussed by many authors, e.g. [5-8,15,21], The idea behind such methods, as suggested by Brandt, [3,6] i s to con-s i d e r the d i s c r e t i z a t i o n s of the problem on a nested sequence of grids and to e x p l o i t the r e l a t i o n s h i p s among the corresponding d i s c r e t e s o l u t i o n s . Thus we can divide the s o l u t i o n process i n t o three elements: (1) removal of non-smooth err o r components from f i n e g r i d s . (2) improvement of smooth e r r o r by corrections computed on a coarser g r i d . (3) the use of nested i t e r a t i o n to provide i n i t i a l estimates f o r each f i n e r g r i d . The removal of non-smooth, i . e . high frequency, errors from f i n e grids i s normally achieved by one of the i t e r a t i v e r e l a x a t i o n tech-niques often used alone f o r these problems, such as the Gauss Se i d e l i t e r a t i o n . These relaxations have been shown to reduce high frequency error components very e f f e c t i v e l y , but low frequency components much le s s e f f i c i e n t l y . In the m u l t i g r i d process, only a few re l a x a t i o n cycles on the f i n e s t g r i d are required, and then coarse g r i d correc-tions handle the rest of the erro r f o r that g r i d . Alternate Multigrid Methods (2) r e l a x a t i o n methods, such as the preconditioning incomplete LU decom-p o s i t i o n methods, have been shown by several authors, [16-18,23,24], to be robust smoothing algorithms. For the coarse g r i d corrections, what bett e r way could be found than to r e c u r s i v e l y apply the m u l t i g r i d algorithm at the lower l e v e l . Communication between grids at d i f f e r e n t l e v e l s i s c a r r i e d out by r e s t r i c t i o n and prolongation operators. Suitable such operators are discussed by Brandt and Mol[5,19]. Eventually, of course, we reach a coarsest g r i d on which the d i s c r e t i z e d problem must be solved exactly. This g r i d i s s u f f i c i e n t l y small that the use of the r e l a x a -t i o n technique for the accurate s o l u t i o n i s an e f f i c i e n t a l t e r n a t i v e , or the s o l u t i o n f o r t h i s g r i d can be determined by d i r e c t methods such as Gaussian elimination. Except for cases where convergence i s a problem, convenience d i c t a t e s the use of the r e l a x a t i o n technique, since i t must already be a v a i l a b l e . The t h i r d element of the f u l l m u l t i g r i d algorithm, nested i t e r a -t i o n , provides an e f f i c i e n t way to get an accurate s t a r t i n g value f o r the f i n e r g r i d s . Once the di f f e r e n c e system has been approximated as well as p o s s i b l e on a given g r i d , the approximation i s c u b i c a l l y interpolated to the next f i n e r g r i d and the m u l t i g r i d process repeated on the new l e v e l . The a v a i l a b i l i t y of t h i s nested set of grids allows us to use extrapolation and defect c o r r e c t i o n methods to accelerate the convergence and improve the accuracy of the approxima-t i o n of the d i f f e r e n t i a l s o l u t i o n [2,10,14] to b e t t e r than that of the algebraic s o l u t i o n of the f i n e s t g r i d [ 3 , 5 ] . The i t e r a t i v e use of nested grids and the use of coarse g r i d corrections i n the m u l t i g r i d algorithm improve the e f f i c i e n c y of the m u l t i g r i d process [6, 20, 22] such that i t becomes an 0(N) algorithm, where N i s the number of Multigrid Methods (3) g r i d points i n the f i n e s t g r i d . This compares to 0 (N log N ) for 3/2 the most favourable cases or 0 (N ) f o r the more conventional methods. M u l t i g r i d methods have the added advantage that they are r e a d i l y amenable to adaptive refinement of grids, i . e . , during the execution of the algorithm the parts of the mesh which are s u f f i c i e n t l y accu-rate do not need to be considered on f i n e r meshes, only those areas of the g r i d which require more work are included i n the more expen-sive f i n e r g r i d s . A disadvantage of the m u l t i g r i d method i s that the implementation of the algorithm i s much more complex than the usual i t e r a t i v e s o l u -t i o n methods. I t i s t h i s disadvantage which makes t h i s t h e s i s a pro-f i t a b l e venture, that i s , to provide a background or framework upon which further experiments using the m u l t i g r i d a r c h i t e c t u r e can be t r i e d . In t h i s work I have leaned h e a v i l y on Hemker's d e s c r i p -t i o n [15] of the recursive character of the m u l t i g r i d algorithm, rather than the i t e r a t i v e d e s c r i p t i o n of Brandt[5,6]. 1.2. Other Implementations A few other groups have presented implementations of the mul-t i g r i d method. Brandt has been the "father" of m u l t i g r i d methods, publis h i n g very p r o l i f i c a l l y [ 3 - 1 0 ] . He and h i s colleagues have d i s -cussed every aspect of the process, both t h e o r e t i c a l l y and for s p e c i f i c problems, and have published a software package claimed to be e f f e c t i v e f o r general problems. This package i s a development t o o l , providing an environment and data structure f o r further mul-t i g r i d experiments. I t s u f f e r s from the fa c t that much of the algo-rithmic content of the package i s devoted to working within the Multigrid Methods (4) confines of a Fortran environment. Another work, by Dendy[ll,12], allows the user to code only the rectangular matrix problem, leaving the g r i d manipulation to a "black box". I t seems, though, that there must be a new "black box" for each c l a s s of problems, the ones dealt with by Dendy being symmetric and nonsymmetric l i n e a r s c a l a r equations. His package i s more concerned with dealing with a r b i t r a r y grids than with general problems. Stuben and Trottenberg[21] have produced a very extensive discus-sion of the theory and a p p l i c a t i o n of m u l t i g r i d methods and present a demonstration of another m u l t i g r i d package. This o f f e r i n g i s an example of the l i b r a r y approach as opposed to the "general package" approach. Foerster and Witsch[13] have offered to d i s t r i b u t e t h i s c o l l e c t i o n of programs for the s o l u t i o n of p a r t i a l d i f f e r e n t i a l equa-tions using m u l t i g r i d methods. 1.3. This Implementation The aim of t h i s work, therefore, has been to produce an implemen-t a t i o n (program) of the m u l t i g r i d method for the s o l u t i o n of p a r t i a l d i f f e r e n t i a l equations and to experiment with a v a r i e t y of aspects a r i s i n g from such attempts. The package which has been produced i s applicable to systems of equations and does ta c k l e the two major problems i n t h i s area. These are the fact that simple d i s c r e t i z a t i o n s do not always lead to systems of equations which can be solved using i t e r a t i v e techniques and that the s i z e of the systems required f o r p a r t i a l d i f f e r e n t i a l equations may tax or exceed the a v a i l a b l e compu-t a t i o n a l space and time, v i r t u a l or otherwise. The former problem requires the use of techniques such as one-sided d i f f e r e n c e s or a r t i f i c i a l v i s c o s i t y [ l ] . Both have been t r i e d Multigrid Methods (5) and t h i s work shows that the use of a r t i f i c i a l v i s c o s i t y i s a p r a c t i -c a l approach. The use of one-sided di f f e r e n c e s has been found to be useful, even superior to the v i s c o s i t y approach, f o r l i n e a r problems, but more work i s needed before any conclusions concerning nonlinear problems can be made. The l a t t e r problem, computer time and space demands, requires that some sort of adaptive mesh s e l e c t i o n be used so that parts of the domain which are "easy" can be solved with a coarse mesh, r e q u i r -ing l i t t l e storage space or computational time, and the "hard" areas of the domain solved with a p r o g r e s s i v e l y f i n e r mesh. An implementa-t i o n presented with t h i s t h e s i s uses a l i s t based data structure which allows completely adaptive meshes to be used. The choice of which language to use to best implement an algo-rithm i s a question of general i n t e r e s t . Numerical analysts have often favoured Fortran, because i t i s f a m i l i a r to most other numeri-c a l analysts, i t i s widely a v a i l a b l e , and i t s compilers are often very e f f i c i e n t . On the contrary side, Fortran's lack of data s t r u c -tures and of recursion severely r e s t r i c t the programming s t y l e and c l a r i t y of the code produced. An algorithmic, structured language l i k e A l g o l allows the production of c l e a r e r code, probably with l e s s programmer e f f o r t , p a r t i c u l a r l y f o r l i s t based and recursive algo-rithms . Pascal i s an A l g o l - l i k e language which i s beginning to approach Fortran i n popularity. With t h i s p o p u l a r i t y and with hardware design beginning to favor recursive programming, the e f f i -ciency of Pascal compilers and t h e i r code has become l e s s of a disad-vantage . M u l t i g r i d Methods (6) For these reasons, Pascal was chosen as the implementation language for t h i s p r o ject. An e a r l y implementation was a c t u a l l y coded i n the C language, but then the d e c i s i o n to switch was made. C has some advantages, p a r t i c u l a r l y i n the dynamic a l l o c a t i o n of arrays, but has disadvantages i n code c l a r i t y and i n the breadth of i t s acceptance. Although the l a t e r versions of the package are applicable only to l i n e a r systems, the adaptation to nonlinear problems should not be d i f f i c u l t . The basic algorithm was implemented i n the e a r l i e s t ver-sion, and the t r a n s f e r of these d e t a i l s to the f i n a l version should be straightforward. No attempt at optimizing performance has been made, but the l i t e r a t u r e contains many examples of high frequency smoothing schemes and low frequency c o r r e c t i o n techniques which might be t r i e d i n the attempt to optimize performance f o r s p e c i f i c d i f f i c u l t problems. In general the l i t e r a t u r e [7, 9, 17, 19] deals with the theory of p e r f o r -mance and admits that p r a c t i c a l problems require i n d i v i d u a l attention to e f f i c i e n c y f o r each case. This package i s designed as a framework upon which to b u i l d the s p e c i a l i z e d systems for r e a l problems; i t i s modular and extensible. Of course a number of classes of problems are solvable using the presented package, but the future user should be able to add h i s own classes of problems with minimal d i f f i c u l t y . As to e f f i c i e n c y , a major design goal for t h i s program was to allow f l e x i b l e use of adaptive g r i d s . This approach to e f f i c i e n c y seems much more c r u c i a l than doing a few l e s s r e l a x a t i o n sweeps on f u l l g r i d s . Multigrid Methods (7) Consideration has not been given to problems with a r b i t r a r y boun-daries, but the data structure used here could be adapted to the i r r e g u l a r edges of a domain i n a straightforward manner. The tech-nique would involve adding f i n e r and f i n e r mesh points as the boun-dary i s approached u n t i l a g r i d point i s s u f f i c i e n t l y close to the boundary, but i s s t i l l part of the organized g r i d . These g r i d points near the boundary would have to be flagged so they could be included i n a l l g r i d s . This t h e s i s i s a c t u a l l y a report on the evolution of a number of experimental codes and the experience I have had with them. The next section presents the general algorithm and discusses the implementa-t i o n d e t a i l s f o r a sequence of programs, each with a d i f f e r e n t emphasis, but each growing from the previous. A f t e r that, some r e s u l t s f o r each code are presented to demonstrate t h e i r performance, or non-performance. In the discussion I have t r i e d to summarize what I have done and o u t l i n e the many avenues which could be traveled i n a continuation of t h i s work. Multigrid Methods ( 8 ) 2. A Multigrid Method and Its Implementation 2.1. The Algorithm The numerical s o l u t i o n of a system of p a r t i a l d i f f e r e n t i a l equa-tio n s i s generally accomplished by solv i n g the system of differe n c e equations a r i s i n g from a d i s c r e t i z a t i o n of the domain in t o a set of g r i d points, i . e . , a mesh (or g r i d ) . One wishes to f i n d an "optimal" mesh where the erro r i n the computed s o l u t i o n i s s a t i s f a c t o r i l y small, and the number of g r i d points i s small enough to keep the com-pu t a t i o n a l cost reasonable. One way of approaching t h i s optimum i s to solve on very coarse grids, use the s o l u t i o n to approximate a s t a r t -ing point f o r the next f i n e r mesh, and repeat the cycle u n t i l the estimate of the erro r i s s u f f i c i e n t l y low. The m u l t i g r i d method presented here does t h i s and i t uses the a v a i l a b i l i t y of multiple g r i d solutions to improve both e f f i c i e n c y and accuracy. This m u l t i g r i d method i s outlined i n Figure 1. In t h i s algorithm //* denotes the dif f e r e n c e operator representing the system of equa-tions on a g r i d with mesh s i z e h . The l e t t e r n denotes the number of d i v i s i o n s i n one dimension of the g r i d , so h = 1/n on the unit square. We w i l l r e f e r to n as the g r i d s i z e . L represents the same operator, but working on the next coarser g r i d , with mesh s i z e H = 2h . The symbol f represents the r i g h t hand side of the sys-tem. Several parameters and d e s c r i p t i v e constants provide the l i m i t s f o r the algorithm. The constants nmin and nmax are the g r i d s i z e s f o r the coarsest and f i n e s t g r i d s . 7j and a and EF are parameters depen-dent on the order of the d i s c r e t i z a t i o n used. EF i s the extrapola-t i o n factor, which i s always unity on c o r r e c t i o n grids, but may be Multigrid Methods (9) procedure m u l t i g r i d ; l e t n = nmin i n i t i a l i z e u and f on the g r i d ; solve exactly = ; while n < nmax and an err o r estimate > an err o r tolerance do begin l e t n = 2 n ; i n t e r p o l a t e u to l e v e l n from l e v e l n/2 ; tl tl tl solve at l e v e l n the problem L u = f ; end ; end m u l t i g r i d ; procedure solve ( at l e v e l n ) ; l e t T = 0 0 n while r n > a r n / 2 do begin while convergence < 7? do rela x at l e v e l n , „, H _ .H . ih.h . H H . l e t r = r - EF I^ ( L u - Ih L u ) ; H H H solve at l e v e l n/2 the problem L u = r ; . . h h A .H H l e t u = u + I, u ; n end; end solve; Figure 1: The M u l t i g r i d Algorithm An algorithmic d e s c r i p t i o n of the m u l t i g r i d method. See the text f o r a discussion of the parameters def i n i n g and l i m i t i n g the process. greater than unity on t o p - l e v e l g r i d s . I t allows the over-correction Multigrid Methods ( I O ) of the r e s i d u a l to obtain higher order accuracy f o r smooth problems. h H The notation /,, and /. r e f e r s to r e s t r i c t i o n and prolongation hi n operators. In one part of the m u l t i g r i d process r e s i d u a l estimates computed using u values from the f i n e r of a p a i r of grids must be added to the u values of the coarser g r i d , and i n another part a c o r r e c t i o n computed on the coarser of the p a i r of grids must be added to the f i n e r g r i d . These two operators compute the u values for a g r i d at the required l e v e l from the u values at the known l e v e l . In t h i s package r e s t r i c t i o n may be e i t h e r i n j e c t i o n , i . e . , f i n e g r i d values are merely copied to the corresponding locations of the coarse g r i d , or an averaging of 9 points. Prolongation i s l i n e a r i n t e r p o l a -t i o n of points i n the coarse g r i d to produce a f i n e g r i d , which i s equivalent to 9-point prolongation. The relaxations used by t h i s implementation i n the solve routine of the algorithm are Gauss S e i d e l i t e r a t i o n s over the system of equa-tions f o r the g r i d . For l i n e a r problems the system of l i n e a r equa-tions at each g r i d point i s solved using Gaussian elimination. I f the version of the program handles nonlinear systems, the relaxations a c t u a l l y do one or more Newton i t e r a t i o n s to r e l a x the nonlinear sys-tem at each g r i d point. The heart of the m u l t i g r i d system i s the set of relaxations, and the f a c t that these relaxations, though they have very poor conver-gence properties o v e r a l l , have an excellent, i . e . f a s t , convergence rate f o r the highest frequency components of the e r r o r . When working with d i s c r e t i z a t i o n s , as we are now, i t i s appropri-ate, under c e r t a i n conditions, to represent the e r r o r as a Fourier s e r i e s , i . e . , a sum of trigonometric terms with d i f f e r i n g Multigrid Methods ( 1 1 ) wavelengths. We can divide the error, then, into three groups of components. The low frequency error components are those with wavelengths less than h/2, the high frequency components are those with wavelengths between h/2 and 2h , and the invis ib le components are those with s t i l l higher frequencies. These invis ib le errors must be avoided, rather than dealt with, as we shal l see later . The idea then, is to do a few relaxation cycles at the finest grid leve l , un t i l the convergence cr i ter ion returned by the relax routine indicates that further relaxations are of l i t t l e use, and then to correct the low frequency errors by solving on the next coarser gr id . The program parameter 7? controls the def init ion of the relaxation convergence cr i t er ion . When the rat io of the changes for two consecutive relaxations is greater than 77, high frequency conver-gence is said to be obtained. This process i s recursively continued back down to the coarsest grid leve l , which i s solved exactly. In this way usually only two or three relaxations are done at the finest level on each cycle and the rest of the work i s done on lower, cheaper levels . The extrapolation factor i s dependent on the order of the d i f f er -ence approximations used. If the grid has been solved accurately, then this extrapolation increases the accuracy of a second order difference solution to fourth order. The value used in the extrapo-lat ion of the residual is an estimate of the residual, evaluated from two grid computations as per Figure 1. This estimate, r , i s used to decide when further coarse grid corrections are useless, and, i f yes, to go on to a finer gr id . For second order differences, i t can be expected that the residual w i l l decrease by a factor of one fourth Multigrid Methods ( 1 2 ) f o r each g r i d l e v e l . The program parameter a i s given the value of t h i s expected r a t i o and i s used i n the computation of the g r i d tolerance f o r the estimated r e s i d u a l . This tolerance i s a times that obtained f o r the previous, coarser, g r i d . Another constant often used i n the m u l t i g r i d analysis, 6, i s the factor by which the error i s expected to decrease f o r each g r i d l e v e l . For a fourth order method t h i s would be one sixteenth, and the e r r o r estimate would be 6 < ° c - u f ) 1-6 My implementations are very conservative, i n that they do not use 6. Instead the c o r r e c t i o n term u -u. i s reported, thus overestimating C T the e r r o r . The reader may note that that the value of 6 changes with the order of d i f f e r e n c i n g , and that 6 has l i t t l e meaning when a r t i f i -c i a l v i s c o s i t y i s being used. Extrapolation i s only performed at the f i n e s t l e v e l . The function of the lower l e v e l solves i s to approximate the s o l u t i o n of the alge-b r a i c system of equations at the next f i n e r g r i d . The function of the extrapolation i s to improve the top l e v e l algebraic s o l u t i o n as an approximation to the s o l u t i o n of the d i f f e r e n t i a l system. The algorithm described i n Figure 1, i s known as the FAS algo-rithm [ 5 , 6 ]. A l l grids i n the algorithm are dealing with a s i m i l a r problem, an approximation of the problem being solved on the f i n e s t g r i d . The co r r e c t i o n to be applied to the f i n e r g r i d i s computed, i n the FAS case, by subtraction from the best f i n e g r i d s o l u t i o n . This contrasts with the r e s i d u a l c o r r e c t i o n approach, where the coarse grids are used to solve a problem r e l a t e d to the erro r i n the f i n e Multigrid Methods (13) g r i d , and the c o r r e c t i o n i s d i r e c t l y c a l c u l a t e d . The r e s i d u a l c o r r e c t i o n approach i s therefore applicable only to l i n e a r problems. 2.2. Grid Implementations One property of m u l t i g r i d processes i s that they perform s i m i l a r operations on grids of varying s i z e . In most of t h i s study, the larges t g r i d s i z e used was 64, and normally the minimum g r i d s i z e was 4. During the course of a computation, only one of the f i n e s t grids, with n = 64, i s required. For the coarse g r i d corrections, though, one g r i d of each s i z e n = 4,8,16, and 32 are required. Interpola-t i o n , prolongation, and r e s t r i c t i o n operations must use these d i f -ferent s i z e d g r i d s . I w i l l discuss three p o s s i b l e data structures for these g r i d s . 2.2.1. Minimum storage The obvious approach, i f we are using only uniform grids, i s to define an array f o r each g r i d which w i l l contain that g r i d and no more. A g r i d with n = 4 would, i n Pascal, be defined: array [0..4] of array [0..4] of r e a l ; and one with n = 8 as: array [0..8] of array [0..8] of r e a l ; This uses a minimum of storage, so long as f u l l grids are used. The disadvantage of t h i s system i s that the t r a n s l a t i o n of the integer indices /,/ to the (x ,y ) values they represent changes from g r i d to g r i d , and so the i n t e r p o l a t i o n and r e s t r i c t i o n functions are more complicated than i s d e s i r a b l e . Multigrid Methods ( 1 4 ) From a s t r i c t l y programming viewpoint, since Pascal was chosen as the implementation language, problems of data typing and dynamic array a l l o c a t i o n make t h i s type of data structure very inconvenient. The Pascal procedures which work with grids w i l l not allow grids of d i f f e r e n t s i z e to have the same type, and i t i s not straightforward to give arrays of d i f f e r e n t s i z e the same type. There are ways around t h i s problem, but they are u n l i k e l y to be portable, and p o r t a b i l i t y was one of the major reasons f o r choosing Pascal as the implementa-t i o n language. 2 . 2 . 2 . Minimum trouble At a considerable s a c r i f i c e i n storage, and given that we are using uniform grids , a l l the disadvantages of the minimum storage approach can be eliminated. Here a l l grids have the same type, and the r e l a t i o n between (x,y) and (/,/') i s the same regardless of the g r i d s i z e . We merely leave empty spaces f o r a l l the unused g r i d points on a given l e v e l , and a l l grids occupy a maximum amount of space. The storage cost can be e a s i l y computed. I f we l e t m =log2(nmax/nmin ) then, using the minimum storage data structure the grids would require „ m , nmax £ o . 5 2 / < - nmax i=o 2 I f each g r i d required the f u l l nmax g r i d l o c a t i o n s , the propor-t i o n a l cost of the minimum trouble g r i d r e l a t i v e to the minimum storage g r i d would be 2 m nmax _ 3m 4 2 ~ 4 — nmax Figures 2 and 3 present p i c t u r e s of the set of grids which would be Multigrid Methods ( 1 5 ) used by a computation with nm/n=4, nmax=32. Every g r i d contains a l l the (x,y) co-ordinates which are present i n the previous g r i d , thus the coarsest g r i d points are present i n a l l g r i d s . With a l l four of these grids superimposed, and with the respective g r i d points represented by "*", "+", and "@", the representation becomes Figure 4. 2.2.3. Adaptive Grids A t h i r d data structure f o r the grids does not use arrays at a l l . I t uses a l i s t oriented structure to hold the g r i d information. With the exception of the overhead for pointers, the storage cost f o r t h i s data structure i s the same as for the minimum storage array method. Each g r i d point i n Figures 1-4 i s represented by a si n g l e node of the data structure. This node contains the current value of u at that point, the indices (/,/), from which the (x,y) values can be com-puted, pointers to the nodes to the r i g h t and above, the value of the ri g h t hand side for that node, and other information which i s used i n the computation, such as erro r and r e s i d u a l estimates. The actu a l pointer to the g r i d i s the pointer to the lower l e f t corner node of the g r i d . I t would have s i m p l i f i e d programming quite a b i t i f each node had four pointers, one for each of i t s neighbors, but a decision to save those eight bytes of memory per g r i d point was made. Hindsight i n d i -cates that, at l e a s t f or the prototype implementation, i t would have been more e f f i c i e n t of my time to work out the implementation using four pointers and rewrite f o r e f f i c i e n c y l a t e r , i f at a l l . This type of data structure allows the f l e x i b i l i t y required f o r the implementation of adaptive g r i d algorithms. For problems which Multigrid Methods ( 1 6 ) Multigrid Methods ( 1 7 ) The n=8 g r i d The n=4 g r i d * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i x-> x-> Figure 3 : Coarsest g r i d examples A diagram representing the g r i d points used at the two coarsest lev-e l s when nmln =4 Multigrid Methods ( 1 8 ) @_ *—|— *—@— *—1-—*—@—*—|— *_@_ *—|— *_@ H—*—|—*—H—*—|—*—|—*—|—*—|—*—|—*—(. @—*—|—*—@—*—|—*_@—*—|—*—@_*—|—*—@ -1—*—|—*—|—*—|—*—|—*—|—*—1_*—|—*—(. @—*—|— *—@— *—|— *—@— *—|— *—@— *—|— *—@ H—*—|—*—|—*—|—*—|—*—|—*—|—*—|—*—(. * *_*_*_*_*_*_*_*_*_*_*_*_*_*_*_*_* l i y @—*—|—*—(§—*—|—*—@—*—|—*_@_*—|—*_@ H—*—1~*—|—*—(-—*—|—*—|—*—|—*—|—*—(-* _ * _ * _ * _ * _ * _ * _ * • _ * _ * _ * _ * • _ * • _ * _ * _ * _ * @ - * — ( - - * — ( § - *—|— *-@— * — 1 ~ *—@— * — 1 ~ *—@ x - > Figure 4: Composite of many grids A representation of a g r i d with s i z e n=32. A l l points are members of the f i n e s t g r i d . Only the "*" points are i n the n =16 g r i d , the "+" points are the n=8 g r i d , and the "@" points are the coarsest g r i d , at n=4. require very f i n e meshes i n selected areas of the g r i d , i t i s essen-t i a l that the fi n e meshes need not encompass the en t i r e domain. A Multigrid Methods (19) l i s t oriented data structure allows us to a l l o c a t e only those regions of the grids f o r which err o r estimates i n d i c a t e that more refinement i s necessary. In the f i n a l version of the program presented with t h i s t h e s i s , the e r r o r estimates at each g r i d point are used by the i n t e r p o l a t i o n procedure to decide which new points are needed at the next l e v e l . Figure 5 shows a hypothetical p a i r of grids, demonstrating how the adaptive i n t e r p o l a t i o n proceeds. On the n =8 g r i d the "<§" g r i d points are those with large e r r o r estimates. These same g r i d points are represented by "@" symbols again i n the n=16 g r i d . On t h i s f i n e r g r i d , there are two types of new points, i n a d d i t i o n to a l l the coarse g r i d points. The points represented by "x" characters are those which w i l l be relaxed on the f i n e g r i d . The other new points are there merely to allow the relaxations to proceed on the "x" points. I f any point has an e r r o r estimate l a r g e r than the tolerance for the g r i d , a l l points within the square defined by the 9 point group around the g r i d point i n question w i l l be placed i n the r e s u l t i n g g r i d . This kind of arrangement w i l l lead to some g r i d points which do not have the neighbors required for r e l a x a t i o n . Such g r i d points, however, w i l l always e i t h e r have s a t i s f a c t o r y e r r o r estimates, or be the cubic i n t e r p o l a t e of two such points. The r e l a x a t i o n sweeps cover the e n t i r e g r i d each time, but the computations are done only on the "x" and "@" g r i d p oints. The reader may notice that there are seemingly unnecessary points on the boundaries. These are not used f o r the computation, but are needed f o r t r a v e r s i n g the g r i d , ensuring o r d e r l y access to a l l the Multigrid Methods (20) The n=8 g r i d The n=16 g r i d * * * * * * * * * * * * * * * X X X * * x @ x * * * * * * * X X X * * * * * * * * * * * * * * * * * * * * * * * * * x x x x x x x * * @ * * * * * * * X X X @ X X X * * * X X X X X X X * I I @ <a @ * * y * * * * x @ x @ x @ x * * * x x x x x x x * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * x-> x-> Figure 5: Adaptive i n t e r p o l a t i o n A representation of a p a i r of grids, demonstrating the adaptive i n -te r p o l a t i o n procedure used to produce the next f i n e r g r i d . g r i d points. Multigrid Methods ( 2 1 ) During the course of t h i s work, mented. The f i r s t three of these method, the f i n a l version used the ture. four program versions were imple-used the Minimum Trouble storage adaptive, l i s t based, data s t r u c -Multigrid Methods (22) 3. Numerical experiments and their results As mentioned before, several d i f f e r e n t m u l t i g r i d implementations were produced, each growing from the previous one. Here i s a d e s c r i p t i o n of the programs I have b u i l t , the problems incurred, the s t r a t e g i e s t r i e d , and the r e s u l t s obtained. 3.1. Simple nonlinear scalar equations This i n i t i a l version of the package attacked the problem: a Uxx + c uyy + d ux + e uy + g = 0 where a , c , d, e , and g are s c a l a r functions of x , y, and u . Since g can be an a r b i t r a r y function of x, y, u, i t was deemed unnecessary to include a r i g h t hand side i n the implementation. The user was required to provide the functions f o r g(x,y,u ) and i t s par-t i a l d e r i v a t i v e with respect to u , as w e l l as a function s p e c i f y i n g the boundary values of the s o l u t i o n . Figure 6 l i s t s some example problems of t h i s prototype which were used to demonstrate t h i s f i r s t implementation. The motivation f o r presenting these p a r t i c u l a r r e s u l t s i s to demonstrate that at l e a s t one version of these implementations was e f f e c t i v e on nonlinear problems. Table I summarizes the r e s u l t s of t h i s study. I t can be seen that the convergence of the process i s close to fourth order for most cases. Example 2 i s an exception, and t h i s i s because the exact s o l u t i o n i s a low order polynomial, hence the coarse g r i d solutions are close to the exact one not j u s t because of small h . In a l l cases the convergence rate dropped o f f at a g r i d -s i z e of 64. This behaviour was caused by the parameter settings i n the program being tuned f o r e f f i c i e n c y , so a minimal number of r e l a x -Multigrid Methods ( 2 3 ) V + uyy ~ u* = - s i n 2 ( x ) e 2 y ( D y Solution: u = sm(x) e u t u -eu = 4 - e*2^2 <2> xx yy 2 2 Solution: u = x + y 2 2 2 u + u - u + 2 u = - s i n (x) s i n ( y ) (3) xx yy \ / \ / y Solution: u = sin (x ) sin(y) 3 3 3 u + u - u + 2 u = - s i n ( x ) s i n ( y ) (4) xx yy v / \/ / u ± „ (sin(x)sin(y)) / C N u u + u . - e + 2 u = - e v v / v / " (5) Solution: u = s/n(x) sin (y ) j uxx yy Solution: u = sin(x) sin(y) Figure 6: Sample problems (1) to (5) Example problems 1 through 5. These were used i n the demonstration of the f i r s t implementation, which could handle nonlinear problems but only f o r equations producing a diagonally dominant d i s c r e t i z a t i o n ma-t r i x . a tion sweeps were done on the f i n e s t g r i d s . A readjustment of the parameters to allow more m u l t i g r i d cycles would improve the conver-gence at the cost of a d d i t i o n a l computation. An attempt at adaptive grids was made with t h i s version, and i t was very educational. The idea was to have two l i m i t s , SMALL and BIG. A l l m u l t i g r i d operations were performed on grids of s i z e SMALL or le s s . When i t became necessary to go from the SMALL g r i d to that of siz e 2*SMALL the g r i d was p a r t i t i o n e d i n t o four grids of si z e SMALL and each g r i d was treated separately. Once the four grids were solved, a region covering the boundaries between the grids was relaxed a few times. I f the er r o r estimate on a given g r i d was Multigrid Methods (24) e r r o r Problem g r i d e r r o r r a t i o 1 8 3.0e-5 16.6 16 2. 4e-6 12.8 32 1.8e-7 13.1 64 2.7e-8 7.0 16 6.4e-8 5.2 32 5.1e-9 12.7 64 l . l e - 9 4.6 8 16 32 64 1.le-5 8.0e-7 6.4e-8 1.4e-8 16.0 13 . 8 12.5 4.6 8 16 32 64 1.le-5 8.0e-7 6.4e-8 1.4e-8 16.0 13.8 12.5 4.6 8 16 32 64 1.le-5 8.0e-7 6.4e-8 1.4e-8 16 .0 13 . 8 12.5 4.6 Table I: Single nonlinear equations A summary of the errors obtained and the e r r o r r a t i o between adjacent grids f o r several problems, using the o r i g i n a l implementation which was applicable to sing l e nonlinear equations. s a t i s f a c t o r y , i t was not necessary that i t be solved. I f the e r r o r was concentrated i n one or more of the corners, t h i s approach might have been f e a s i b l e , but i t was found that e r r o r on the Multigrid Methods (25) boundaries between subgrids did not improve with decreasing mesh s i z e . In retrospect t h i s could have been expected, since the over-lapping relaxations can remove only the high frequency er r o r from the subgrid boundaries, and there i s no coarse g r i d c o r r e c t i o n to remove the low frequency e r r o r s . This approach was abandoned f o r a number of reasons: I t i s not applicable to systems of equations, the use of adaptive meshes was not successful, and no p r o v i s i o n i s made f o r problems which y i e l d a system of equations which i s not diagonally dominant, which leaves out many i n t e r e s t i n g problems. Also, the fineness of the mesh i s l i m i t e d due to the very i n e f f i c i e n t use of storage space f o r the grid s . This storage method was chosen f o r ease of programming, and i t severely l i m i t s p r a c t i c a l g r i d s i z e . There were some accomplishments from t h i s part of the project, besides experience. I t was demonstrated that the use of FAS and a newton i t e r a t i o n i n r e l a x a t i o n i s a f e a s i b l e way to approach non-l i n e a r problems. Further, the accuracy obtained using the r extrapo-l a t i o n method was demonstrated to be fourth order f o r the smooth t e s t problems. One approach to adaptive grids was eliminated from con-s i d e r a t i o n . 3.2. Linear scalar equations with mild stiffness Proceeding i n a stepwise fashion, the next program implementation addressed the question of robustness. Several d i f f e r e n t approaches were used f o r problems which, with the usual centered d i f f e r e n c i n g , did not y i e l d a system of equations which had a diagonally dominant matrix. To s i m p l i f y programming, the same storage i n e f f i c i e n t data structure was used, and problems were l i m i t e d to a si n g l e l i n e a r Multigrid Methods (26) equation. Thus the c l a s s of equations was 9 "xx + C V + d Ux + 6 Uy + ' U = 9 Again the user was required to provide the functions d e f i n i n g the problem, but the function g(x,y,u) described i n the previous section was s p l i t i n t o g{x,y,u) = f (x ,y) u - rhs(x,y). Since the problems must be l i n e a r , no d e r i v a t i v e s were necessary. The problems approached with t h i s package are shown i n Figure 7. Examples (8) to (12) have the same representation, using ft value of 1, 8, 16, 32, and 64. To avoid overflow during computation ft Q the s o l u t i o n i n the l a s t two cases was scaled by d i v i d i n g by e With the values of mesh s i z e h used here, the l a t t e r three problems yielded systems of equations which were not diagonally dominant on one or more of the coarsest grids when the usual centered differences were used. u + u - 2(7 = 0 (6) . , , xx yy Solution: u = e U + U = 0 (7) xx yy Solution: u = sin(x) e 7 r + i r " ux + uy = 0 <8-12> e e 7 xft -yft Solution: u = (1-e ) ( l _ e )• Figure 7: Sample problems (6) to (12) The sample problems used to t e s t the implementations which could handle s i n g l e l i n e a r equations. Examples (8) to (12) use ft values of 1. 8. 16. 32. and 64. e Multigrid Methods (27) Several methods were u t i l i z e d i n order to solve t h i s convergence problem. Two new r e l a x a t i o n routines were implemented. One group of methods u t i l i z e d one-sided di f f e r e n c e s at g r i d points which do not s a t i s f y the diagonal dominance c r i t e r i o n . Another method u t i l i z e d only second order d i f f e r e n c e s , but used a r t i f i c i a l v i s c o s i t y on the g r i d points which required i t . In the constant c o e f f i c i e n t problems examined here, where any part of a g r i d requires a r t i f i c i a l v i s c o s -i t y , the e n t i r e g r i d would require i t . As a side issue i n t h i s study, two types of c o r r e c t i o n g r i d solve routines were examined. The o r i g i n a l such routine i s d i f f e r e n t from the top l e v e l solve i n that relaxations are done both before and a f t e r the coarse g r i d c o r r e c t i o n step. The other approach was to use i d e n t i c a l solve routines at both stages, but avoiding the r extrapo-l a t i o n i n the c o r r e c t i o n routines. The l a t t e r involves fewer relaxa-t i o n cycles i n t o t a l , i s l e s s complicated i n programming, and allows more c o n t r o l of the c o r r e c t i o n c y c l e s . My findings were that the second approach was equivalent, i f not s l i g h t l y better, i n accuracy, than the f i r s t , so i t was adopted for general use. Using f i r s t order differences i n the relaxations created some s i g n i f i c a n t problems. When, during the course of m u l t i g r i d execu-t i o n , i t became necessary to solve on a f i n e r g r i d , t h i s being the f i r s t g r i d which i s f i n e enough to use centered d i f f e r e n c i n g , the jump from a f i r s t order g r i d to a second order g r i d did not proceed smoothly. I t was desirable to use the FAS algorithm, since i t i s needed f o r extension of the method to nonlinear problems but, with FAS, the problem solved on each c o r r e c t i o n g r i d i s a close approxima-t i o n to the o r i g i n a l problem and the e r r o r on a f i r s t order correc-t i o n g r i d was found to be too large to be useful as a s t a r t i n g point Multigrid Methods (28) f o r the next g r i d . For l i n e a r problems, the use of the corrections to solve only the r e s i d u a l problem allowed f i r s t order grids to be u s e f u l . In order to preserve the fourth order convergence properties of the process, i t was necessary to use an extrapolated r e s i d u a l value on the f i n e s t g r i d . example n cgerr[n] err[n/2] r a t i o 6) 8 1.6e-5 4 .9e-4 31 16 1.le-6 1 .7e-5 15 7) 8 1.le-5 3 .2e-4 29 16 8.6e-7 1 .2e-5 14 8) 8 5.0e-6 8 . 3e-5 17 R =1 e 9) 8 3.8e-2 8 .2e-l 22 R =8 e 16 4.8e-3 1 .oe-l 21 32 5. 0e-4 5 .5e-3 11 64 5.3e-5 5 .Oe-4 9.4 10) 8 16 2.8e-l 2.4e-l 1 2 .3e-l . 8e-l 0.5 1.1 R =16 e 32 2.4e-l 2 . 4e-l 1.0 64 2.4e-l 2 . 4e-l 1.0 Table I I : One-sided d i f f e r e n c e s with FAS fo r l i n e a r s i n g l e equations Results f o r some of the sample problems from the implementation f o r l i n e a r s i n g l e equations which used one sided d i f f e r e n c i n g when neces-sary, using the FAS co r r e c t i o n scheme and using e i t h e r centered or one sided di f f e r e n c e s when computing T values. Presented i n the table are the example number, the g r i d l e v e l n , the maximum erro r on the coarse g r i d at l e v e l n, the maximum erro r on the en t i r e previous g r i d , and t h e i r r a t i o . At the l e f t are some d e s c r i p t i v e comments about some of the examples. Multigrid Methods ( 2 9 ) Tables II to VI present summaries of the r e s u l t s from some of the combinations t r i e d . In each of the tables are presented the erro r found at a given g r i d l e v e l , the erro r found at the previous g r i d l e v e l , and t h e i r r a t i o . Comparisons were made only at g r i d points which occurred i n both grids, since the newest gridpoints on the n-l e v e l g r i d have no counterparts f o r comparison. A r a t i o of sixteen would i n d i c a t e fourth order convergence, four would indi c a t e second order convergence, and two f i r s t order. The r e s u l t s i n Table II are from a program which used f i r s t order r e l a x a t i o n when necessary, but used centered d i f f e r e n c i n g e x c l u s i v e l y when computing the value of r . This being the case, extrapolation of 1.333 was used at the top l e v e l i n a l l cases. Results were not too encouraging, i n d i c a t i n g that coarse g r i d errors i n c o r r e c t i o n cycles were propagated up to the f i n e g r i d l e v e l s . Results f o r the easy problems, those which do not require the use of f i r s t order d i f -ferencing, were good, the others were not. The problems (11) and (12), with flg of 32 and 64, were both very s i m i l a r to (10), with s l i g h t l y higher e r r o r s . I d e n t i c a l r e s u l t s were obtained from the implementation where f i r s t order differences were used when necessary, and the computation of r used f i r s t order differences when the f i n e g r i d required i t . Extrapolation of 1.333 was used f o r g r i d points which are f u l l y second order, and a value of 2.0 was used f o r g r i d points which were f i r s t order, at the top l e v e l only. Again the easy problems were fi n e , the more d i f f i c u l t problems had erro r which was not removed on fi n e g r i d s , and the r e s u l t s were e s s e n t i a l l y i d e n t i c a l to those presented above. Multigrid Methods (30) A t h i r d implementation used 1st order d i f f e r e n c i n g only i n the relaxations and only when necessary to obtain diagonal dominance, but the coarse g r i d c o r r e c t i o n was a r e s i d u a l c o r r e c t i o n only. This i s not the FAS algorithm, but corresponds more c l o s e l y to Brandt's cycle C algorithm [5,10]. At the top l e v e l the r e s i d u a l computed i s an extrapolated one, which improves the order of the top l e v e l accuracy from second to fourth order. The r e s u l t s , presented i n Table I I I , are somewhat more encouraging. Example problem (10) requires f i r s t order d i f f e r e n c i n g only on the g r i d at n =4, (11) at n = 4 and n = 8 and (12) when n = 4,8, and 16. Results are much improved over the f i r s t attempts. The lowest g r i d i s solved exactly, and t h i s seems to y i e l d a b e t t e r er r o r than expected. A f t e r t h i s g r i d , we seem to be ge t t i n g approximately f i r s t order convergence u n t i l both grids involved i n the t r a n s f e r are second order. I t i s unfortunate that t h i s r e s i d u a l c o r r e c t i o n i s applicable only to l i n e a r problems. The behaviour noted here i s reasonable, considering that the er r o r generated on the coarse grids i s 0(h) i n the r e s i d u a l , which i t s e l f i s 0(h), so the t o t a l c o r r e c t i o n should 2 be 0(h ). In the FAS case, each c o r r e c t i o n cycle solves a problem very close to the o r i g i n a l , and a subtraction i s performed to get the c o r r e c t i o n f a c t o r . I f an 0(h) e r r o r i s generated on a coarse g r i d , t h i s e n t i r e e r r o r would be brought up with the c o r r e c t i o n . In Table IV r e s u l t summaries are presented f o r a method which uses f i r s t order d i f f e r e n c i n g when necessary, i n both the relaxations and i n the computation of the r e s i d u a l . Comments about a p p l i c a b i l i t y to l i n e a r problems only are relevant here as w e l l . Multigrid Methods ( 3 1 ) example n cgerr[n] err[n/2] r a t i o 6) 8 1.6e-5 4. 9e-4 31 16 1.le-6 1. 7e-5 15 32 7.3e-8 1. le-6 15 7) 8 l , l e - 5 3 . 2e-4 29 16 8.5e-7 1. 2e-5 14 32 5,7e-8 8. 6e-7 15 8) 8 5.Oe-6 8. 3e-5 17 R =1 16 4.7e-7 6. 2e-6 13 e 9) 8 3.8e-2 8. 2e-l 22 R =8 16 4.8e-3 1. Oe-1 21 e 32 5.2e-4 5. 5e-3 11 64 5.5e-5 5. 2e-4 9.4 10) 8 3.2e-2 1. 3e-l 4.0 R =16 16 4.9e-3 8 . le-2 17 e 32 3.6e-4 1. 0e-2 28 64 4.5e-5 9. 9e-4 22 11) 8 1.5e-l 9. Oe-2 0.6 R =32 16 2.7e-2 1. 5e-l 5.6 e 32 5.Oe-3 8. 5e-2 17 64 3.6e-4 1. 0e-2 28 12) 8 4.0e-2 5 . 0e-2 1.2 R =64 16 2.0e-2 1. 3e-l 6.5 e 32 2.9e-2 1. 6e-l 5.5 64 5.Oe-3 8. 7e-2 17 Table I I I : One-sided dif f e r e n c e s with r e s i d u a l corrections only i n the relaxations Results f o r some of the sample problems from the implementation using f i r s t order differences when necessary, using coarse g r i d corrections to estimate res i d u a l s only. Multigrid Methods (32) example n 10) 8 16 32 64 cgerr[n] 3.2e-2 4.9e-3 3.6e-4 4.5e-5 err[n/2] l . 3 e - l 8.1e-2 1.0e-2 9.9e-4 r a t i o 4.0 17 28 22 R =16 e 11) 8 16 32 64 1.9e-2 2.7e-2 5.0e-3 3.6e-4 9.0e-2 1.3e-l 8.5e-2 1.0e-2 4.7 4.8 17 28 R =32 e 12) 8 16 32 64 2.4e-2 1.7e-2 2.7e-2 5.0e-3 5.0e-2 8.le-2 1.4e-l 8.5e-2 2.1 4.8 5.2 17 R =64 e Table IV: One-sided di f f e r e n c e s with r e s i d u a l corrections i n relaxations and i n coarse g r i d corrections. Results f o r some of the sample problems from the implementation using f i r s t order differences when necessary, i n both the relaxations and the coarse g r i d corrections. Here again the co r r e c t i o n was a r e s i d u a l c o r r e c t i o n only. Problems (6) to (9) were i d e n t i c a l to the above, since no f i r s t order d i f f e r e n c i n g i s used, and the r e s u l t s f o r the other three prob-lems are i n Table IV. For problems (10 ),(11),and (12 ) the err o r was smoothed out over the f i r s t order grid s , but the end r e s u l t when the second order grids were reached was the same. This may imply that the extra e f f o r t of doing f i r s t order coarse g r i d corrections may not be worth the t r o u -b l e . Enough experimentation has not been done to indic a t e whether t h i s approach would be preferable f o r even higher values f o r R . As f a r as the number of relaxations was concerned, t h i s scheme was observed to be s l i g h t l y more expensive than the corresponding scheme without the matching coarse g r i d corrections, and i t would Multigrid Methods ( 3 3 ) have more overhead, i n the t e s t i n g and recomputing at f i r s t order g r i d p oints. Another implementation makes use of a r t i f i c i a l v i s c o s i t y to han-dle the diagonal dominance requirement. Table V summarizes the r e s u l t s f o r the problems of i n t e r e s t . Again r e s u l t s f o r examples (6) to (9) were i d e n t i c a l to those f o r the rest of the schemes. For examples (10) to (12) the r e s u l t s i n Table V show that the er r o r i s not improving at a l l on the grids which require the use of a r t i f i c i a l v i s c o s i t y , but convergence bounces back very w e l l once that con-s t r a i n t i s gone. This i s a c t u a l l y reasonable, since at those g r i d points d i f f e r e n t problems are being solved on the two g r i d s . For completeness, another experiment was done where relaxations used a r t i f i c i a l v i s c o s i t y , as above, but the r computation and coarse g r i d c o r r e c t i o n did not take t h i s i n t o account. As usual, r e s u l t s f o r examples (6) to (9) were i d e n t i c a l to previous r e s u l t s , but for the other examples, r e s u l t s were t e r r i b l e , the coarse g r i d c o r r e c t i o n obviously s p o i l i n g things. This contrasts with the r e s i d u a l correc-t i o n methods discussed above, where the matching of the coarse g r i d c o r r e c t i o n with the relaxations seemed to be of l i t t l e importance. 3 . 3 . Linear systems with mild stiffness A f t e r completing the experiments outlined i n the previous sec-t i o n , i t was time to modify the package to operate on systems of d i f -f e r e n t i a l equations. The c l a s s of problems handled by t h i s version i s : a u + c u + d u + e u + f u = g xx yy x y a Here, however, the c o e f f i c i e n t s are matrices, the r i g h t hand side and Multigrid Methods (34) example n cgerrfn] err[n/2] r a t i o 10) 8 1.8-2 1.5e-2 0.83 R =16 16 5.0-3 l . l e - 1 22 e 32 3.6e-4 1.0e-2 28 64 4.5e-5 1.0e-3 22 11) 8 16 32 64 3.4e-4 1.8e-2 5.0e-3 3.6e-4 2 1. 1. 1. 9e-4 8e-2 2e-l 0e-2 0.85 1.0 25 28 R =32 e 12) 8 16 32 64 1, 3. 1. 5, le-7 4e-4 8e-2 le-3 9, 3. 1. 1. 8e-8 3e-4 8e-2 2e-l 0.89 0.97 1.0 23.5 R =64 e Table V: A r t i f i c i a l v i s c o s i t y with FAS i n both relaxations and coarse g r i d c o r r e c t i o n Results f o r some of the sample problems from the implementation using a r t i f i c i a l v i s c o s i t y to ensure that the l i n e a r system representing the problem i s diagonally dominant. Centered d i f f e r e n c e s are used i n both relaxations and i n the computation of T f o r the coarse g r i d c o r r e c t i o n . s o l u t i o n are vectors. The example problems upon which t h i s version was tested are presented i n Figure 8. Problem sample (14) used R = 1 , sample (15) R =16, and sample (16 ) R = 32. 6 G © The only version of t h i s program which has been implemented i s one using a r t i f i c i a l v i s c o s i t y when a g r i d point y i e l d s a system of equations which does not meet the diagonal dominance c r i t e r i o n . At each g r i d point the diagonal dominance t e s t used was Multigrid Methods (35) U, + U, = 0 (13) lxx lyy u + u -u - 2u = o 2xx 2yy lx 2 Solution: u = 4 cos(x) sinh(y) u2 = 2x s/n (x ) s/n/7 (y ) U, + U, - 2U, = 0 (14-16) lxx lyy l 2xx ^ 2yy R R 2x 2y e e 7 Solution: = e xR -yR u2 = ( l - e e ) ( l - e e ) Figure 8: Sample problems (13 ) to (16 ) The sample problems used to t e s t the implementations which could handle systems of l i n e a r equations. Examples (14) to (16 ) used R values of 1, 16, and 32. e for each / . Table VI presents the errors and convergences. Examples (13) and (14) require no a r t i f i c i a l v i s c o s i t y and r e s u l t s are excellent. Exam-ples (15) and (16) give r e s u l t s very s i m i l a r to the si n g l e equation examples discussed above. The erro r i s not improved on the grids r e q u i r i n g a r t i f i c i a l v i s c o s i t y , but the s t a b i l i t y of the grids i s maintained u n t i l s u f f i c i e n t l y f i n e grids are reached. Multigrid Methods (36) example n cgerrfn] err[n/2] r a t i o 13) 8 1.6e-5 4.7e-4 29.4 16 l . l e - 6 1.6e-5 14.5 32 6.9e-8 l . l e - 6 15.9 16 4.4e-9 6.9e-8 15.7 14) 8 16 32 64 1.6e-5 1.le-6 7.4e-8 4.8e-9 4.9e-4 1.7e-5 1.le-6 7.4e-8 30.6 15.5 14.9 15.7 Ft =1 e 15) 8 16 32 64 1.8e-2 5,0e-3 3.6e-4 4.5e-5 1.5e-2 1.2e-l 1.0e-2 1.0e-3 0.83 24 28 22 R =16 e 16) 8 16 32 64 4e-4 8e-2 le-3 6e-4 4. 1. 1, 1, 9e-4 8e-2 2e-l 0e-2 1.4 1.0 23.5 28 R =32 e Table VI: A r t i f i c i a l v i s c o s i t y and FAS f o r systems of equations Results are presented f o r some sample problems from the implementa-t i o n which w i l l handle systems of equations, and uses a r t i f i c i a l v i s c o s i t y . This implementation i s not adaptive and uses the i n e f f i -c i e n t storage method discussed i n the te x t . 3.4. Linear systems, with adaptive meshes One of the conclusions from the r e s u l t s so f a r obtained i s that grids f i n e r than those used above w i l l be necessary f o r many prob-lems . For t h i s i t i s necessary that adaptive grids be implemented as part of an adaptive program which w i l l check e r r o r estimates on a g r i d point basis, rather than on a whole g r i d b a s i s , and make the g r i d f i n e r only i n the regions which require t h i s . The l a t e s t ver-sion of t h i s m u l t i g r i d package has the a b i l i t y to adaptively u t i l i z e f i n e r meshes only i n regions where i t i s necessary. Multigrid Methods (37) The same set of problems was treated as i n the previous section, with a tolerance s u f f i c i e n t l y low that f u l l grids were used at a l l l e v e l s . The r e s u l t s are tabulated i n Table VII. Problem (13) i s an example of a problem which has a non-diagonal term i n one of the c o e f f i c i e n t matrices. This example i s smooth and t h i s i s r e f l e c t e d i n the exc e l l e n t convergence r e s u l t . Results f o r problem (14) are approximately as expected, f o r (15) and (16) there seems to be a "catch up" g r i d at the second g r i d which does not use a r t i f i c i a l v i s c o s i t y . As long as the problem requires example n 13) 8 16 32 64 cgerr[n] 9.le-6 8.le-7 5.4e-8 3,5e-9 err[n/2] 4.7e-4 1.6e-5 1.le-6 6.8e-8 r a t i o 51.6 19. 8 20.4 19.4 14) 8 16 32 64 7e-5 le-6 9e-8 8e-9 4. 1. ,9e-4 ,0e-5 9.6e-7 5.9e-8 28 . 8 9.1 16.3 15.5 Ft =1 e 15) 8 16 32 64 1, 4. 1. 8e-2 7e-3 Oe-4 4.5e-5 1, 1. 1. 9. 5e-2 2e-l 0e-2 9e-4 0. 83 25.5 100 22.0 R =16 e 16) 8 16 32 64 3.4e-4 1.8e-2 4.7e-3 l . l e - 4 4.9e-4 1.8e-2 1.2e-l 1.0e-2 1.4 1.0 25.5 90.9 R =32 e Table VII: The use of the adaptive g r i d implementation with a very low tolerance Results from the adaptive implementation, but with the tolerance set low enough that no adaptive g r i d formations were used. Multigrid Methods (38) a r t i f i c i a l v i s c o s i t y , we cannot r e a l l y expect an improvement i n the maximum error. What happens i s that we get c l o s e r and c l o s e r to the s i n g u l a r i t y and keep the e r r o r there reasonably small. When the same problems are attacked with a low tolerance l i m i t , with the idea of s a c r i f i c i n g high accuracy to achieve b e t t e r computa-t i o n time, the accuracy r e s u l t s f o r example (14) are as i n Table VIII. Runs f o r the other examples terminated at g r i d s i z e s of 8 or 16 because t h e i r e r r o r estimates were s u f f i c i e n t l y small. As can be seen, once we reach the tolerance l e v e l , we see no improvement i n the actual e r r o r . As a matter of f a c t , there i s some backward movement i n Problem (15) : R =16 n cgerr[n] err[n/2] r a t i o i n t e r n a l points not solved proportional g r i d s i z e 8 1.8e-2 1.5e-2 1.2 49 43 43/49 = 0, .88 16 4.8e-3 1.2e-l 25 .0 279 130 130/256 = 0. .51 32 l . l e - 4 1.Oe-2 90.9 789 368 368/1024 = 0, .36 64 6.6e-3 9.9e-4 0.16 2158 1169 1169/4096 = 0, ,29 128 1.8e-2 1.8e-2 1.0 6402 4887 4887/16384 = 0. 30 Table VIII: Actual use of adaptive grids(high e r r o r tolerance) Results from the adaptive implementation, using a tolerance low enough that incomplete grids were used i n some cases. Presented here are the g r i d s i z e n, the maximum er r o r found on the coarse g r i d at l e v e l n , the maximum er r o r found on the g r i d at l e v e l n/2, which corresponds to the coarse g r i d at l e v e l n , the r a t i o of these two er-rors, the t o t a l number of g r i d points which e x i s t i n t h i s g r i d , the number of these g r i d points which are accurate enough that they w i l l not be interpolated to the f i n e r g r i d , and the proportional g r i d s i z e r e l a t i v e to a f u l l g r i d . Multigrid Methods ( 3 9 ) the e r r o r obtained. This may be due to some leverage of e r r o r from points considered to be solved s a t i s f a c t o r i l y . This behaviour i s the s a c r i f i c e made for the computational savings. The l a s t three columns show the number of g r i d points e x i s t i n g i n the grids at the various l e v e l s , the number of g r i d points being relaxed, and the r a t i o of the l a t t e r to the maximum number of points i n the g r i d at that point. This l a s t value i s approximately the proportional work being done at that g r i d l e v e l r e l a t i v e to the non-adaptive method. Since each suc-cessive g r i d i s four times the s i z e of the previous, the r e l a t i v e values on the f i n e s t g r i d are the most important. In Figures 9 to 12 we give g r i d diagrams representing the g r i d points at the four l e v e l s which were a c t u a l l y used i n the adaptive s o l u t i o n of example (14). The a s t e r i s k s represent the g r i d points where a s o l u t i o n was found to an estimated accuracy b e t t e r than the tolerance and the "@" characters represent those which were not. These l a t t e r points are interpolated to provide the f i n e r g r i d points for the next l e v e l . As can be seen, only about a quarter of the n = 64 g r i d needed i n t e r p o l a t i o n , which has allowed the program to proceed to even f i n e r grids for the more d i f f i c u l t areas of the domain. Multigrid Methods (40) * * * * * * * * * * * @ @ @ @ @ @ * * @ @ @ @ @ @ @ * * @ @ @ @ @ @ @ * @ @ @ @ @ @ < a * * @ @ @ @ @ @ @ * * * @ @ @ @ @ @ * * * * * * @ @ @ * * * * * * * * * * Figure 9: Level 8 g r i d , example (15) The g r i d at l e v e l n=8 used on example (15). This i s a f u l l g r i d , and the e r r o r estimates on the a s t e r i s k points are les s than the t o l e r -ance, while those on the "@" points are greater than the tolerance. Multigrid Methods ( 4 1 ) * * * * * * * * * * * * * * * * * * * * * * * * @ @ < a @ ( a @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ ( a @ ( a * * * * * * * * @ @ ( a @ ( a @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * @ ( a @ @ @ ( a ( a @ @ * * * * * * * * @ @ @ ( a @ ( a ( a @ ( a * * * * * * * @ @ @ @ @ @ @ @ @ * * * * * * * * * @ @ @ @ @ @ @ * * * @ @ @ @ @ @ @ * * * * * * * * * * * * * * * * * * Figure 10: Level 16 gr i d , example (15) The g r i d at l e v e l n =16 used on example (15). This i s not a f u l l g r i d , and the erro r estimates on the a s t e r i s k points are le s s than the tolerance, while those on the "@" points are greater than the tolerance. The blank areas of the g r i d have no g r i d points a l l o c a t e d . Multigrid Methods ( 4 2 ) * * * * * * ********************* * *********@@@@@@@<§@* * * * * * ***********@@@@@@@@@* * ***********@@@@@@@@@* * * * * * * **********@@@@@@@@@@* * **********@@@@@@@@<§@* * * * * * * *********@@@@@@@@@(§@* * *********@@@@@@@@@@@* * * * * * * *********@@@@@@@@@@@* * *********@@@@@@@@@@@* * * * * * * *********@@@@<a(a@<§<§(§@* * *********@@@@@@@@@@@* * * * * * * *********@@@@@@@@@@<a* * *********@@@@@@@@@@@* * * * * * * *********@@@@@@@@@@@* * *********@@@@@@@@@@@* * * * * * * *********<a@@@@@@@<§(a@* * *********@@@@@@@@@@@* * * * * * * *********@@<a<a<a@@@<§<§@* * *********@@@@@@@@@@@* * * * * * * *********@@@@@@@@@@@* * *********@@@@@@@@@@@* * * * * * * *********@@@@@@@@<§@@* * ***@@****@@@@@@@@@@@* * * * * * * ***@@@***@@@@@@@@@@@* * ***@@@@@*@@@@@@@@@@@* * * * * * **@@@@@@*@@@@@@@@@@@* * ***@@@@@@@@@@@@@@@@@* * * * * * *****@@@@@@@@@@@@@@@* * *@@@@@<§@@@@@@@@@* * *@@@@@@@@@@@@@@@* * **@@**@@@@@@@@@@* * * * * * * ********************* Figure 11: Level 32 grid, example (15) The grid at level n=32 used on example (15). This i s not a f u l l gr id , and the error estimates on the asterisk points are less than the tolerance, while those on the "@" points are greater than the tolerance. The blank areas of the grid have no grid points allocated. Multigrid Methods (43) * * * * * * * *************************************** * * * * * * * * * ***********@(a@@@@@@** * * * * * * * * * * * * * * * ****************@@@** * * * * * * * * * * ****************@@@@@@* * * * * * * * * * * * * * * * *@**@**********@@@@@@@* * * * * * * * * * *****************@@@@@@($j* * * * * * * * * * * * * * * ****************@@<§@<§@@@* * * * * * * * * * ****************@@@@@@@(§* * * * * * * * * * * * * * * *@**************@@@@@(§<3(§* * * * * * * * * * *<§*************@@@@@@(§<§(§* * * * * * * * * * * * * * * *@*************@@@@@@@@@* * * * * * * * * * *@*************@@@@@@@@@* * * * * * * * * * * * * * * *@**@(§(a********@@@@@@@@@* * * * * * * * * * ***@@@@@@******@@@@@@@@@* * * * * * * * * * * * * * * **@@@@@@@@@****@(a@(a@@@@(a* * * * * * * * * * **@@@@@@@@@@***@@@@@@@@@* * * * * * * * * * * * * * * **@@@@@@@@@@***@@@@@@@@@* * * * * * * * * * *@@@@@@@@@@@@*@@@@@@@@@@* * * * * * * * * * * * * * * *@@g@<a@@@@(a@(a*@@@@@@(a@@(a* * * * * * * * * * *@(§(§<§(§(a@(a@<§(a(a*(a@@@@@(§(§(§(§* * * * * * * * * * * * * * * *(a@@@@@@@@@@@@@@@@@@@@@@* * * * * * * * * * *@@@@@@@@@@@@@@@@@@@@@@@* * * * * * * * * ******* * * *@@@@@@@@@@<a@@@@(a@@@@(a@(a* * * * *@@@@@*** * *@@@@@@@@@@@@@@@@@@@@@@@* * * * * * * * * *@@@@@@@*****@@@@@@@@@@@@@@@@@@@@@@@* * * ***@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@(a* * * * * * * *@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@* * * ***@@@@@@@@@@@@@@@i§@@@@@@@@@@@@@@@@@@@* * * * * * * * *****@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@* * *@@@@@@@@@@@@@@@@@*@@@@@@@@@@@@@ * * *@@@@@@@@@@@@@@@ * * *@@@@@@@@@<§@@(a * * * * *@@@@@ * * *@@@ * * * * * *@@@@@@@(a@@(a@ * * * * * * * * *************************************** Figure 12: Level 64 gr i d , example (15) The g r i d at l e v e l n=64 used on example (15). Only h a l f of the rows i n t h i s g r i d are displayed here, due to space l i m i t a t i o n s . The erro r estimates on the a s t e r i s k points are le s s than the tolerance, while those on the "@" points are more than the tolerance. The blank areas of the g r i d have no gr i d points a l l o c a t e d . The width scale of the diagram has been doubled f o r t h i s diagram, to accommodate the larger number of points displayed here r e l a t i v e to the previous diagrams. Multigrid Methods ( 4 4 ) 4. Discussion 4.1. What has been done Several implementations of the m u l t i g r i d method f o r the s o l u t i o n of systems of p a r t i a l d i f f e r e n t i a l equations have been produced. The family of problems to which the l a t e s t of these implementations may be applied i s written as a u + c u + d u + e u + f u = g xx yy x y a where a ,c ,d ,e , and f are matrices, g i s a vector, and a l l are func-t i o n s of x and y. Some of these implementations use a data structure which i s very wasteful of storage, but i s simple to u t i l i z e . The l a t e s t version uses a l i s t based, space e f f i c i e n t data structure which allows adaptive g r i d i n t e r p o l a t i o n s . A v a r i e t y of combinations of relaxations and coarse g r i d c orrec-tions were t r i e d . Some of these were: i ) 2nd order d i f f e r e n c i n g , FAS coarse g r i d c o r r e c t i o n (CGC) i i ) 1 s t order d i f f e r e n c i n g when necessary, FAS CGC i i i ) 1 s t order d i f f e r e n c i n g when necessary, r e s i d u a l CGC i v ) 2nd order d i f f e r e n c i n g , with a r t i f i c i a l v i s c o s i t y The l a s t three options were used f o r m i l d l y s t i f f problems. We can define such a problem as one which, with the usual centered d i f -ferencing, would y i e l d a d i f f e r e n c e system which does not meet the diagonal dominance c r i t e r i o n on one or more of the grids , but which does meet the c r i t e r i o n on at l e a s t the f i n e s t g r i d . I t was found that, though i t i s advantageous to use the FAS method ( i t allows r extrapolation and easy adaptation to nonlinear problems), d i f f i c u l t i e s a r i s e when 1 s t order d i f f e r e n c i n g i s needed M u l t i g r i d Methods (45) on the c o r r e c t i o n grids, viz. the extra e r r o r i s propagated from the coarse grids up to the f i n e s t g r i d , r u i n i n g the convergence. I f r e s i d u a l c o r r e c t i o n i s used instead of FAS, the e r r o r on the coarse grids does improve according to theory, and i t does not propagate and r u i n the f i n e g r i d convergence. The error which propagates i n the FAS case i s 1st order with respect to u, while the propagating er r o r i n the r e s i d u a l c o r r e c t i o n case i s 1st order with respect to the error i n u. Fourth order accuracy i s attainable with r e s i d u a l c o r r e c t i o n , but adaptation to nonlinear problems w i l l not be straightforward. The use of a r t i f i c i a l v i s c o s i t y on the coarse grids was found to work we l l , but the extrapolation cannot be expected to work at the g r i d points using the v i s c o s i t y , so the approximation and e r r o r improvement on those gridpoints i s poor. Whether t h i s w i l l become a s i g n i f i c a n t d i f f i c u l t y when e n t i r e grids need v i s c o s i t y has yet to be examined. One c e r t a i n l y must worry about the "leverage" of error, where the e r r o r estimate on a given g r i d point may be s u f f i c i e n t l y low and may even be a correct estimate, but may cause a dispropor-t i o n a t e l y large e r r o r i n nearby g r i d points which are s t i l l being processed. In t h i s case the algorithm may proceed n i c e l y but y i e l d an i n v a l i d r e s u l t . 4.2. Take Home Message The i n d i c a t i o n s from t h i s study are that the use of a r t i f i c i a l v i s c o s i t y i s a p r a c t i c a l method f o r ensuring convergence i n the s o l u -t i o n of m i l d l y s t i f f boundary value problems by the m u l t i g r i d method. It i s easier to program than the one sided d i f f e r e n c e s approach and when used i n conjunction with the FAS c o r r e c t i o n scheme, i s more Multigrid Methods (46) robust with respect to the propagation of e r r o r from coarse to f i n e g r i d . This p o s s i b i l i t y of using the FAS algorithm i s a requirement for extending the implementation to nonlinear systems. For l i n e a r systems, the use of r e s i d u a l coarse g r i d corrections seems to be preferable. 4.3. Omissions from this work This p r o j e c t can be extended i n a number of ways which, using hindsight, are beyond the scope of t h i s t h e s i s : (1) Tune up the adaptive strategy used by the package. I f we l i m i t e d the proportion of g r i d points which could be i n t e r p o l a t e d at each l e v e l , the t o t a l work to be done would be more predictable than i the present system of using a simple tolerance on the e r r o r e s t i -mate. Also, the strategy should be more f l e x i b l e , so the system could change i t s mind about points that i t once considered to be solved. (2) Implement a version using 1st order d i f f e r e n c i n g i n combination with the FAS algorithm, or demonstrate that t h i s combination can-not work. There should be some way of preventing the 1st order coarse grids from i n t e r f e r i n g with the f i n e r 2nd order g r i d s . I t may even be that I have not examined the d i s t r i b u t i o n of e r r o r over the breadth of the g r i d i n s u f f i c i e n t d e t a i l . (3) I f (1) i s impractical, implement a r e s i d u a l c o r r e c t i o n version f o r use on l i n e a r problems. For some l i n e a r problems a one-sided algorithm may be advantageous, and i t i s always preferable to have the choice. Multigrid Methods (47) ( 4 ) Implement a v e r s i o n f o r n o n l i n e a r problems. T h i s would p o s s i b l y use t h e FAS c o r r e c t i o n scheme, and Newton i t e r a t i o n a t each g r i d p o i n t , as d i d t h e f i r s t i m p l e m e n t a t i o n d i s c u s s e d i n t h i s t h e s i s . ( 5 ) Implement some o f t h e r e l a x a t i o n t e c h n i q u e s s u g g e s t e d i n t h e l i t e r a t u r e , such as r e d - b l a c k o r d e r i n g Gauss S e i d e l (RBGS), l i n e Gauss S e i d e l ( L G S ) , and i n c o m p l e t e LU d e c o m p o s i t i o n ( I L U ) . One might i n c l u d e i n t h i s c a t e g o r y e x p e r i m e n t s i n t y p e s o f r e s t r i c -t i o n s , p r o l o n g a t i o n s , i n t e r p o l a t i o n s , and c o a r s e g r i d c o r r e c -t i o n s , a l l components o f t h e m u l t i g r i d method. ( 6) T r y some r e a l i s t i c problems, as opposed t o ones I have c o n -s t r u c t e d from known s o l u t i o n s . Some examples might be t h e Cauchy Riemann problem, t h e Steady S t a t e S t o k e s system, and t h e Steady S t a t e I n c o m p r e s s i b l e N a v i e r S t o k e s system. ( 7 ) T r y 9 p o i n t a v e r a g i n g r e s t r i c t i o n f o r problems w i t h o s c i l l a t o r y b e h a v i o u r . T h i s c o u l d be c o n s i d e r e d t o be a s u b s e t o f t h e a s p e c t o f t r y i n g out d i f f i c u l t problems t o gauge t h e power o f t h i s i m p l e m e n t a t i o n . T h i s t y p e o f r e s t r i c t i o n was implemented i n between t h e f i r s t two v e r s i o n s , b u t no c o n c l u s i v e r e s u l t s were o b t a i n e d , due t o t h e n a t u r e o f t h e t e s t p roblems. ( 8 ) E x t e n d t h e a l g o r i t h m t o h a n d l e a r b i t r a r y b o u n d a r i e s . T h i s would i n v o l v e s e t t i n g up boundary g r i d p o i n t s by b i n a r y s e a r c h u n t i l t h e boundary was s u f f i c i e n t l y c l o s e t o a p o i n t , and t h e n i n c l u d -i n g t h e s e boundary p o i n t s i n t h e c o a r s e g r i d s . Multigrid Methods (48) 5. References 1. E. J . A s s e l t , The m u l t i g r i d method and a r t i f i c i a l v i s c o s i t y , i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, W. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982 . 2. W. Auzinger and H. J . Ste t t e r , Defect Correction and m u l t i g r i d i t e r a t i o n s , i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, W. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 3. A. Brandt, M u l t i - l e v e l adaptive technique (MLAT) f o r f a s t numerical s o l u t i o n to boundary value problems. , i n Proceedings Third International Conference on Numerical Methods in Fluid Mechanics, Paris, 1972, v o l . 18, H. Cabannes and R. Teman (ed.), Springer-Verlag, B e r l i n , 1973, 82-89. 4. A. Brandt, M u l t i - l e v e l adaptive techniques (MLAT) I. The M u l t i -g r i d Method, Research Report RC 6026, IBM T.J. Watson Research Center, Yorktown Heights, NY, 1976. 5. A. Brandt, M u l t i - l e v e l adaptive solutions to p a r t i a l d i f f e r e n t i a l equations - Ideas and Software, i n Proceedings of Symposium on Mathematical Software, J. Rice (ed.), Academic Press, New York, March 1977, 277-318. 6. A. Brandt, M u l t i - l e v e l adaptive solutions to boundary value problems, J. Math. Comp. 31, (1977), 333-390. 7. A. Brandt, M u l t i - l e v e l Adaptive Techniques (MLAT) f o r Singular Perturbation Problems, i n Numerical Analysis of Singular Perturbation Problems, P. W. Hemker and J . J . H. M i l l e r (ed.), Academic Press, London, 1979. 8. A. Brandt, Stages i n Developing M u l t i g r i d Solutions, i n Numerical Methods for Engineering, E. Absi, R. Glowinski and P. Lascaux (ed.), Dunod, P a r i s , 1980, 23-44. 9. A. Brandt and S. Ta'asan, M u l t i - g r i d methods f o r hig h l y o s c i l l a t o r y problems, Research Report, Dept. of Applied Mathematics, Weizmann I n s t i t u t e of Science, Rehovot, 1981. 10. A. Brandt, A Guide to M u l t i g r i d Development, i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, w. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 11. J . E. Dendy, Black Box M u l t i g r i d , J. Comp. Phys. 48, 3 (1982), 366-386. 12. J. E. Dendy, Black Box M u l t i g r i d f o r NonSymmetric Problems, J. Appl. Math, and Comp. 13, (1983), 261. 13. H. Foerster and K. Witsch, M u l t i g r i d Software f o r the Solution of E l l i p t i c Problems on Rectangular Domains, i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, w. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 14. W. Hackbusch, On M u l t i g r i d I t e r a t i o n s with Defect Correction, i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, W. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 15. P. W. Hemker, On the Structure of an Adaptive M u l t i - l e v e l Algorithm. BIT, , 1980. Multigrid Methods ( 4 9 ) 16. P. w. Hemker, The Incomplete LU Decomposition as a Relaxation Method in Multigrid Algorithms Boundary and Interior Layers -Computational and Asymptotic Methods, Boole Press, Dublin, 1980. 17. R. K e t t l e r , Analysis and comparison of r e l a x a t i o n schemes i n robust m u l t i g r i d and preconditioned conjugate gradient methods., i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, W. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 18. W. J . A. Mol, Smoothing and Coaarse Grid Approximation Properties of M u l t i g r i d Methods, NW 110/81, Department of Numerical Mathematics, S t i c h t i n g Mathematisch Centrum, Amsterdam, 1981. 19. w. J . A. Mol, On the Choice of Suitable Operators and Parameters i n M u l t i g r i d Methods, NW 107/81, Department of Numerical Mathematics, S t i c h t i n g Mathematisch Centrum, Amsterdam, 1981. 20. R. A. Nicolaides, On the observed rate of convergence of an i t e r a t i v e method applied to a model e l l i p t i c defference equation, Math. Comp. 32, (1978), 127-133.. 21. K. Stuben and U. Trottenberg, M u l t i g r i d Methods: Fundamental Algorithms, Model Problem Analysis, and Applications, i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, w. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. 22. P. Wesseling, A convergence proof f or a multiple g r i d method, NA-21, D e l f t Technical University, D e l f t , The Netherlands., 1978. 23. P. Wesseling, T h e o r e t i c a l and p r a c t i c a l aspects of a m u l t i g r i d method, NA-37, Department of Mathematics, D e l f t U n i v e r s i t y of Technology, D e l f t , 1980. 24. P. Wesseling, A Robust and E f f i c i e n t M u l t i g r i d Method, i n Multigrid Methods Proceedings, Koln-Porz, November 1981, v o l . 960, W. Hackbush and U. Trottenberg (ed.), Springer-Verlag, B e r l i n , 1982. Multigrid Methods (50) 6. APPENDICES Appendix A: User Instructions The mgrid program i s used t o n u m e r i c a l l y s o l v e systems o f l i n e a r e l l i p t i c p a r t i a l d i f f e r e n t i a l e q u a t i o n s , as: a u + c u + d u + e u + f u = g xx yy x y y where a,c,d,e,f a r e m a t r i c e s , g i s a v e c t o r , and a l l a r e f u n c t i o n s o f x and y, The u s e r must p r o v i d e t h e s e f u n c t i o n s as p a s c a l f u n c t i o n s and l i n k them t o g e t h e r w i t h t h e r e s t o f t h e package( See appendix I I f o r sample s o u r c e ) . The m a k e f i l e and s e v e r a l f u n c t i o n examples i n the package s o u r c e code can a c t as examples. A sample ru n command on t h e u n i x system used f o r development i s : mgrida t i t l e < p r m f i l e > r e s u l t s f i l e The t i t l e i s p r i n t e d on t h e o u t p u t , and i s o p t i o n a l . The s t a n d a r d i n p u t c o n t a i n s t h e e x e c u t i o n p a r ameters, one p e r l i n e , i n t h e f o l l o w i n g o r d e r : nmin minimum g r i d s i z e used ( a m u l t i p l e o f 2) nmax maximum g r i d s i z e used ( a l a r g e r m u l t i p l e o f 2) a l p h a a, t h e r a t i o e x p e c t e d o f t a u from g r i d t o g r i d ( u s u a l l y 0.25 ) d e l t a 6, t h e r a t i o e x p e c t e d o f e r r o r from g r i d t o g r i d ( u s u a l l y use 0.125) e t a slow convergence c r i t e r i o n 7?, u s u a l l y use 0.5, (0 < e t a < 1) t o l t o l e r a n c e on e r r o r e s t i m a t e , an e x i t c r i t e r i o n . t a u t o l t o l e r a n c e on t h e change i n r , an e x i t c r i t e r i o n w i t h i n t h e s o l v e r o u t i n e . c y c l e l i m i t a n o t h e r e x i t c r i t e r i o n w i t h i n t h e s o l v e r o u t i n e , i f r r a t i o does not r e a c h a l p h a w i t h i n t h i s number o f c y c l e s , s o l v e s t o p s t r y i n g , a c c e p t i n g t h e v a l u e o f r so f a r o b t a i n e d . e x t r a p t a u t h e r e x t r a p o l a t i o n f a c t o r , u s u a l l y s e t t o 1.33 l b x t h e lower bound o f t h e x component o f t h e domain ubx t h e upper bound o f t h e x component o f t h e domain l b y t h e lower bound o f t h e y component o f t h e domain uby t h e upper bound o f t h e y component o f t h e domain In t h e parameter f i l e , comments may be p l a c e d on t h e l i n e a f t e r t h e v a l u e . U n l i k e t h i s example, a d d i t i o n a l comment l i n e s may not be i n s e r t e d . S t a n d a r d o u t p u t c o n t a i n s t h e l e v e l by l e v e l summaries o f t h e p r o c e s s o f t h e computation. The maximum v a l u e s o f r and e r r o r e s t i m a t e s a r e p r i n t e d each time t h e y a r e computed a t t h e f i n e g r i d l e v e l . A t t h e c o n c l u s i o n o f each g r i d l e v e l c o m putation i n t h e n e s t e d s e t , a d e t a i l e d o u t p u t f i l e i s c r e a t e d w i t h t h e name " g r i d . n n n " , where t h e "nnn" i s t h e g r i d l e v e l . T h i s f i l e c o n t a i n s a r e c o r d f o r each g r i d p o i n t i n t h e g r i d , w i t h i t s x,y, and u v a l u e s , as w e l l as t h e e r r o r e s t i m a t e s and t h e a c t u a l e r r o r s . T h i s i m p l e m e n t a t i o n computes t h e a c t u a l e r r o r u s i n g a s u p p l i e d f u n c t i o n which computes Multigrid Methods (51) the actual solution, a feature which would probably not be carried over into a production model, since in real l i f e an analytical solution i s rarely known. Appendix B contains sample source for the required functions, appendix C a sample parameter f i l e , and appendix 4 sample output for that particular run. M u l t i g r i d Methods (52) Appendix B: Sample Source Code Below i s a l i s t i n g o f an example o f the s o u r c e which must be s u p p l i e d by t h e u s e r . Of c o u r s e t h e p r o b l e m must f i t i n t o t h e f a m i l y o f problems t o which t h i s package may be a p p l i e d , l i n e a r second o r d e r systems o f e q u a t i o n s . The r o u t i n e s i n c l u d e a boundary v a l u e f u n c t i o n ( f b c ), a s o l u t i o n f u n c t i o n ( u s o l ) , and s i x f u n c t i o n s d e f i n i n g t h e problem, a , c , d , e , f , and g . # i n c l u d e " d e f s . h " ttinclude " e x t . h " p r o c e d u r e f b c ; { ( x , y : r e a l ; v a r v : v e c t o r ) ; } c o n s t Re = 16; b e g i n v [ l ] := ( e x p ( - R e ) - e x p ( R e * ( x - l ) ) ) * ( l - e x p ( - R e * y ) ) ; v [ 2 ] := exp(x+y); end; p r o c e d u r e u s o l ; { ( x , y : r e a l ; v a r v : v e c t o r ) ; } c o n s t Re = 16; b e g i n v [ l ] := ( e x p ( - R e ) - e x p ( R e * ( x - l ) ) ) * ( l - e x p ( - R e * y ) ) ; v [ 2 ] := exp(x+y); end; p r o c e d u r e a; { ( x , y : r e a l ; v a r v : m a t r i x ) ; } c o n s t Re = 16; b e g i n v [ l ] [ l ] := 1/Re; v [ l ] [ 2 ] := 0; v [ 2 ] [ l ] := 0; v [ 2 ] [ 2 ] •= 1; end; p r o c e d u r e c; { ( x , y : r e a l ; v a r v : m a t r i x ) ; } c o n s t Re = 16; b e g i n v [ l ] [ l ] := 1/Re; v [ l ] [ 2 ] := 0; v [ 2 ] [ l ] := 0; v [ 2 ] [ 2 ] := 1; end; Multigrid Methods procedure d; {(x,y:real; var vsmatrix);} begin v [ l ] [ l ] := -1; v [ l ] [ 2 ] : = 0; v [ 2 ] [ l ] := 0; v[2][2] != 0; end; procedure e; {(x,y:real; var v:matrix);} begin v [ l ] [ l ] := 1; v [ l ] [ 2 ] := 0; v [ 2 ] [ l ] := 0; v[2][2] := 0; end; procedure f; {(x,y:real; var v:matrix);} begin v [ l ] [ l ] := 0; v [ l ] [ 2 ] := 0; v [ 2 ] [ l ] := 0; v[2][2] := -2; end; procedure g; {(x,y:real; var v:vector);} begin v [ l ] := 0; v[2] := 0; end; M u l t i g r i d Methods (54) Appendix C: A Sample Parameter F i l e Below i s a sample o f a parameter f i l e . T h i s parameter f i l e and the s o u r c e i n t h e p r e v i o u s appendix were used t o produce t h e r e s u l t s i n t h e next appendix. 4 nmin 32 nmax 0.25 a l p h a 0.125 d e l t a 0.5 e t a 1.0e-8 t o l 0.01 t o l e r a n c e on t h e change i n t a u v a l u e 3 c y c l e l i m i t 0.1 lowerbndx 0.9 upperbndx 0.1 lowerbndy 0.9 upperbndy Multigrid Methods (55) Appendix D: A Sample Output F i l e Here i s t h e d e s c r i p t i v e o u t p u t f o r t h e r u n s p e c i f i e d by t h e p r e v i o u s 2 a p p e n d i c e s : t h e p a r a m e t e r s were: d e l t a = 1.250e-01, alpha= 2.500e-01, Eta= 5.000e-01 e r r o r and t a u t o l e r a n c e s : 1.000e-04 1.000e-02, l i m i t on # c y c l e s 3 nmin= 4, nmax= 32 The upper and lower bounds o f x: 0.000e+00 1.000e+00 The upper and lower bounds o f y: 0.000e+00 1.000e+00 p r i n t e r r : n= 4, max e r r o r 1.535e-02, a t x= 7.500e-01, and y= 2.500e-p r i n t e r r : n= 4, max eg e r r o r 4.478e-04, a t x= 5.000e-01, and y= 5.000e s o l v e : n= 8 maximum t a u v a l u e [ 1 ] 3 . 2 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : n= 8 maximum t a u v a l u e [ 2 ] 3 . 4 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 7 s o l v e : n= 8 maximum t a u v a l u e [ l ] 3 . 2 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : n= 8 maximum t a u v a l u e [ 2 ] 3 . 4 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 9 s o l v e : n= 8 maximum t a u v a l u e [ l ] 3 . 2 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : n= 8 maximum t a u v a l u e [ 2 ] 3 . 4 e - 0 2 , a t x= 7.5e-01 y= 7.5e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 11 s o l v e : t a u t o l e r a n c e a t t a i n e d p r i n t e r r e s t : f o r g r i d n= 4, max e r r o r c o r r e c t i o n i s 3.959e-03 p r i n t e r r e s t : f o r g r i d n= 2, max e r r o r c o r r e c t i o n i s 3.959e-03 p r i n t e r r : n= 8, max e r r o r 1.150e-01, a t x= 8.750e-01, and y= 2.500e-p r i n t e r r : n= 8, max eg e r r o r 1.772e-02, a t x= 7.500e-01, and y= 2.500e s o l v e : n= 16 maximum t a u v a l u e [ l ] 8 , 7 e - 0 1 , a t x= 8.8e-01 y= 3.8e-01 s o l v e : n= 16 maximum t a u v a l u e [ 2 ] 1 . 8 e - 0 2 , a t x= 8.8e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 2 s o l v e : n= 16 maximum t a u v a l u e [ l ] 7 . 7 e - 0 1 , a t x= 8.8e-01 y= 5.0e-01 s o l v e : n= 16 maximum t a u v a l u e [ 2 ] l . l e - 0 2 , a t x= 8.8e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 5 Multigrid Methods (56) s o l v e : n= 16 maximum t a u v a l u e [ l ] 7 . 4 e - 0 1 , a t x= 8.8e-01 y= 8.8e-01 s o l v e : n= 16 maximum t a u v a l u e [ 2 ] l . l e - 0 2 , a t x= 8.8e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 7 s o l v e : r e a c h e d l i m i t o f number o f c y c l e s p r i n t e r r e s t : f o r g r i d n= 8, max e r r o r c o r r e c t i o n i s 1.025e-01 p r i n t e r r e s t : f o r g r i d n= 4, max e r r o r c o r r e c t i o n i s 9.661e-02 p r i n t e r r : n= 16, max e r r o r 1.007e-02, a t x= 9.375e-01, and y= 8.125e-01 p r i n t e r r : n= 16, max eg e r r o r 4.789e-03, a t x= 8.750e-01, and y= 5.000e-01 s o l v e : n= 32 maximum t a u v a l u e [ l ] 4 . l e - 0 1 , a t x= 9.4e-01 y= 6.3e-01 s o l v e : n= 32 maximum t a u v a l u e [ 2 ] l . 2 e + 0 0 , a t x= 5.0e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 4 s o l v e : n= 32 maximum t a u v a l u e [ l ] 4 . 0 e - 0 1 , a t x= 9.4e-01 y= 8.8e-01 s o l v e : n= 32 maximum t a u v a l u e [ 2 ] l . l e + 0 0 , a t x= 5.0e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 6 s o l v e : n= 32 maximum t a u v a l u e [ l ] 4 . 0 e - 0 1 , a t x= 9.4e-01 y= 6.3e-01 s o l v e : n= 32 maximum t a u value[2]1.2e+00, a t x= 5.0e-01 y= 8.8e-01 s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) 8 s o l v e : r e a c h e d l i m i t o f number o f c y c l e s p r i n t e r r e s t : f o r g r i d n= 16, max e r r o r c o r r e c t i o n i s 8.108e-03 p r i n t e r r e s t : f o r g r i d n= 8, max e r r o r c o r r e c t i o n i s 8.069e-03 p r i n t e r r : n= 32, max e r r o r 9.936e-04, a t x= 9.688e-01, and y= 9.063e-01 p r i n t e r r : n= 32, max eg e r r o r 5.366e-04, a t x= 5.000e-01, and y= 8.750e-01 mgrid: e x e c u t i o n f i n i s h e d , u s i n g f i n e s t g r i d 32 Multigrid Methods (57) Appendix E: Source for the central routines solve, relax, taufunc I f t h e u s e r wishes t o a p p l y a new d i f f e r e n t i a l o p e r a t o r , o r use a d i f f e r e n t c o a r s e g r i d c o r r e c t i o n a l g o r i t h m , new v e r s i o n s o f t h e r e l a x a t i o n r o u t i n e and t h e c o a r s e g r i d c o r r e c t i o n r o u t i n e would have t o be w r i t t e n . P o s s i b l y some changes would be r e q u i r e d i n t h e s o l v e r o u t i n e as w e l l , i n t h e p a r t s c oncerned w i t h t h e c o a r s e g r i d c o r r e c t i o n r o u t i n e . Here i s the s o u r c e f o r t h o s e t h r e e r o u t i n e s . p r e c e d e d by t h e d e f i n i t i o n s f i l e used by a l l r o u t i n e s . l a b e l 999; c o n s t TINY = 1.0e-30; BIG = l . O e l O ; GRIDSIZE = 16384; { max g r i d d i v i s i o n : 2~14 : keep i , j t o 2 b y t e s } VECTORSIZE = 2 ; { s h o u l d be t h e o n l y change f o r more e q u a t i o n s } FIELDSPERPAGE = 8 ; { j u s t used f o r f o r m a t t i n g o u t p u t } EXTRAPTAU = 1.3333333333333333; { e x t r a p o l a t i o n f a c t o r } UNKNOWN = 1.0E35; { use t h i s t o t e s t f o r unknown t a u value} t y p e p o s i t i v e = l . . m a x i n t ; n o n n e g a t i v e = 0..maxint; g r i d r a n g e = 0..GRIDSIZE; v e c t o r = array[1..VECTORSIZE] o f r e a l ; i v e c t o r = a r ray[l..VECTORSIZE] o f i n t e g e r ; m a t r i x = array[1..VECTORSIZE] o f v e c t o r ; a t v l i s t = ~ v l i s t ; v l i s t = r e c o r d u : v e c t o r ; l i n k : a t v l i s t end; a t g r i d = " g r i d t y p ; g r i d t y p = r e c o r d i , j : g r i d r a n g e ; u, g r h s , t a u , e r r e s t : v e c t o r ; smrate, { g r i d i s made o f t h e s e nodes} { t h e r e s u l t v e c t o r } { r i g h t hand s i d e f o r t h i s g r i d p o i n t } { t a u v a l u e a t l a s t CGC, v a l i d on c o a r s e g r i d , V a l u e from the CCG ( c o a r s e c o a r s e g r i d , a r e t o r i g h t o f t h e CCG g r i d p o i n t s } { e r r o r e s t i m a t e , c o u l d be f l a g o r s c a l a r i f we wanted t o save memory?} { smoothing r a t e , see i f f u r t h e r r e l a x i n g i s needed, c o u l d be a f l a g ? } Multigrid Methods (58) r l x c h a n g e : r e a l ; { l a s t change which produced u i n t h e r e l a x a t i o n s , a l s o used as a f l a g i f n e g a t i v e , t h i s p o i n t has not been i n t e r p o l a t e d and s h o u l d not be used f o r f u r t h e r i n t e r p o l a t i o n } t a u f l a g : b o o l e a n ; { i s t a u s u f f i c i e n t a t t h i s p o i n t ? } x l i n k , y l i n k : a t g r i d { p o i n t e r s t o r i g h t ( x ) and up(y) d i r e c t i o n s } end; l o n g s t r i n g = a r r a y [ 1 . . 1 0 0 ] o f c h a r ; v a r { t h e e x t e r n a l d e c l a r a t i o n s , s i g n i f i e d by c a p i t a l f i r s t l e t t e r as opposed t o c o n s t a n t s which a r e a l l caps. } Nmin,Nmax : 1..GRIDSIZE; C y c l e l i m i t : n o n n e g a t i v e ; Alpha,Delta,Eta,Tol,Tautol,Lbndx,Ubndx,Lbndy,Ubndy,Hx,Hy : r e a l ; Z ero : v e c t o r ; { used f o r z e r o i n g v e c t o r s } M a i n g r i d : a t g r i d ; { p o i n t s t o g r i d p o i n t a t x=Lbx,y=Iiby} T i t l e : l o n g s t r i n g ; { R o u t i n e which does t h e main work o f t h e m u l t i g r i d method f o r s o l v i n g a d i s c r e t i z e d form o f a boundary v a l u e o r d i n a r y d i f f e r e n t i a l e q u a t i o n } # i n c l u d e " d e f s . h " # i n c l u d e " e x t . h " p r o c e d u r e s o l v e ; { ( n : g r i d r a n g e ; g r : a t g r i d ; c o r r g r i d f l a g : b o o l e a n ) ; } v a r x,y,change,conv,t : r e a l ; i n t 2 , i n t 4 , 1 , c o u n t , c y c l e c o u n t : n o n n e g a t i v e ; maxi,maxj : i v e c t o r ; s u b g r , n e x t x , n e x t y , s u b n e x t x , s u b n e x t y , b e l o w , a b o v e , r i g h t , l e f t : a t g r i d ; u l i s t , s a v e d g r i d : a t v l i s t ; f i r s t , t f l a g , r e t u r n f l a g , t a u f l a g : b o o l e a n ; maxtau,prevtau : v e c t o r ; Multigrid Methods (59) b e g i n { 0 } i f n=Nmin t h e n s o l v e d ( g r ) { s o l v e d i r e c t l y and t o h i g h a c c u r a c y } e l s e b e g i n { 1 } c o p y t a u ( n , g r ) ; { copy t a u e x p e c t a t i o n on CCG t o a d j a c e n t g r i d p o i n t s and i n i t i a l i z e g r i d f l a g s } f o r 1:= 1 t o VECTORSIZE do p r e v t a u [ l ] := BIG; c y c l e c o u n t := 1; f i r s t : = t r u e ; change:=BIG; { s i g n a l f o r f i r s t sweep } i n t 2 != 2*(GRIDSIZE d i v n ) ; i n t 4 : = 2 * i n t 2 ; count:=0; r e t u r n f l a g := f a l s e ; w h i l e not r e t u r n f l a g do b e g i n { 2 } change:=BIG; { s i g n a l f o r f i r s t sweep } conv:=0; w h i l e (conv < E t a ) do b e g i n { 3 } r e l a x ( g r , n, change, c o n v ) ; count:=count+l; end; { 3 } t a u f u n c ( n , g r , f i r s t ) ; { g e t t a u f o r each g r i d p o i n t } f i r s t : = f a l s e ; r e s t r i c t ( n , g r , s u b g r ) ; { i n c l u d i n g t h e b o u n d a r i e s ; c r e a t e subgr} f o r 1 := 1 t o VECTORSIZE do m a x t a u [ l ] := 0; t a u f l a g := t r u e ; n e x t y : = g r * . y l i n k ; s u b n e x t y : = s u b g r ~ . y l i n k ; below:=gr; w h i l e s u b n e x t y * . y l i n k <> n i l do { boundary v a l u e s not used l a t e r } b e g i n { 3 } w h i l e n e x t y * . i <> s u b n e x t y * . i do b e g i n b e l o w : = n e x t y ; n e x t y : = n e x t y * . y l i n k ; { 4 } end; { 4 } a b o v e : = n e x t y * . y l i n k ; i f a b o v e * . i - b e l o w * . i = i n t 2 t h e n b e g i n { 4 } l e f t : = n e x t y ; nextx := n e x t y * . x l i n k ; subnextx := s u b n e x t y * . x l i n k ; w h i l e s u b n e x t x * . x l i n k <> n i l do b e g i n { 5 } w h i l e n e x t x * . j <> s u b n e x t x * . j do b e g i n l e f t : = n e x t x ; n e x t x : = n e x t x * . x l i n k ; { 6 } end; { 6 } w h i l e b e l o w * . j < n e x t x * . j do below:=below*.xlink; i f b e l o w * . j = n e x t x * . j t h e n b e g i n { 6 } ab o v e : = n e x t x * . y l i n k ; r i g h t : = n e x t x * . x l i n k ; Multigrid Methods (60) i f ( above".! - b e l o w " . i = i n t 2 ) and ( r i g h t " . j - l e f t " . j = i n t 2 ) t h e n b e g i n { we a r e a t a g r i d p o i n t i n a f i n e g r i d r e g i o n } { 7 } i f t a u f l a g t h e n i f ( n e x t x " . j mod i n t 4 = 0) and ( n e x t x " . i mod i n t 4 = 0) t h e n f o r 1 := 1 t o VECTORSIZE do i f ( a b s ( n e x t x " . t a u [ l ] ) > a b s ( n e x t x " . x l i n k " . t a u [ l ] )) t h e n t a u f l a g : = f a l s e ; i f c o r r g r i d f l a g t h e n f o r 1 := 1 t o VECTORSIZE do s u b n e x t x " . g r h s [ l ] := n e x t x " . t a u [ l ] + n e x t x " . g r h s [ 1 ] e l s e f o r 1 := 1 t o VECTORSIZE do s u b n e x t x " . g r h s [ 1 ] := EXTRAPTAU*nextx".tau[1] + n e x t x " . g r h s [ 1 ] ; f o r 1;=1 t o VECTORSIZE do b e g i n { 8 } t := a b s ( n e x t x " . t a u [ l ] ) ; i f t > m a x t a u f l ] t h e n b e g i n { 9 } m a x t a u [ l ] := t ; maxi[1] := n e x t x " . i ; maxj[1] := n e x t x " . j ; end; { 9 } end; { 8 } end; { 7 } end; { 6 } subnextx := s u b n e x t x " . x l i n k ; end; { 5 } end; { 4 } s u b n e x t y : = s u b n e x t y " . y l i n k ; end; { 3 } s a v e g r i d ( s u b g r , s a v e d g r i d ) ; { save uH f o r ops a f t e r s o l v i n g } { n e c e s s a r y i f r e s t r i c t i o n i s not i n j e c t i o n } { m e r e l y c o p i e s g r i d u v a l u e s i n t o a l i s t t o be r e a d back i n a g a i n , t h e g r i d s t r u c t u r e w i l l not change b e f o r e i t i s used again} i f not c o r r g r i d f l a g t h e n f o r 1 := 1 t o VECTORSIZE do b e g i n { 3 } x := Lbndx + maxj[1] * Hx; y := Lbndy + m a x i f l ] * Hy; w r i t e I n ; w r i t e l n ( ' s o l v e : n= ' , n : l , ' maximum t a u v a l u e [ ' , 1 : 0 , ' ] * , m a x t a u f l ] : 4 , ' , a t x= ',x:4, ' y= ' , y : 4 ) ; end; { 3 } s o l v e ( n d i v 2 , s u b g r , t r u e ) ; Multigrid Methods (61) nexty:=subgr; u l i s t : = s a v e d g r i d ; w h i l e n e x t y < > n i l do b e g i n { p u t t h e c o r r e c t i o n i n t o t h e g r i d } { 3 } nextx:=nexty; { i n s t e a d o f t h e r e s i d u a l s o l u t i o n } w h i l e n e x t x o n i l do b e g i n { 4 } f o r 1 := 1 t o VECTORSIZE do nextx*.u[1]:=nextx*.u[1] - u l i s t * . u [ l ] ; u l i s t : = u l i s t " . l i n k ; n e x t x : = n e x t x * . x l i n k ; end; { 4 } n e x t y : = n e x t y * . y l i n k ; end; { 3 } u l i s t : = s a v e d g r i d ; { g e t r i d o f t h e save s t o r a g e } w h i l e u l i s t < > n i l do b e g i n { 3 } s a v e d g r i d : = s a v e d g r i d " . l i n k ; d i s p o s e ( u l i s t ) ; u l i s t : = s a v e d g r i d ; end; { 3 } p r o l o n g r i d ( g r , s u b g r ) ; { p r o l o n g a t e t h e g r i d t o l e v e l n } { g r i s used t o ensure i d e n t i c a l s t r u c t u r e o f b o t h g r i d s f o r e f f i c i e n c y c o u l d a v o i d expanding subgr by d o i n g t h e c o r r e c t i o n a t t h e same s t e p , i . e . compute t h e f i n e g r i d c o r r e c t i o n s as one goes. T h i s would i n v o l v e combining p r o l o n g r i d and t h e f o l l o w i n g c o r r e c t i o n code} n e x t y : = g r " . y l i n k ; s u b n e x t y : = s u b g r " . y l i n k ; w h i l e n e x t y " . y l i n k o n i l do { uH = 0 a t b o u n d a r i e s , so don't add them } b e g i n { 3 } n e x t x : = n e x t y " . x l i n k ; s u b n e x t x : = s u b n e x t y " . x l i n k ; w h i l e nextx* . x l i n k o n i l do b e g i n { 4 } i f n e x t x * . r l x c h a n g e >= 0 th e n { s i g n a l t h a t t h i s p o i n t } { i s b e i n g r e l a x e d } f o r 1 := 1 t o VECTORSIZE do b e g i n { 5 } n e x t x * . u [ l ] := n e x t x * . u [ l ] + subnextx*.u[1]; i f not c o r r g r i d f l a g t h e n n e x t x * . e r r e s t [ l ] := n e x t x * . e r r e s t [ l ] + a b s ( s u b n e x t x * . u [ 1 ] ) / ( l + a b s ( n e x t x * . u [ 1 ] ) ) ; { r e l a t i v e e r r o r f o r u > > 1 } end; { 5 } s u b n e x t x : = s u b n e x t x * . x l i n k ; n e x t x : = n e x t x * . x l i n k ; end; { 4 } s u b n e x t y : = s u b n e x t y * . y l i n k ; n e x t y : = n e x t y * . y l i n k ; end; { 3 } d i s p o s e g r i d ( s u b g r ) ; Multigrid Methods (62) i f not c o r r g r i d f l a g t h e n w r i t e l n ( ' s o l v e : r e l a x a t i o n c o u n t ( t o p l e v e l t o t a l ) ', c o u n t : 1 ) ; i f t a u f l a g t h e n { s u c c e s s f u l l y a t t a i n e d t a u g o a l } b e g i n { 3 } i f not c o r r g r i d f l a g t h e n w r i t e l n ( ' s o l v e : t a u g o a l a t t a i n e d ' ) ; r e t u r n f l a g := t r u e ; end { 3 } e l s e b e g i n { 3 } t f l a g : = t r u e ; f o r 1:= 1 t o VECTORSIZE do i f abs( l - a b s ( m a x t a u [ l ] / p r e v t a u [ l ] ) ) > T a u t o l t h e n t f l a g : = f a l s e ; i f t f l a g t h e n b e g i n { 4 } i f not c o r r g r i d f l a g t h e n w r i t e l n ( ' s o l v e : t a u t o l e r a n c e a t t a i n e d ' ) ; r e t u r n f l a g := t r u e ; end { 4 } e l s e i f c y c l e c o u n t = C y c l e l i m i t t h e n b e g i n { 4 } i f not c o r r g r i d f l a g t h e n w r i t e l n ( ' s o l v e : r e ached l i m i t o f number o f c y c l e s ' ) ; r e t u r n f l a g : = t r u e ; end { 4 } e l s e c y c l e c o u n t := c y c l e c o u n t + 1 ; prevtau:=maxtau; end; { 3 } end; { 2 } end; { 1 } end; Multigrid Methods (63) # i n c l u d e " d e f s . h " ttinclude " e x t . h " { T h i s i s a r o u t i n e t o c a l c u l a t e t h e d i f f e r e n c e between one element o f t h e v e c t o r r e s u l t i n g from the a p p l i c a t i o n o f t h e d i f f e r e n c e o p e r a t o r on g r i d l e v e l n/2 and g r i d l e v e l n, on t h e v e c t o r u. n i s t h e f i n e g r i d l e v e l . T h i s f u n c t i o n c a l c u l a t e s t h e d i f f e r e n c e between t h e g r i d o p e r a t o r v e c t o r p r o d u c t a t a g i v e n g r i d p o i n t , between l e v e l n and l e v e l n/2. h h h H h h So t a u := - I L u + L I u H H The o p e r a t o r i s t h e d i s c r e t i z e d o p e r a t o r c o r r e s p o n d i n g t o t h e p e r t u r b e d l a p l a c e e q u a t i o n a*u + c*u + d*u + e*u + f * u = g xx y y x y The o p e r a t o r does not use g, and i n t h e d i f f e r e n c e t h e c o e f f i c i e n t o f u, i . e . t h e f f u n c t i o n , c a n c e l , so t h e y a r e not used h e r e . T h i s r o u t i n e presumes t h a t t h e r e s t r i c t i o n o p e r a t o r i s i n j e c t i o n , so i t j u s t uses t h e c o a r s e g r i d v a l u e s p r e s e n t , not computing them. In t h i s r o u t i n e o n l y f i n e g r i d r e g i o n s a r e p r o c e s s e d , i . e . , c o a r s e g r i d p o i n t s which a r e surrounded by f i n e g r i d p o i n t s . } p r o c e d u r e t a u f u n c ; { ( n : g r i d r a n g e ; g r : a t g r i d ; f i r s t : b o o l e a n ) ; } v a r n e x t y , c b e l o w , b e l o w , a b o v e , c a b o v e , l e f t , c l e f t , r i g h t , c r i g h t , h e r e : a t g r i d ; 1 , i n t , i n t 2 : i n t e g e r ; a m a t r i x , c a m a t r i x , c m a t r i x , c c m a t r i x , d m a t r i x , e m a t r i x : m a t r i x ; t v e c t , t v e c t l , t v e c t 2 , t x x , t x x c , t y y , t y y c , t x , t x c , t y , t y c : v e c t o r ; t a , t c , t d , t e , t : r e a l ; x,y,hx,hy,hhx,hhy : r e a l ; b e g i n { 0 } i n t := GRIDSIZE d i v n; i n t 2 := 2 * i n t ; hx := i n t * H x ; hhx:=hx*hx; hy := i n t * H y ; hhy:=hy*hy; { need 5 i n a f i n e g r i d row on t h e l e f t boundary} nexty:=gr; cbelow:=nexty; below:=cbelow".ylink; c l e f t : = b e l o w " . y l i n k ; a b o v e : = c l e f t " . y l i n k ; cabove:=above".ylink; Multigrid Methods (64) w h i l e cabove<>nil do b e g i n { 1 } i f ( above".! - b e l o w " . i = i n t 2 ) t h e n { f i n e g r i d region} b e g i n l e f t : = c l e f t " . x l i n k ; h e r e : = l e f t " . x l i n k ; r i g h t : = h e r e " . x l i n k ; { 2 } c r i g h t : = r i g h t " . x l i n k ; { x l i n e i s s e t up} w h i l e c r i g h t < > n i l do b e g i n { 3 } i f ( r i g h t " . j - l e f t " . j = i n t 2 ) t h e n { 5 i n a row} b e g i n { 4 } w h i l e c b e l o w " . j < h e r e " . j do cbelow:=cbelow".xlink; i f ( c b e l o w " . j = h e r e " . j ) and ( c b e l o w " . y l i n k <> h e r e ) t h e n b e g i n a b o v e : = h e r e " . y l i n k ; below:=cbelow".ylink; { 5 } i f ( a b o v e " . i - b e l o w " . i = i n t 2 ) t h e n b e g i n cabove:=above".ylink; { 6 } y := here".i*Hy+Lbndy; x := h e r e " . j *Hx+Lbndx; a ( x , y , a m a t r i x ) ; { f o r m a t r i x computations p r o c e d u r e s } c ( x , y , c m a t r i x ) ; { used below have p r e f i x e s : } d ( x , y , d m a t r i x ) ; { mt s i g n i f i e s m a t r i x o p e r a t i o n , } e ( x , y , e m a t r i x ) ; { v t s i g n i f i e s v e c t o r o p e r a t i o n } m t c o p y ( a m a t r i x , c a m a t r i x ) ; m t c o p y ( c m a t r i x , c c m a t r i x ) ; f o r 1 := 1 t o VECTORSIZE do b e g i n { 7 } t a := a m a t r i x [ l ] [ 1 ] ; t c := c m a t r i x [ l ] [ 1 ] ; t d := d m a t r i x [ l ] [ 1 ] ; t e := e m a t r i x f 1 ] [ 1 ] ; t := a b s ( t d * h x / 2 ) ; i f t a < 2*t t h e n c a m a t r i x [ l ] [ 1 ] := 2 * t ; i f t a < t t h e n a m a t r i x [ l ] [ 1 ] := t ; t := a b s ( t e * h y / 2 ) ; i f t c < 2*t t h e n c c m a t r i x [ l ] [ 1 ] := 2 * t ; i f t c < t t h e n c m a t r i x [ 1 ] [ 1 ] : = t ; t x x [ l ] := ( l e f t " . u [ l ] - 2 * h e r e " . u [ l ] + r i g h t " . u [ 1 ] ) / h h x ; t x x c [ l ] := ( c l e f t " . u [ l ] - 2 * h e r e " . u [ l ] + c r i g h t " . u [ l ] )/(4*hhx); t y y [ l ] := ( b e l o w " . u [ l ] - 2 * h e r e " . u [ l ] + above".u[1])/hhy; t y y c [ l ] := ( c b e l o w " . u [ l ] - 2 * h e r e " . u [ l ] + c a b o v e " . u [ l ] )/(4*hhy); t x [ l ] := ( r i g h t " . u [ l ] - l e f t " . u [ l ] ) / ( 2 * h x ) ; t x c [ l ] := ( c r i g h t " . u [ l ] - c l e f t " . u [ l ] )/(4*hx); t y [ l ] := ( a b o v e " . u [ l ] - b e l o w " . u [ l ] ) / ( 2 * h y ) ; t y c [ l ] := ( c a b o v e " . u [ l ] - c b e l o w " . u [ l ] ) / ( 4 * h y ) ; end; { 7 } Multigrid Methods (65) i f f i r s t then right".tau := here".tau; mtvtmult(camatrix,txxc,tvectl); mtvtmult(amatrix,txx,tvect2); vtsub(tvectl,tvect2,here".tau); mtvtmult(ccmatrix,tyyc,tvectl); mtvtmult(cmatrix,tyy,tvect2); vtsub(tvectl,tvect2,tvect); vtadd(tvect,here".tau,here".tau); vtsub(txc,tx,tvect); mtvtmult(dmatrix,tvect,tvectl); vtadd( tvectl,here".tau,here". tau); vtsub(tyc,ty,tvect); mtvtmult(ematrix,tvect,tvectl); vtadd(tvectl,here".tau,here".tau); end; { 6 } end; { 5 } end; { 4 } i f ( l e f t " . j - c l e f t " . j = int) then begin cleft:=here; left:=right; here:=cright; right:=cright".xlink; i f right <> n i l then cright:=right".xlink else cright:=nil; end else begin cleft:=left; left:=here; here:=right; right:=cright; cright:=cright".xlink; end; end; { 3 } end; { 2 } i f ( nexty".ylink".i - nexty".i = int) then begin cbelow := nexty".ylink".ylink; nexty:=cbelow; below:=cbelow".ylink; cleft:=below".ylink; above:=cleft".ylink; i f above<> n i l then cabove:=above".ylink else cabove:=nil; end else begin cbelow := nexty".ylink; nexty:=cbelow; below:=cbelow".ylink; cleft:=below".ylink; above:=cleft".ylink; cabove:=above".ylink; end; end; { 1 } end; { 0 } Multigrid Methods (66) # i n c l u d e " d e f s . h " # i n c l u d e " e x t . h " { r o u t i n e which does one G a u s s - S e i d e l r e l a x a t i o n sweep, a c t u a l l y s o l v i n g t h e n e c e s s a r y system o f n e q u a t i o n s f o r t h e g r i d p o i n t f o r t h e n s i m u l t a n e o u s pdes. d e s i g n e d t o be used i n s o l v i n g t h e g e n e r a l p r o b l e m au + cu + d*u + e*u + f * u = g xx y y x y where a , c , d , e , f , and g a r e f u n c t i o n s o f ( x , y ) , } p r o c e d u r e r e l a x ; { ( g r : a t g r i d ; n : g r i d r a n g e ; v a r change,conv : r e a l ) ; } v a r n e x t x , n e x t y , r i g h t , l e f t , a b o v e , b e l o w : a t g r i d ; 1 : 1..VECTORSIZE; i n t , i n t 2 : g r i d r a n g e ; t m p , n e w c h g , t h i s c h a n g e , h x , h y , h h x , h h y , x , y , t t a , t t c , t t d , t t e , h 2 x , h 2 y : r e a l ; t a , t c , t d , t e , a m a t r i x , c m a t r i x , d m a t r i x , e m a t r i x , f m a t r i x : m a t r i x ; b e t a 0 , b e t a l , b e t a 2 , b e t a 3 , b e t a 4 : m a t r i x ; rhs,uprev,wl,w2,w3,w4 : v e c t o r ; b e g i n newchg:=0; int:=GRIDSIZE d i v n; hx:=int*Hx; hy:=int*Hy; hhx:=hx*hx; hhy:=hy*hy; h2x:=2*hx; h2y:=2*hy; i n t 2 : = 2 * i n t ; n e x t y : = g r " . y l i n k ; below:=gr; w h i l e n e x t y " . y l i n k < > n i l do b e g i n n e x t x : = n e x t y * . x l i n k ; l e f t : = n e x t y ; w h i l e n e x t x " . x l i n k < > n i l do b e g i n r i g h t : = n e x t x " . x l i n k ; a b o v e : = n e x t x " . y l i n k ; below:=below".xlink; i f b e l o w " . y l i n k <> nextx t h e n b e l o w : = b e l o w l i n k ( g r , n e x t x ) ; i f ( a b o v e " . i - b e l o w " . i = i n t 2 ) and ( r i g h t " . j - l e f t " . j = i n t 2 ) t h e n b e g i n v t c o p y ( n e x t x " . u , u p r e v ) ; x:=nextx".j *Hx+Lbndx; y : = n e x t x " . i *Hy+Lbndy; a ( x , y , a m a t r i x ) ; c ( x , y , c m a t r i x ) ; d ( x , y , d m a t r i x ) ; e ( x , y , e m a t r i x ) ; f ( x , y , f m a t r i x ) ; Multigrid Methods (67) { p u t i n t h e a r t i f i c i a l v i s c o s i t y } f o r 1 := 1 t o VECTORSIZE do b e g i n t t a := a m a t r i x [ l ] [ l ] / h h x ; t t d := d m a t r i x [ l ] [ l ] / h 2 x ; t t c := c m a t r i x [ l ] [ l ] / h h y ; t t e := e m a t r i x [ l ] [ l ] / h 2 y ; i f t t a < a b s ( t t d ) t h e n a m a t r i x [ l ] [ 1 ] := h h x * a b s ( t t d ) ; i f t t c < a b s ( t t e ) t h e n c m a t r i x [ l ] [ 1 ] := h h y * a b s ( t t e ) ; end; s c m t m u l t ( - 2 / h h x , a m a t r i x , t a ) ; scintmult( -2/hhy, c m a t r i x , t c ) ; m t a d d ( t a , t c , b e t a O ) ; m t a d d ( b e t a O , f m a t r i x , b e t a O ) ; { betaO i s th e m a t r i x f o r our g r i d p o i n t } s c m t m u l t ( 0 . 5 / h x , d m a t r i x , t d ) ; s c m t m u l t ( 1 / h h x , a m a t r i x , t a ) ; m t a d d ( t a , t d , b e t a l ) ; m t s u b ( t a , t d , b e t a 2 ) ; s c m t m u l t ( 1 / h h y , c m a t r i x , t c ) ; s c m t m u l t ( 0 . 5 / h y , e m a t r i x , t e ) ; m t a d d ( t c , t e , b e t a 3 ) ; m t s u b ( t c , t e , b e t a 4 ) ; { b e t a [ 1 - 4 ] a r e f o r t h e r h s o f t h i s g r i d p o i n t } m t v t m u l t ( b e t a l , r i g h t " . u , w l ) ; m t v t m u l t ( b e t a 2 , l e f t " . u , w 2 ); mtvtmult(beta3,above".u,w3 ); mtvtmult(beta4,below".u,w4); v t s u b ( n e x t x " . g r h s , w l , r h s ) ; v t s u b ( r h s , w 2 , r h s ) ; v t s u b ( r h s , w 3 , r h s ) ; v t s u b ( r h s , w 4 , r h s ) ; { c a l l t h e g a u s s i a n e l i m i n a t i o n r o u t i n e t o s o l v e t h e system } gauss ( betaO,rhs, nextx".u ); t h i s c h a n g e := 0; f o r 1 := 1 t o VECTORSIZE do b e g i n tmp := a b s ( n e x t x " . u [ l ] - u p r e v f l ] ) / ( l + a b s ( n e x t x " . u [ 1 ] ) ) ; i f tmp > t h i s c h a n g e t h e n thischange:=tmp; end; Multigrid Methods (68) i f nextx".rlxchange > 0 then nextx".smrate;=thischange/nextx".rlxchange; i f nextx".rlxchange >= 0 then { don't muck up the boundary flag} nextx".rlxchange:=thischange; i f thischange > newchg then newchg:=thischange; end; l e ft:=nextx; nextx:=nextx".xlink; end; below:=nexty; nexty:=nexty".ylink; end; i f newchg=0 then begin conv:=1; w r i t e l n ( ' r e l a x : newchg = 0'); end else i f (change > 1) or (newchg < 1 ) then conv:=newchg/change else i f (change/newchg = 0 ) then begin conv:=l; w r i t e l n ( ' r e l a x : conv would have overflowed'); end e l s e conv:=newchg/change; change := newchg; end;
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use.
Graduate
On the implementation of multigrid methods for the numerical solution of partial differential equations
Text