A Multigrid Tutorial, Second Edition

The front matter includes the title page, copyright page, dedicaton, TOC, preface to second edition and preface to first edition

Multigrid methods were originally applied to simple boundary value problems that arise in many physical applications. For simplicity and for historical reasons, these problems provide a natural introduction to multigrid methods. As an example, consider the two-point boundary value problem that describes the steady-state temperature distribution in a long uniform rod. It is given by the second-order boundary value problem − u″ (x)+σu (x) = ƒ (x), 0<x<1,σ≥0, 1.1 u (0)=u (1) = 0. 1.2

While this problem can be handled analytically, our present aim is to consider numerical methods. Many such approaches are possible, the simplest of which is a finite difference method (finite element formulations will be considered in Chapter 10). The domain of the problem {x:0≤x≤1} is partitioned into n subintervals by introducing the grid points xj =jh , where h=1/n is the constant width of the subintervals. This establishes the grid shown in Fig. 1.1, which we denote Ωh .

At each of the n − 1 interior grid points, the original differential equation (1.1) is replaced by a second-order finite difference approximation. In making this replacement, we also introduce υj as an approximation to the exact solution u ( xj ) . This approximate solution may now be represented by a vector v= ( υ1 ,… , υn−1 )T , whose components satisfy the n − 1 linear equations − υj−1 +2 υj − υj+1 h2 +σ υj = ƒ ( xj ), 1≤j≤n−1, υ0 = υn = 0. 1.3

We now consider how model problems (1.3) and (1.5) might be treated using conventional iterative or relaxation methods. We first establish the notation for this and all remaining chapters. Let Au=f denote a system of linear equations such as (1.3) or (1.6). We always use u to denote the exact solution of this system and v to denote an approximation to the exact solution, perhaps generated by some iterative method. Bold symbols, such as u and v, represent vectors, while the jth components of these vectors are denoted by uj and υj . In later chapters, we need to associate u and v with a particular grid, say Ωh . In this case, the notation uh and vh is used.

Suppose that the system Au=f has a unique solution and that v is a computed approximation to u. There are two important measures of v as an approximation to u. One is the error (or algebraic error) and is given simply by e=u−v. The error is also a vector and its magnitude may be measured by any of the standard vector norms. The most commonly used norms for this purpose are the maximum (or infinity) norm and the Euclidean or 2-norm, defined, respectively, by ‖e‖∞ = max 1≤j≤n | ej | and ‖e‖2 = ∑ j=1 n e j 2 1/2 .

Unfortunately, the error is just as inaccessible as the exact solution itself. However, a computable measure of how well v approximates u is the residual, given by r=f−Av. The residual is simply the amount by which the approximation v fails to satisfy the original problem Au=f . It is also a vector and its size may be measured by the same norm used for the error. By the uniqueness of the solution, r=0 if and only if e=0 . However, it may not be true that when r is small in norm, e is also small in norm.

Through analysis and experimentation, we have examined some of the basic iterative methods. Our discoveries have formed the beginnings of what we might call a spectral (or Fourier mode) picture of relaxation schemes. As we proceed, more essential details of this picture will become clear. So far we have established that many standard iterative methods possess the smoothing property. This property makes these methods very effective at eliminating the high-frequency or oscillatory components of the error, while leaving the low-frequency or smooth components relatively unchanged. The immediate issue is whether these methods can be modified in some way to make them effective on all error components.

One way to improve a relaxation scheme, at least in its early stages, is to use a good initial guess. A well-known technique for obtaining an improved initial guess is to perform some preliminary iterations on a coarse grid. Relaxation on a coarse grid is less expensive because there are fewer unknowns to be updated. Also, because the convergence factor behaves like 1−O ( h2 ) , the coarse grid will have a marginally improved convergence rate. This line of reasoning at least suggests that coarse grids might be worth considering.

With the coarse grid idea in mind, we can think more carefully about its implications. Recall that most basic relaxation schemes suffer in the presence of smooth components of the error. Assume that a particular relaxation scheme has been applied until only smooth error components remain. We now ask what these smooth components look like on a coarser grid. Figure 3.1 shows the answer. A smooth wave with k=4 on a grid Ωh with n=12 points has been projected directly to the grid Ω2h with n=6 points. On this coarse grid, the original wave still has a wavenumber of k=4 . We see that a smooth wave on Ωh looks more oscillatory on Ω2h .

