f19-6

862 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).standardtridiagonalalgorithm. Given un,onesolves(19.5.36)for un+1/2,substitutes on the right-hand side of (19.5.37), and then solves for un+1. The key question is how to choose the iteration parameter r, the analog of a choice of timestep for an initial value problem. As usual, the goal is to minimize the spectral radius of the iteration matrix. Although it is beyond our scope to go into details here, it turns out that, for the optimal choice of r, the ADI method has the same rate of convergence as SOR. The individual iteration steps in the ADI method are much more complicated than in SOR, so the ADI method would appear to be inferior. This is in fact true if we choose the same parameter rfor every iteration step. However, it is possible to choose a different r for each step. If this is done optimally, then ADI is generally more efficient than SOR. We refer you to the literature [1-4]for details. Our reason for not fully implementing ADI here is that, in most applications, it has been superseded by the multigrid methods described in the next section. Our advice is to use SOR for trivial problems (e.g., 20×20), or for solving a larger problem once only, where ease of programming outweighs expense of computertime. Occasionally, the sparse matrix methods of §2.7 are useful for solving a set of difference equations directly. For production solution of large elliptic problems, however, multigrid is now almost always the method of choice. CITED REFERENCES AND FURTHER READING: Hockney, R.W., and Eastwood, J.W. 1981, Computer Simulation Using Particles (New York: McGraw-Hill), Chapter 6. Young, D.M. 1971, Iterative Solution of LargeLinearSystems (NewYork: Academic Press). [1] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §§8.3–8.6. [2] Varga, R.S. 1962, Matrix Iterative Analysis (Englewood Cliffs, NJ: Prentice-Hall). [3] Spanier, J. 1967, in Mathematical Methods for Digital Computers, Volume 2 (New York: Wiley), Chapter 11. [4] 19.6 Multigrid Methods for Boundary Value Problems Practicalmultigridmethodswerefirstintroducedinthe1970sbyBrandt. These methods can solve elliptic PDEs discretized on Ngrid points in O(N)operations. The “rapid” direct elliptic solvers discussed in §19.4 solve special kinds of elliptic equations in O(NlogN)operations. The numerical coefficients in these estimates are such that multigrid methods are comparable to the rapid methods in executionspeed. Unlike the rapid methods, however,the multigridmethods can solve general elliptic equations with nonconstant coefficients with hardly any loss in efficiency. Even nonlinear equations can be solved with comparable speed. Unfortunately there is not a single multigrid algorithm that solves all elliptic problems. Rather there is a multigrid technique that provides the framework for solvingtheseproblems. Youhaveto adjustthevariouscomponentsofthealgorithm within this framework to solve your specific problem. We can only give a brief 19.6MultigridMethodsforBoundaryValueProblems 863Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).introduction to the subject here. In particular, we will give two sample multigrid routines, one linear and one nonlinear. By following these prototypes and byperusing the references [1-4], you should be able to develop routines to solve your own problems. Therearetworelated,butdistinct,approachestotheuseofmultigridtechniques. Thefirst,termed“themultigridmethod,”isameansforspeedinguptheconvergence of a traditional relaxation method, as defined by you on a grid of pre-specified fineness. Inthis case, youneed defineyourproblem(e.g.,evaluateits sourceterms) only on this grid. Other, coarser, grids defined by the method can be viewed as temporary computational adjuncts. The second approach,termed (perhapsconfusingly)“the full multigrid(FMG) method,” requires you to be able to define your problem on grids of various sizes (generallybydiscretizingthesameunderlyingPDEintodifferent-sizedsetsoffinite-difference equations). In this approach, the method obtains successive solutions on finer and finer grids. You can stop the solution either at a pre-specified fineness, or you can monitor the truncation error due to the discretization, quitting only whenit is tolerably small. Inthissectionwewillfirstdiscussthe“multigridmethod,”thenusetheconcepts developed to introduce the FMG method. The latter algorithm is the one that we implement in the accompanying programs. FromOne-Grid,throughTwo-Grid,toMultigrid The key idea of the multigrid method can be understood by considering the simplest case of a two-grid method. Suppose we are trying to solve the linear elliptic problem Lu=f (19.6.1 ) where Lissomelinearellipticoperatorand fisthesourceterm. Discretizeequation (19.6.1) on a uniform grid with mesh size h. Write the resulting set of linear algebraic equations as Lhuh=fh (19.6.2 ) Let/tildewideuhdenote some approximate solution to equation (19.6.2). We will use the symbol uhto denote the exact solution to the difference equations (19.6.2). Then theerrorin/tildewideuhor thecorrection is vh=uh−/tildewideuh (19.6.3 ) Theresidualordefectis dh=Lh/tildewideuh−fh (19.6.4 ) (Beware: someauthorsdefineresidualasminusthedefect,andthereisnotuniversal agreement about which of these two quantities 19.6.4 defines.) Since Lhis linear, the error satisfies Lhvh=−dh (19.6.5 ) 864 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).At this point we need to make an approximationto Lhin orderto find vh. The classical iteration methods, such as Jacobi or Gauss-Seidel, do this by finding, ateach stage, an approximate solution of the equation /hatwideL h/hatwidevh=−dh (19.6.6 ) where /hatwideLhis a “simpler” operator than Lh. For example, /hatwideLhis the diagonal part of Lhfor Jacobi iteration, or the lower triangle for Gauss-Seidel iteration. The next approximation is generated by /tildewideunew h=/tildewideuh+/hatwidevh (19.6.7 ) Now consider, as an alternative, a completely different type of approximation forLh, one in which we “coarsify” rather than “simplify.” That is, we form some appropriate approximation LHofLhon a coarser grid with mesh size H(we will always take H=2h,butotherchoicesarepossible). Theresidualequation(19.6.5) is now approximated by LHvH=−dH (19.6.8 ) Since LHhas smaller dimension, this equationwill be easier to solve than equation (19.6.5). To define the defect dHon the coarse grid, we need a restriction operator Rthat restricts dhto the coarse grid: dH=Rdh (19.6.9 ) The restriction operator is also called the fine-to-coarse operator or theinjection operator. Oncewe havea solution /tildewidevHto equation(19.6.8),we needa prolongation operator Pthat prolongates or interpolates the correction to the fine grid: /tildewidevh=P/tildewidevH (19.6.10 ) The prolongation operator is also called the coarse-to-fine operator or theinter- polation operator . Both RandPare chosen to be linear operators. Finally the approximation /tildewideuhcan be updated: /tildewideunew h=/tildewideuh+/tildewidevh (19.6.11 ) One step of this coarse-grid correction scheme is thus: Coarse-GridCorrection •Compute the defect on the fine grid from (19.6.4). •Restrict the defect by (19.6.9). •Solve (19.6.8) exactly on the coarse grid for the correction. •Interpolate the correction to the fine grid by (19.6.10). 19.6MultigridMethodsforBoundaryValueProblems 865Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).•Compute the next approximation by (19.6.11). Let’scontrasttheadvantagesanddisadvantagesofrelaxationandthecoarse-grid correctionscheme. Considertheerror vhexpandedintoadiscreteFourierseries. Call the componentsin the lower half of the frequencyspectrumthe smoothcomponents and the high-frequencycomponents the nonsmoothcomponents . We have seen that relaxationbecomesveryslowlyconvergentinthelimit h→0,i.e.,whentherearea largenumberofmeshpoints. Thereasonturnsouttobethatthesmoothcomponentsare only slightly reduced in amplitude on each iteration. However, many relaxation methods reduce the amplitude of the nonsmooth components by large factors on each iteration: They are good smoothing operators . For the two-grid iteration, on the other hand, components of the error with wavelengths <∼2Hare not even representable on the coarse grid and so cannot be reducedto zero on this grid. But it is exactly these high-frequencycomponentsthat can be reduced by relaxation on the fine grid! This leads us to combine the ideas of relaxation and coarse-grid correction: Two-GridIteration •Pre-smoothing: Compute ¯u hby applying ν1≥0steps of a relaxation method to /tildewideuh. •Coarse-grid correction: As above, using ¯uhto give ¯unew h. •Post-smoothing: Compute /tildewideunew hbyapplying ν2≥0stepsoftherelaxation method to ¯unew h. It is only a short step from the above two-grid method to a multigrid method. Instead of solving the coarse-grid defect equation (19.6.8) exactly, we can get an approximatesolutionofitbyintroducinganevencoarsergridandusingthetwo-griditeration method. If the convergencefactor of the two-gridmethodis small enough, we will need only a few steps of this iteration to get a good enough approximate solution. We denote the number of such iterations by γ. Obviously we can apply this idea recursively down to some coarsest grid. There the solution is found easily, for example by direct matrix inversion or by iterating the relaxation schemeto convergence. One iteration of a multigrid method, from finest grid to coarser grids and back to finest grid again, is called a cycle. The exact structure of a cycle depends on the value of γ, the number of two-grid iterations at each intermediate stage. The case γ=1iscalledaV-cycle,while γ=2iscalledaW-cycle(seeFigure19.6.1). These are the most important cases in practice. Note that once more than two grids are involved,the pre-smoothingsteps after the first one on the finest grid need an initial approximation for the error v. This should be taken to be zero. Smoothing,Restriction,and ProlongationOperators The most popular smoothing method, and the one you should try first, is Gauss-Seidel,sinceitusuallyleadstoagoodconvergencerate. Ifweorderthemesh points from 1 to N, then the Gauss-Seidel scheme is ui=−/parenleftBig N/summationdisplay j=1 j /negationslash=iLijuj−fi/parenrightBig1 Liii=1,...,N (19.6.12 ) 866 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).E γ = 2 γ = 12-grid 3-grid 4-gridSS SSS S ESS SS ESS S ES SS ESSS SESS ES SS S ES SS ESSS Figure 19.6.1. Structure of multigrid cycles. S denotes smoothing, while E denotes exact solution on the coarsest grid. Each descending line \denotes restriction ( R) and each ascending line /denotes prolongation ( P). Thefinest grid is at the top level of each diagram. For the V-cycles ( γ=1) the E step is replaced by one 2-grid iteration each time the number of grid levels is increased by one. For the W-cycles ( γ=2), each E step gets replaced by two 2-grid iterations. wherenewvaluesof uareusedontheright-handsideastheybecomeavailable. The exactformoftheGauss-Seidelmethoddependsontheorderingchosenforthemeshpoints. For typical second-orderelliptic equations like our model problem equation (19.0.3), as differenced in equation (19.0.8), it is usually best to use red-black ordering,makingonepassthroughthemeshupdatingthe “even”points(likethered squares of a checkerboard) and another pass updating the “odd”points (the black squares). When quantities are more strongly coupled along one dimension thananother, one should relax a whole line along that dimension simultaneously. Line relaxation for nearest-neighborcoupling involves solving a tridiagonal system, and so is still ef ficient. Relaxing oddand evenlines on successive passes is called zebra relaxation and is usually preferred over simple line relaxation. Notethat SOR should notbe usedas a smoothingoperator. Theoverrelaxation destroysthe high-frequencysmoothingthat is so crucialforthe multigridmethod. A succint notationforthe prolongationandrestrictionoperatorsis to givetheir symbol. The symbol of Pis found by considering v Hto be 1 at some mesh point (x, y), zero elsewhere, and then asking for the values of PvH. The most popular prolongationoperatoris simple bilinearinterpolation. It gives nonzerovalues at the 9 points (x, y),(x+h, y),..., (x−h, y−h), where the values are 1,1 2,...,1 4. 19.6MultigridMethodsforBoundaryValueProblems 867Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Its symbol is therefore  1 41 21 4 1 211 2 1 41 21 4  (19.6.13 ) The symbolof Ris definedby considering vhto be definedeverywhereonthe fine grid, and then asking what is Rvhat(x, y)as a linear combination of these values. Thesimplestpossiblechoicefor Risstraightinjection ,whichmeanssimply filling each coarse-grid point with the value from the corresponding fine-grid point. Its symbol is “[1].”However, dif ficulties can arise in practice with this choice. It turnsoutthatasafechoicefor Ristomakeittheadjointoperatorto P.T od efinethe adjoint,de finethescalar productoftwo gridfunctions uhandvhformeshsize has /angbracketleftuh|vh/angbracketrighth≡h2/summationdisplay x,yuh(x, y)vh(x, y)( 19.6.14 ) Then the adjoint of P, denoted P†,i sd efined by /angbracketleftuH|P†vh/angbracketrightH=/angbracketleftPuH|vh/angbracketrighth (19.6.15 ) Nowtake Ptobebilinearinterpolation,andchoose uH=1at(x, y),zeroelsewhere. SetP†=Rin (19.6.15) and H=2h. You will find that (Rvh)(x,y )=1 4vh(x, y)+1 8vh(x+h, y)+1 16vh(x+h, y+h)+··· (19.6.16 ) so that the symbol of Ris  1 161 81 16 1 81 41 8 1 161 81 16  (19.6.17 ) Note the simple rule: The symbolof Ris1 4the transposeof the matrixde finingthe symbolof P,equation(19.6.13). Thisruleisgeneralwhenever R=P†andH=2h. Theparticularchoiceof Rin(19.6.17)iscalled fullweighting . Anotherpopular choice for Rishalf weighting ,“halfway”between full weighting and straight injection. Its symbol is  01 80 1 81 21 8 01 80  (19.6.18 ) A similar notation can be used to describe the difference operator Lh.F o r example, the standard differencing of the model problem, equation (19.0.6), isrepresented by the five-point difference star L h=1 h2 01 0 1−41 01 0  (19.6.19 ) 868 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).If you are confronted with a new problem and you are not sure what PandR choices are likely to work well, here is a safe rule: Suppose mpis the order of the interpolation P(i.e.,itinterpolatespolynomialsofdegree mp−1exactly). Suppose mristheorderof R,andthat Ris theadjointofsome P(notnecessarilythe Pyou intend to use). Then if mis the order of the differential operator Lh, you should satisfy the inequality mp+mr>m. For example, bilinear interpolation and its adjoint, full weighting, for Poisson ’s equationsatisfy mp+mr=4>m =2. Of course the PandRoperators should enforce the boundary conditions for your problem. The easiest way to do this is to rewrite the difference equation to have homogeneousboundary conditions by modifying the source term if necessary(cf.§19.4). Enforcing homogeneous boundary conditions simply requires the P operator to produce zeros at the appropriate boundary points. The corresponding Ris then found by R=P †. FullMultigridAlgorithm So far we have described multigrid as an iterative scheme, where one starts with some initial guess on the finest grid and carries out enough cycles (V-cycles, W-cycles, ...) to achieve convergence. This is the simplest way to use multigrid: Simply apply enough cycles until some appropriate convergence criterion is met.However,ef ficiencycanbeimprovedbyusingthe FullMultigridAlgorithm (FMG), also known as nested iteration . Instead of starting with an arbitrary approximation on the finest grid (e.g., u h=0), thefirst approximation is obtained by interpolating from a coarse-grid solution: uh=PuH (19.6.20 ) Thecoarse-gridsolutionitselfis foundbya similarFMG processfromevencoarser grids. Atthecoarsestlevel,youstartwiththeexactsolution. Ratherthanproceedas inFigure19.6.1,then,FMGgetstoitssolutionbyaseriesofincreasinglytall “N’s,” each taller one probing a finer grid (see Figure 19.6.2). Note that Pin (19.6.20) need not be the same Pused in the multigrid cycles. It should be at least of the same order as the discretization Lh, but sometimes a higher-order operator leads to greater ef ficiency. It turns out that you usually need one or at most two multigrid cycles at each level before proceeding down to the next finer grid. While there is theoretical guidance on the required number of cycles (e.g., [2]), you can easily determine it empirically. Fix the finest level and study the solution values as you increase the numberofcyclesperlevel. Theasymptoticvalueofthesolutionistheexactsolution of the difference equations. The difference between this exact solution and thesolution for a small number of cycles is the iteration error. Now fix the number of cyclestobelarge,andvarythenumberoflevels,i.e.,thesmallestvalueof hused. In thiswayyoucanestimatethetruncationerrorforagiven h. Inyourfinalproduction code, there is no point in using more cycles than you need to get the iteration error down to the size of the truncation error. The simple multigrid iteration (cycle) needs the right-hand side fonly at the finestlevel. FMGneeds fatalllevels. Iftheboundaryconditionsarehomogeneous, 19.6MultigridMethodsforBoundaryValueProblems 869Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).4-grid ncycle = 1 4-grid ncycle = 2S SS SSS S SSS S E ESSS S S E EESSSSS EESSS ESSS S S EES SSS ESSS S S ESSS EESSS SSS S S E Figure 19.6.2. Structure of cycles for the full multigrid (FMG) method. This method starts on the coarsest grid, interpolates, and then re fines (by“V’s”), the solution onto grids of increasing fineness. you can use fH=Rfh. This prescription is not always safe for inhomogeneous boundaryconditions. In that case it is better to discretize fon each coarse grid. NotethattheFMGalgorithmproducesthesolutiononalllevels. Itcantherefore be combined with techniques like Richardson extrapolation. We now give a routine mglinthat implements the Full Multigrid Algorithm for a linear equation, the model problem (19.0.6). It uses red-black Gauss-Seidel as the smoothing operator, bilinear interpolation for P, and half-weighting for R. To change the routine to handle another linear problem, all you need do is modify the subroutines relax,resid, and slvsmlappropriately. A feature of the routine is the dynamical allocation of storage for variables de fined on the various grids. The subroutine malocemulates the Cfunction malloc. It allows you to write subroutines that operate on two-dimensionalarrays in the usual way, but to allocatestorage for these arrays in the calling program “on thefly”out of a single long one-dimensional array. SUBROUTINE mglin(u,n,ncycle) INTEGER n,ncycle,NPRE,NPOST,NG,MEMLEN DOUBLE PRECISION u(n,n)PARAMETER (NG=5,MEMLEN=13*2**(2*NG)/3+14*2**NG+8*NG-100/3)PARAMETER (NPRE=1,NPOST=1) C USES addint,copy,fill0,interp,maloc,relax,resid,rstrct,slvsml Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem(19.0.6). On input u(1:n,1:n) contains the right-hand side ρ, while on output it returns the solution. The dimension nis related to the number of grid levels used in the solution, NGbelow, by n= 2** NG+1.ncycle is the number of V-cycles to be used at each level. Parameters: NGis the number of grid levels used; MEMLEN is the maximum amount of memory that can be allocated by calls to maloc ;NPRE andNPOST are the number of relaxation sweeps before and after the coarse-grid correction is computed. INTEGER j,jcycle,jj,jpost,jpre,mem,nf,ngrid,nn,ires(NG), * irho(NG),irhs(NG),iu(NG),maloc DOUBLE PRECISION z 870 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).COMMON /memory/ z(MEMLEN),mem Storage for grid functions is allocated by maloc from array z. mem=0 nn=n/2+1 ngrid=NG-1irho(ngrid)=maloc(nn**2) Allocate storage for r.h.s. on grid NG−1, call rstrct(z(irho(ngrid)),u,nn) and fill it by restricting from the fine grid. 1 if (nn.gt.3) then Similarly allocate storage and fill r.h.s. on all coarse grids. nn=nn/2+1 ngrid=ngrid-1 irho(ngrid)=maloc(nn**2) call rstrct(z(irho(ngrid)),z(irho(ngrid+1)),nn) goto 1endif nn=3 iu(1)=maloc(nn**2)irhs(1)=maloc(nn**2) call slvsml(z(iu(1)),z(irho(1))) Initial solution on coarsest grid. ngrid=NGdo 16j=2,ngrid Nested iteration loop. nn=2*nn-1 iu(j)=maloc(nn**2) irhs(j)=maloc(nn**2)ires(j)=maloc(nn**2) call interp(z(iu(j)),z(iu(j-1)),nn) Interpolate from coarse grid to next finer grid. if (j.ne.ngrid) then call copy(z(irhs(j)),z(irho(j)),nn) Set up r.h.s. else call copy(z(irhs(j)),u,nn) endifdo 15jcycle=1,ncycle V-cycle loop. nf=nn do12jj=j,2,-1 Downward stoke of the V. do11jpre=1,NPRE Pre-smoothing. call relax(z(iu(jj)),z(irhs(jj)),nf) enddo 11 call resid(z(ires(jj)),z(iu(jj)),z(irhs(jj)),nf) nf=nf/2+1call rstrct(z(irhs(jj-1)),z(ires(jj)),nf) Restriction of the residual is the next r.h.s. call fill0(z(iu(jj-1)),nf) Zero for initial guess in next relaxation. enddo 12 call slvsml(z(iu(1)),z(irhs(1))) Bottom of V: solve on coarsest grid. nf=3do 14jj=2,j Upward stroke of V. nf=2*nf-1 call addint(z(iu(jj)),z(iu(jj-1)),z(ires(jj)),nf) Use res for temporary storage inside addint . do13jpost=1,NPOST Post-smoothing. call relax(z(iu(jj)),z(irhs(jj)),nf) enddo 13 enddo 14 enddo 15 enddo 16 call copy(u,z(iu(ngrid)),n) Return solution in u. return END SUBROUTINE rstrct(uc,uf,nc) INTEGER nc DOUBLE PRECISION uc(nc,nc),uf(2*nc-1,2*nc-1) Half-weighting restriction. ncis the coarse-grid dimension. The fine-grid solution is input inuf(1:2*nc-1,1:2*nc-1) , the coarse-grid solution is returned in uc(1:nc,1:nc) . INTEGER ic,if,jc,jf 19.6MultigridMethodsforBoundaryValueProblems 871Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do12jc=2,nc-1 Interior points. jf=2*jc-1 do11ic=2,nc-1 if=2*ic-1uc(ic,jc)=.5d0*uf(if,jf)+.125d0*(uf(if+1,jf)+ * uf(if-1,jf)+uf(if,jf+1)+uf(if,jf-1)) enddo 11 enddo 12 do13ic=1,nc Boundary points. uc(ic,1)=uf(2*ic-1,1) uc(ic,nc)=uf(2*ic-1,2*nc-1) enddo 13 do14jc=1,nc uc(1,jc)=uf(1,2*jc-1) uc(nc,jc)=uf(2*nc-1,2*jc-1) enddo 14 returnEND SUBROUTINE interp(uf,uc,nf) INTEGER nfDOUBLE PRECISION uc(nf/2+1,nf/2+1),uf(nf,nf)INTEGER ic,if,jc,jf,nc Coarse-to-fine prolongation by bilinear interpolation. nfis the fine-grid dimension. The coarse-grid solution is input as uc(1:nc,1:nc) ,w h e r e nc =nf/2+1 .T h efi n e - g r i d solution is returned in uf(1:nf,1:nf) . nc=nf/2+1 do12jc=1,nc Do elements that are copies. jf=2*jc-1do 11ic=1,nc uf(2*ic-1,jf)=uc(ic,jc) enddo 11 enddo 12 do14jf=1,nf,2 Do odd-numbered columns, interpolating ver- tically. do13if=2,nf-1,2 uf(if,jf)=.5d0*(uf(if+1,jf)+uf(if-1,jf)) enddo 13 enddo 14 do16jf=2,nf-1,2 Do even-numbered columns, interpolating hor- izontally. do15if=1,nf uf(if,jf)=.5d0*(uf(if,jf+1)+uf(if,jf-1)) enddo 15 enddo 16 returnEND SUBROUTINE addint(uf,uc,res,nf) INTEGER nf DOUBLE PRECISION res(nf,nf),uc(nf/2+1,nf/2+1),uf(nf,nf) C USES interp Does coarse-to-fine interpolation and adds result to uf.nfis the fine-grid dimension. The coarse-grid solution is input as uc(1:nc,1:nc) ,w h e r e nc =nf/2+1 .T h efi n e - g r i d solution is returned in uf(1:nf,1:nf) .res(1:nf,1:nf) is used for temporary storage. INTEGER i,j call interp(res,uc,nf) do12j=1,nf do11i=1,nf uf(i,j)=uf(i,j)+res(i,j) enddo 11 enddo 12 returnEND 872 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE slvsml(u,rhs) DOUBLE PRECISION rhs(3,3),u(3,3) C USES fill0 Solution of the model problem on the coarsest grid, where h=1 2. The right-hand side is input in rhs(1:3,1:3) and the solution is returned in u(1:3,1:3) . DOUBLE PRECISION h call fill0(u,3) h=.5d0u(2,2)=-h*h*rhs(2,2)/4.d0 return END SUBROUTINE relax(u,rhs,n) INTEGER n DOUBLE PRECISION rhs(n,n),u(n,n) Red-black Gauss-Seidel relaxation for model problem. The current value of the solution u(1:n,1:n) is updated, using the right-hand side function rhs(1:n,1:n) . INTEGER i,ipass,isw,j,jsw DOUBLE PRECISION h,h2h=1.d0/(n-1)h2=h*h jsw=1 do 13ipass=1,2 Red and black sweeps. isw=jsw do12j=2,n-1 do11i=isw+1,n-1,2 Gauss-Seidel formula. u(i,j)=0.25d0*(u(i+1,j)+u(i-1,j)+u(i,j+1) * +u(i,j-1)-h2*rhs(i,j)) enddo 11 isw=3-isw enddo 12 jsw=3-jsw enddo 13 return END SUBROUTINE resid(res,u,rhs,n) INTEGER nDOUBLE PRECISION res(n,n),rhs(n,n),u(n,n) Returns minus the residual for the model problem. Input quantities are u(1:n,1:n) and rhs(1:n,1:n) , while res(1:n,1:n) is returned. INTEGER i,jDOUBLE PRECISION h,h2i h=1.d0/(n-1) h2i=1.d0/(h*h)do 12j=2,n-1 Interior points. do11i=2,n-1 res(i,j)=-h2i*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)- * 4.d0*u(i,j))+rhs(i,j) enddo 11 enddo 12 do13i=1,n Boundary points. res(i,1)=0.d0res(i,n)=0.d0 res(1,i)=0.d0 res(n,i)=0.d0 enddo 13 returnEND 19.6MultigridMethodsforBoundaryValueProblems 873Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE copy(aout,ain,n) INTEGER n DOUBLE PRECISION ain(n,n),aout(n,n) Copies ain(1:n,1:n) toaout(1:n,1:n) . INTEGER i,j do12i=1,n do11j=1,n aout(j,i)=ain(j,i) enddo 11 enddo 12 returnEND SUBROUTINE fill0(u,n) INTEGER n DOUBLE PRECISION u(n,n) Fills u(1:n,1:n) with zeros. INTEGER i,j do12j=1,n do11i=1,n u(i,j)=0.d0 enddo 11 enddo 12 returnEND FUNCTION maloc(len) INTEGER maloc,len,NG,MEMLEN PARAMETER (NG=5,MEMLEN=13*2**(2*NG)/3+14*2**NG+8*NG-100/3) formglin C PARAMETER (NG=5,MEMLEN=17*2**(2*NG)/3+18*2**NG+10*NG-86/3) formgfas , N.B.! INTEGER mem DOUBLE PRECISION z COMMON /memory/ z(MEMLEN),mem Dynamical storage allocation. Returns integer pointer to the starting position for len array elements in the array z. The preceding array element is filled with the value of len,a n d the variable mem is updated to point to the last element of zthat has been used. if (mem+len+1.gt.MEMLEN) pause ’insufficient memory in maloc’z(mem+1)=len maloc=mem+2 mem=mem+len+1returnEND The routine mglinis written for clarity, not maximum ef ficiency, so that it is easy to modify. Several simple changes will speed up the executiontime: •Thedefect dhvanishesidenticallyatallblackmeshpointsafterared-black Gauss-Seidel step. Thus dH=Rdhforhalf-weightingreduces to simply copyinghalfthedefectfromthe finegridtothecorrespondingcoarse-grid point. The calls to residfollowed by rstrctin thefirst part of the V-cycle can be replaced by a routine that loops only over the coarse grid, filling it with half the defect. •Similarly, the quantity /tildewideunew h=/tildewideuh+P/tildewidevHneed not be computed at red mesh points, since they will immediately be rede fined in the subsequent Gauss-Seidel sweep. This means that addintneed only loop over black points. 874 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).•You can speed up relaxin several ways. First, you can have a special form when the initial guess is zero, and omit the routine fill0. Next, youcanstore h2fhonthevariousgridsandsaveamultiplication. Finally, it is possible to save an additionin the Gauss-Seidel formulaby rewriting it with intermediate variables. •On typical problems, mglinwith ncycle =1will return a solution with the iteration error bigger than the truncation error for the given size of h. To knock the error down to the size of the truncation error, you have to setncycle =2or, more cheaply, npre =2. A more ef ficient way turns out to be to use a higher-order Pin (19.6.20)than the linear interpolation used in the V-cycle. Implementing all the above features typically gives up to a factor of two improvementin executiontime and is certainlyworthwhilein a productioncode. NonlinearMultigrid: The FAS Algorithm Now turn to solving a nonlinear elliptic equation, which we write symbolically as L(u)=0 ( 19.6.21 ) Any explicit source term has been moved to the left-hand side. Suppose equation (19.6.21) is suitably discretized: Lh(uh)=0 ( 19.6.22 ) We will see below that in the multigrid algorithm we will have to consider equations where a nonzero right-hand side is generated during the course of the solution: Lh(uh)=fh (19.6.23 ) OnewayofsolvingnonlinearproblemswithmultigridistouseNewton ’smethod,which produces linear equations for the correction term at each iteration. We can then use linearmultigrid to solve these equations. A great strength of the multigrid idea, however, is that itcan be applied directlyto nonlinear problems. All we need is a suitable nonlinear relaxation method tosmooth theerrors,plus aprocedure for approximating corrections on coarsergrids.This direct approach is Brandt ’s Full Approximation Storage Algorithm (FAS). No nonlinear equations need be solved, except perhaps on the coarsest grid. To develop the nonlinear algorithm, suppose we have a relaxation procedure that can smooth theresidualvector aswedidinthelinearcase. Thenwecanseek asmooth correctionv hto solve (19.6.23): Lh(/tildewideuh+vh)=fh (19.6.24 ) Tofindvh, note that Lh(/tildewideuh+vh)−L h(/tildewideuh)=fh−L h(/tildewideuh) =−dh(19.6.25 ) The right-hand side is smooth after a few nonlinear relaxation sweeps. Thus we can transfer the left-hand side to a coarse grid: LH(uH)−L H(R/tildewideuh)=−Rdh (19.6.26 ) that is, we solve LH(uH)=LH(R/tildewideuh)−Rdh (19.6.27 ) on the coarse grid. (This is how nonzero right-hand sides appear.) Suppose the approximate solution is /tildewideuH. Then the coarse-grid correction is /tildewidevH=/tildewideuH−R/tildewideuh (19.6.28 ) 19.6MultigridMethodsforBoundaryValueProblems 875Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).and /tildewideunew h =/tildewideuh+P(/tildewideuH−R/tildewideuh)( 19.6.29 ) Notethat PR /negationslash =1ingeneral,so /tildewideunew h/negationslash=P/tildewideuH. Thisisakey point: Inequation (19.6.29) the interpolation error comes only from the correction, not from the full solution /tildewideuH. Equation (19.6.27) shows that one is solving for the full approximation uH, not just the error as in the linear algorithm. This is the origin of the name FAS. The FAS multigrid algorithm thus looks very similar to the linear multigrid algorithm. The only differences are that both the defect dhand the relaxed approximation uhhave to be restricted to the coarse grid, where now it is equation (19.6.27) that is solved by recursiveinvocation of the algorithm. However, instead of implementing the algorithm this way, wewillfirst describe the so-called dual viewpoint , which leads to a powerful alternative way of looking at the multigrid idea. The dual viewpoint considers the local truncation error ,d efined as τ≡L h(u)−fh (19.6.30 ) where uis the exact solution of the original continuum equation. If we rewrite this as Lh(u)=fh+τ (19.6.31 ) we see that τcan be regarded as the correction to fhso that the solution of the fine-grid equation will be the exact solution u. Now consider the relative truncation error τh, which is de fined on the H-grid relative to the h-grid: τh≡L H(Ruh)−R L h(uh)( 19.6.32 ) SinceLh(uh)=fh, this can be rewritten as LH(uH)=fH+τh (19.6.33 ) In other words, we can think of τhas the correction to fHthat makes the solution of the coarse-grid equation equal to the fine-grid solution. Of course we cannot compute τh, but we do have an approximation to it from using /tildewideuhin equation (19.6.32): τh/similarequal/tildewideτh≡L H(R/tildewideuh)−R L h(/tildewideuh)( 19.6.34 ) Replacing τhby/tildewideτhin equation (19.6.33) gives LH(uH)=LH(R/tildewideuh)−Rdh (19.6.35 ) which is just the coarse-grid equation (19.6.27)! Thus we see that there are two complementary viewpoints for the relation between coarse and fine grids: •Coarse grids are used to accelerate the convergence of the smooth components of thefine-grid residuals. •Fine grids are used to compute correction terms to the coarse-grid equations, yieldingfine-grid accuracy on the coarse grids. Onebene fitofthisnewviewpointisthatitallowsustoderiveanaturalstoppingcriterion for a multigrid iteration. Normally the criterion would be /bardbldh/bardbl≤/epsilon1 (19.6.36 ) and the question is how to choose /epsilon1. There is clearly no bene fit in iterating beyond the point when the remaining error is dominated by the local truncation error τ. The computable quantity is /tildewideτh. What is the relation between τand/tildewideτh? For the typical case of a second-order accurate differencing scheme, τ=Lh(u)−L h(uh)=h2τ2(x, y )+··· (19.6.37 ) 876 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Assume the solution satis fiesuh=u+h2u2(x, y )+···. Then, assuming Ris of high enough order that we can neglect its effect, equation (19.6.32) gives τh/similarequalL H(u+h2u2)−L h(u+h2u2) =LH(u)−L h(u)+h2[L/prime H(u2)−L/prime h(u2)] +··· =(H2−h2)τ2+O(h4)(19.6.38 ) For the usual case of H=2hwe therefore have τ/similarequal1 3τh/similarequal1 3/tildewideτh (19.6.39 ) The stopping criterion is thus equation (19.6.36) with /epsilon1=α/bardbl/tildewideτh/bardbl,α ∼1 3(19.6.40 ) We have one remaining task before implementing our nonlinear multigrid algorithm: choosing a nonlinear relaxation scheme. Once again, your first choice should probably be the nonlinear Gauss-Seidel scheme. If the discretized equation (19.6.23) is written withsome choice of ordering as L i(u1,...,u N)=fi,i =1,...,N (19.6.41 ) then the nonlinear Gauss-Seidel schemes solves Li(u1,...,u i−1,unew i,u i+1,...,u N)=fi (19.6.42 ) forunew i. Asusualnew u’sreplaceold u’sassoonastheyhavebeencomputed. Oftenequation (19.6.42)islinearin unew i,sincethenonlineartermsarediscretizedbymeansofitsneighbors. If this is not the case, we replace equation (19.6.42) by one step of a Newton iteration: unew i =uold i−Li(uold i)−fi ∂L i(uold i)/∂u i(19.6.43 ) For example, consider the simple nonlinear equation ∇2u+u2=ρ (19.6.44 ) In two-dimensional notation, we have L(ui,j)=( ui+1 ,j+ui−1,j+ui,j+1+ui,j−1−4ui,j)/h2+u2 i,j−ρi,j=0 (19.6.45 ) Since ∂L ∂u i,j=−4/h2+2ui,j (19.6.46 ) the Newton Gauss-Seidel iteration is unew i,j =ui,j−L(ui,j) −4/h2+2ui,j(19.6.47 ) Hereisaroutine mgfasthatsolvesequation(19.6.44)usingtheFullMultigridAlgorithm and the FAS scheme. Restriction and prolongation are done as in mglin. We have included theconvergencetestbasedonequation(19.6.40). Asuccessfulmultigridsolutionofaproblemshould aim to satisfy this condition with the maximum number of V-cycles, maxcyc, equal to 1 or 2. The routine mgfasuses the same subroutines copy,interp,maloc, and rstrct asmglin, but with a larger storage requirement MEMLENinmaloc(be sure to change the PARAMETER statement in that routine, as indicated by the commented line). 19.6MultigridMethodsforBoundaryValueProblems 877Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE mgfas(u,n,maxcyc) INTEGER maxcyc,n,NPRE,NPOST,NG,MEMLEN DOUBLE PRECISION u(n,n),ALPHA PARAMETER (NG=5,MEMLEN=17*2**(2*NG)/3+18*2**NG+10*NG-86/3)PARAMETER (NPRE=1,NPOST=1,ALPHA=.33d0) C USES anorm2,copy,interp,lop,maloc,matadd,matsub,relax2,rstrct,slvsm2 Full Multigrid Algorithm for FAS solution of nonlinear elliptic equation, here equation(19.6.44). On input u(1:n,1:n) contains the right-hand side ρ, while on output it re- turns the solution. The dimension nis related to the number of grid levels used in the solution, NGbelow, by n= 2** NG +1.maxcyc is the maximum number of V-cycles to be used at each level.Parameters: NGis the number of grid levels used; MEMLEN is the maximum amount of memory that can be allocated by calls to maloc ;NPRE andNPOST are the number of relaxation sweeps before and after the coarse-grid correction is computed; ALPHA relates the estimated truncation error to the norm of the residual. INTEGER j,jcycle,jj,jm1,jpost,jpre,mem,nf,ngrid,nn,irho(NG), * irhs(NG),itau(NG),itemp(NG),iu(NG),maloc DOUBLE PRECISION res,trerr,z,anorm2COMMON /memory/ z(MEMLEN),mem Storage for grid functions is allocated by maloc from array z. mem=0 nn=n/2+1 ngrid=NG-1irho(ngrid)=maloc(nn**2) Allocate storage for r.h.s. on grid NG−1, call rstrct(z(irho(ngrid)),u,nn) and fill it by restricting from the fine grid. 1 if (nn.gt.3) then Similarly allocate storage and fill r.h.s. on all coarse grids. nn=nn/2+1 ngrid=ngrid-1 irho(ngrid)=maloc(nn**2) call rstrct(z(irho(ngrid)),z(irho(ngrid+1)),nn) goto 1 endif nn=3iu(1)=maloc(nn**2)irhs(1)=maloc(nn**2) itau(1)=maloc(nn**2) itemp(1)=maloc(nn**2)call slvsm2(z(iu(1)),z(irho(1))) Initial solution on coarsest grid. ngrid=NG do 16j=2,ngrid Nested iteration loop. nn=2*nn-1iu(j)=maloc(nn**2) irhs(j)=maloc(nn**2) itau(j)=maloc(nn**2)itemp(j)=maloc(nn**2)call interp(z(iu(j)),z(iu(j-1)),nn) Interpolate from coarse grid to next finer grid. if (j.ne.ngrid) then call copy(z(irhs(j)),z(irho(j)),nn) Set up r.h.s. else call copy(z(irhs(j)),u,nn) endifdo 15jcycle=1,maxcyc V-cycle loop. nf=nn do12jj=j,2,-1 Downward stoke of the V. do11jpre=1,NPRE Pre-smoothing. call relax2(z(iu(jj)),z(irhs(jj)),nf) enddo 11 call lop(z(itemp(jj)),z(iu(jj)),nf) Lh(/tildewideuh). nf=nf/2+1jm1=jj-1 call rstrct(z(itemp(jm1)),z(itemp(jj)),nf) RL h(/tildewideuh). call rstrct(z(iu(jm1)),z(iu(jj)),nf) R/tildewideuh. call lop(z(itau(jm1)),z(iu(jm1)),nf) LH(R/tildewideuh)stored temporarily in /tildewideτh. call matsub(z(itau(jm1)),z(itemp(jm1)),z(itau(jm1)),nf) Form /tildewideτh. if(jj.eq.j)trerr=ALPHA*anorm2(z(itau(jm1)),nf) Estimate truncation error τ. 878 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).call rstrct(z(irhs(jm1)),z(irhs(jj)),nf) fH. call matadd(z(irhs(jm1)),z(itau(jm1)),z(irhs(jm1)),nf) fH+/tildewideτh. enddo 12 call slvsm2(z(iu(1)),z(irhs(1))) Bottom of V: Solve on coarsest grid. nf=3 do14jj=2,j Upward stroke of V. jm1=jj-1call rstrct(z(itemp(jm1)),z(iu(jj)),nf) R/tildewideu h. call matsub(z(iu(jm1)),z(itemp(jm1)),z(itemp(jm1)),nf) /tildewideuH−R/tildewideuh. nf=2*nf-1 call interp(z(itau(jj)),z(itemp(jm1)),nf) P(/tildewideuH−R/tildewideuh)stored in /tildewideτh. call matadd(z(iu(jj)),z(itau(jj)),z(iu(jj)),nf) Form /tildewideunew h. do13jpost=1,NPOST Post-smoothing. call relax2(z(iu(jj)),z(irhs(jj)),nf) enddo 13 enddo 14 call lop(z(itemp(j)),z(iu(j)),nf) Form residual /bardbldh/bardbl. call matsub(z(itemp(j)),z(irhs(j)),z(itemp(j)),nf)res=anorm2(z(itemp(j)),nf)if(res.lt.trerr)goto 2 No more V-cycles needed if residual small enough. enddo 15 2 continue enddo 16 call copy(u,z(iu(ngrid)),n) Return solution in u. returnEND SUBROUTINE relax2(u,rhs,n) INTEGER n DOUBLE PRECISION rhs(n,n),u(n,n) Red-black Gauss-Seidel relaxation for equation (19.6.44). The current value of the solution u(1:n,1:n) is updated, using the right-hand side function rhs(1:n,1:n) . INTEGER i,ipass,isw,j,jsw DOUBLE PRECISION foh2,h,h2i,res h=1.d0/(n-1)h2i=1.d0/(h*h)foh2=-4.d0*h2i jsw=1 do 13ipass=1,2 Red and black sweeps. isw=jsw do12j=2,n-1 do11i=isw+1,n-1,2 res=h2i*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)- * 4.d0*u(i,j))+u(i,j)**2-rhs(i,j) u(i,j)=u(i,j)-res/(foh2+2.d0*u(i,j)) Newton Gauss-Seidel formula. enddo 11 isw=3-isw enddo 12 jsw=3-jsw enddo 13 return END SUBROUTINE slvsm2(u,rhs) DOUBLE PRECISION rhs(3,3),u(3,3) C USES fill0 Solution of equation (19.6.44) on the coarsest grid, where h=1 2. The right-hand side is input in rhs(1:3,1:3) and the solution is returned in u(1:3,1:3) . DOUBLE PRECISION disc,fact,h call fill0(u,3) h=.5d0fact=2.d0/h**2 disc=sqrt(fact**2+rhs(2,2)) 19.6MultigridMethodsforBoundaryValueProblems 879Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).u(2,2)=-rhs(2,2)/(fact+disc) return END SUBROUTINE lop(out,u,n) INTEGER n DOUBLE PRECISION out(n,n),u(n,n) Given u(1:n,1:n) , returns Lh(/tildewideuh)for equation (19.6.44) in out(1:n,1:n) . INTEGER i,j DOUBLE PRECISION h,h2i h=1.d0/(n-1)h2i=1.d0/(h*h) do 12j=2,n-1 Interior points. do11i=2,n-1 out(i,j)=h2i*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)- * 4.d0*u(i,j))+u(i,j)**2 enddo 11 enddo 12 do13i=1,n Boundary points. out(i,1)=0.d0 out(i,n)=0.d0out(1,i)=0.d0out(n,i)=0.d0 enddo 13 return END SUBROUTINE matadd(a,b,c,n) INTEGER nDOUBLE PRECISION a(n,n),b(n,n),c(n,n) Adds a(1:n,1:n) tob(1:n,1:n) and returns result in c(1:n,1:n) . INTEGER i,jdo 12j=1,n do11i=1,n c(i,j)=a(i,j)+b(i,j) enddo 11 enddo 12 returnEND SUBROUTINE matsub(a,b,c,n) INTEGER nDOUBLE PRECISION a(n,n),b(n,n),c(n,n) Subtracts b(1:n,1:n) from a(1:n,1:n) and returns result in c(1:n,1:n) . INTEGER i,j do12j=1,n do11i=1,n c(i,j)=a(i,j)-b(i,j) enddo 11 enddo 12 return END DOUBLE PRECISION FUNCTION anorm2(a,n) INTEGER n DOUBLE PRECISION a(n,n) Returns the Euclidean norm of the matrix a(1:n,1:n) . INTEGER i,j DOUBLE PRECISION sum sum=0.d0do 12j=1,n do11i=1,n 880 Chapter19. PartialDifferentialEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).sum=sum+a(i,j)**2 enddo 11 enddo 12 anorm2=sqrt(sum)/nreturn END CITED REFERENCES AND FURTHER READING: Brandt, A. 1977, Mathematics of Computation , vol. 31, pp. 333 –390. [1] Hackbusch, W. 1985, Multi-Grid Methods and Applications (New York: Springer-Verlag). [2] Stuben, K., and Trottenberg, U. 1982, in Multigrid Methods , W. Hackbusch and U. Trottenberg, eds.(SpringerLectureNotesinMathematics No.960)(NewYork:Springer-Verlag),pp.1 – 176. [3] Brandt,A.1982,in MultigridMethods ,W.HackbuschandU.Trottenberg,eds.(SpringerLecture Notes in Mathematics No. 960) (New York: Springer-Verlag). [4] Baker, L. 1991, More C Tools for Scientists and Engineers (New York: McGraw-Hill). Briggs, W.L. 1987, A Multigrid Tutorial (Philadelphia: S.I.A.M.). Jespersen, D. 1984, Multigrid Methods for Partial Differential Equations (Washington: Mathe- matical Association of America). McCormick,S.F.(ed.)1988, MultigridMethods:Theory,Applications,andSupercomputing (New York: Marcel Dekker). Hackbusch, W., andTrottenberg, U. (eds.) 1991, Multigrid Methods III (Boston: Birkhauser). Wesseling, P. 1992, An Introduction to Multigrid Methods (New York: Wiley).