To be more precise, note that the grid points of the coarse grid Ω2h are the even-numbered grid points of the fine grid Ωh . Consider the kth mode on the fine grid evaluated at the even-numbered grid points. if 1≤k< n 2 , its components may be written as w k,2j h = sin 2jkπ n = sin jkπ n/2 = w k,j 2h , 1≤k < n 2 . Notice that superscripts have been used to indicate the grids on which the vectors are defined.

The preceding chapter was devoted to the development of several multigrid schemes. We now turn to the practical issues of writing multigrid programs and determining whether they work. This will lead us to issues such as data structures, complexity, predictive tools, diagnostic tools, and performance.

Complexity

Writing multigrid programs can be both fun and challenging. The experience of many practitioners suggests that such programs should be highly modular. This allows them to evolve from simple relaxation programs and makes them much easier to check and debug. Also, the various components of the program (for example, relaxation, interpolation, and restriction subroutines) can be replaced individually.

Choosing a manageable data structure for a multigrid program is essential. Modern programming languages are replete with devices that make data management easy. For example, in most languages, one can declare a structure that groups together all the associated information for each grid level. In a structured language, for instance, a V-cycle could be written along the lines of the following pseudocode:

declare structure:

 grid = { double Ddim_array f %% the right hand side

 double Ddim_array v %% the current approximation }

declare Grid: array of structure grid

for j = 0 to coarsest - 1

 Grid[j].v <- relax(Grid[j].v, Grid[j].f, num_sweeps_down)

 Grid[j+1].f <- restrict(Grid[j].f - apply_operator(Grid[j].v))

endfor

Grid[coarsest].v = direct_solve(Grid[coarsest].v, Grid[coarsest].f)

for j = coarsest-1 to 1

 Grid[j].v <- Grid[j].v + interpolate(Grid[j+1].v)

 Grid[j].v <- relax(Grid[j].v, Grid[j].f, num_sweeps_down)

endfor

In the previous chapter, we considered some practical questions concerning the implementation, cost, and performance of multigrid methods. The arguments and experiments of that chapter offer good reason to believe that multigrid methods can be extremely effective. Now we must confront some questions on the theoretical side. The goal of this chapter is to present multigrid in a more formal setting and offer an explanation of why these methods work so well. In the first part of this chapter, we sketch the ideas that underlie the convergence theory of multigrid. In the second section, we present what might be called the subspace picture of multigrid. While the terrain in this chapter may seem a bit more rugged than in previous chapters, the reward is an understanding of why multigrid methods are so remarkably effective.

Variational Properties

We first return to a question left unanswered in previous chapters. In expressing the coarse-grid problem, we wrote A2h u2h = f2h and said that A2h is the Ω2h version of the original operator Ah . Our first goal is to define the coarse-grid operator A2h precisely.

The argument that follows assumes we are working with the model problem, − u″ (x) =ƒ (x) , and the corresponding discrete operator Ah . We adopt the notation that Ωph represents not only the grid with grid spacing ph, but also the space of vectors defined on that grid. In keeping with our former notation, we let vh be a computed approximation to the exact solution uh . For the purpose of this argument, assume that the error in this approximation, eh = uh − vh , lies entirely in the range of interpolation, which will be denoted ℛ ( I 2h h ) . This means that for some vector u2h ∈ Ω2h , eh = I 2h h u2h . Therefore, the residual equation on Ωh may be written Ah eh = Ah I 2h h u2h = rh . 5.1 In this equation, Ah acts on a vector that lies entirely in the range of interpolation. How does Ah act on ℛ ( I 2h h ) ? Figure 5.1 gives the answer. An arbitrary vector u2h ∈ Ω2h is shown in Fig. 5.1(a); I 2h h u2h , which is the interpolation of u2h up to Ωh , is shown in Fig. 5.1(b); and the effect of Ah acting pointwise on I 2h h u2h is shown in Fig. 5.1(c).

The preceding chapters presented multigrid methods applied to linear problems. Because of their iterative nature, multigrid ideas should be effective on nonlinear problems, where some form of iteration is usually imperative. Just as multigrid methods were developed to improve classical linear relaxation methods, we will see that they can also guide us in the use of nonlinear methods.

We begin by clarifying the most significant formal difference between linear and nonlinear systems of equations. Consider a system of nonlinear algebraic equations, A (u) =f , where u,f ∈ Rn . (The notation A(u), rather than Au, signifies that the operator is nonlinear.) Suppose that v is an approximation to the exact solution u. It is possible to define the error and the residual for this problem, much as we did earlier: the error is simply e=u−v and the residual is r=f−A (v) . If we now subtract the original equation, A (u) =f , from the definition of the residual, we find that A (u) −A (v) =r. 6.1 We are one step from the residual equation that we derived for linear systems in Chapter 2. However, because A is nonlinear, even though u−v=e , we cannot conclude that A (u) −A (v) =A (e) . This means that we no longer have a simple linear residual equation, which necessitates changes in the methods that we devise for nonlinear problems. We must now use (6.1) as the residual equation.

It makes sense to review quickly a classical relaxation method for nonlinear systems; it will have a part to play in multigrid methods. Working with the system of nonlinear algebraic equations A (u) =f , one of the most frequently used methods is nonlinear Gauss-Seidel relaxation [16]. Suppose the jth equation of the system can be solved for the jth variable in terms of the other variables ( u1 , u2 ,…, uj−1 , uj+1 ,…, un ) . We write the result in the abstract form uj = Mj ( u1 , u2 ,…, uj−1 , uj+1 ,…, un ), 1≤j≤n. Just as with linear systems, Gauss—Seidel relaxation takes the form υj ← Mj ( υ1 , υ2 ,…, υj−1 , υj+1 ,…, υn ) ,1≤j≤n, where the current value of each component is used on the right side.

Our model problem for highlighting basic multigrid concepts has been Poisson's equation in one and two dimensions. To see how multigrid can be applied to other problems, we now consider three variations on our model problem, each, we hope, of fairly general interest: Neumann boundary conditions, anisotropic equations, and variable-coefficient/mesh problems. In each case, we present a typical, but simplified, problem that we treat with a specially tailored multigrid technique. The method we propose is not necessarily the most robust, and we ignore complicating issues for the sake of clarity. The purpose here is to illustrate the application of multigrid to new problems and to gain an entry-level understanding of the broad scope of multigrid techniques.

Neumann Boundary Conditions

Dirichlet boundary conditions specify values of the unknown function on the boundary of the domain. In most cases, these so-called essential boundary conditions are easily treated in the discretization and in the multigrid solver by simply eliminating the corresponding unknowns in the equations. On the other hand, Neumann conditions involve derivatives of the unknown function. This relatively small change in the boundary conditions can be treated in a variety of ways; it also introduces a few important subtleties.

Consider the following two-point boundary value problem with homogeneous Neumann conditions at both end points (the inhomogeneous case is left for Exercise 8): − u″ (x) =ƒ (x) ,0<x<1, u′ (0) = u′ (1) =0. 7.1 With n an odd positive integer, denote the points of a uniform grid on the interval [0,1] by xj =jh , for 0≤j≤n+1 , where h= 1 n+1 . A common way to discretize this problem involves the ghost points x−1 =− 1 h and xn+2 =1+ 1 h . These points are outside the problem domain and are used only temporarily to produce approximations at the boundary points.

With this extended grid, we can use central differences for the differential equation at the interior points and boundary points.

A natural question arises: Can we apply multigrid techniques when there is no grid? Suppose we have relationships among the unknowns that are similar to those in the model problem, but the physical locations of the unknowns are themselves unknown (or immaterial). Can we hope to apply the tools we have developed? A related question is: Can we apply multigrid in the case where grid locations are known but may be highly unstructured or irregular, making the selection of a coarse grid problematic? These are the problems that are addressed by a technique known as algebraic multigrid, or AMG [6]. We describe AMG using many of the concepts and principles developed in [18]. For a theoretical foundation of general algebraic methods, of which AMG is a part, see [5].

For any multigrid algorithm, the same fundamental components are required. There must be a sequence of grids, intergrid transfer operators, a relaxation (smoothing) operator, coarse-grid versions of the fine-grid operator, and a solver for the coarsest grid.

Let us begin by deciding what we mean by a grid. Throughout this chapter, we look to standard multigrid (which we refer to as the geometric case) to guide us in defining AMG components. In the geometric case, the unknown variables ui are defined at known spatial locations (grid points) on a fine grid. We then select a subset of these locations as a coarse grid. As a consequence, a subset of the variables ui is used to represent the solution on the coarse grid. For AMG, by analogy, we seek a subset of the variables ui to serve as the coarse-grid unknowns. A useful point of view, then, is to identify the grid points with the indices of the unknown quantities.

Numerical problems often exhibit special features in small local regions that require resolution and accuracy well beyond what is required in the rest of the domain. In numerical weather modeling, isolated phenomena such as tornados or storm fronts may demand substantially enhanced accuracy. The numerical simulation of the flight of an aircraft may require especially accurate approximations around the fuselage or wings. In these and many other cases, it is wasteful to let local accuracy requirements dictate the global discretization and solution process. The goal of this chapter is to understand how to treat local demands in a multilevel context without overburdening the overall computation.

There is a wide variety of methods that effectively treat local demands. Here we consider the so-called fast adaptive composite grid method (FAC) [14]. Its distinctive features are that it always works with uniform grids and subgrids, and it is in tune with the variational theme that appears throughout this book. We begin with a one-dimensional example.

Consider the two-point boundary value problem − u″ (x) =ƒ (x), 0<x<1,u (0)=u (1) =0. To keep matters simple, assume that a local fine grid with mesh size h= 1 8 is needed on the interval ( 1 2 , 1) to resolve some special feature near x= 3 4 , but that a mesh size of 2h= 1 4 is deemed adequate for the rest of the domain. This need might arise, for example, in the presence of a source term ƒ that is nonzero only near x= 3 4 (Fig. 9.1).

Imagine first that we solve the problem on a global grid, Ωh , with mesh size h= 1 8 , using a two-grid method based on linear interpolation, full weighting, and the variational properties as described in Chapter 5. We write the discrete problem as Ah uh = fh . 9.1 In component form, we have − u i−1 h + 2u i h − u i+1 h h2 = f i h ,1≤i≤n, u 0 h = u n+1 h =0, where n=7 for this case.

Until now, our development of multigrid methods has relied on finite differences. However, there are several other discretization methods, all of which have their place in computation. We now turn to the finite element method. This approach may seem more abstract and less direct at the outset, but its use of variational properties provides a powerful framework that is well suited to developing and analyzing multigrid methods. It is impossible to cover finite elements completely in this short chapter, so we focus instead on surveying the critical ideas, leaving the important details to the exercises and further study.

Said simply, finite difference methods replace the problem domain by a grid and produce a vector whose components are approximations to the solution at the grid points. On the other hand, finite element methods partition the problem domain into subregions and produce a simple function in each subregion that approximates the solution. Finite element methods can be applied to two major categories of problems: differential equations and functional minimization. When the differential equation is defined in terms of a self-adjoint linear operator and the functional is defined appropriately, these two categories coincide. To exploit this fortunate duality between differential equations and functional minimization, we focus on the self-adjoint case.

Let Ω be a suitably nice bounded domain in the plane with boundary ∂Ω. Consider the following Poisson equation with homogeneous Dirichlet boundary conditions: Lu≡− uxx − uyy =ƒin Ω,u=0on ∂Ω,

We first need to introduce some new tools and notation. Most of our time will be spent in the space L2 (Ω) of square integrable functions on Ω (functions that satisfy ∫Ω u2 dΩ<∞ equipped with the L2 (Ω) inner product (⋅, ⋅) defined by (u,v)= ∫Ω uv dΩ. Because we will be taking partial derivatives, it is not enough for the functions we consider just to be L2 (Ω) integrable. We need to know that they are differentiable.

The back matter includes bibliography and index