A Complex Analysis Based Derivation of Multigrid Smoothing Factors for Lexicographic Gauss-Seidel by Laird Robert Hocking B.Sc., The University of British Columbia, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Computer Science) The University Of British Columbia (Vancouver) August 2011 c Laird Robert Hocking, 2011 Abstract This thesis aims to present a unified framework for deriving analytical formulas for multigrid smoothing factors in arbitrary dimensions, under certain simplifying assumptions. To derive these expressions we rely on complex analysis and geometric considerations, using the maximum modulus principle and M¨obius transformations. We restrict our attention to pointwise and block lexicographic Gauss-Seidel smoothers on a d-dimensional uniform mesh, where the computational molecule of the associated discrete operator forms a 2d + 1 point star. In the pointwise case the effect of a relaxation parameter, as well as different choices of mesh ratio, are analyzed. The results apply to any number of spatial dimensions, and are applicable to high-dimensional versions of a few common model problems with constant coefficients, including the Poisson and anisotropic diffusion equations as well as the convection-diffusion equation. We show that in most cases our formulas, exact under the simplifying assumptions of Local Fourier Analysis, form tight upper bounds for the asymptotic convergence of geometric multigrid in practice. We also show that there are asymmetric cases where lexicographic Gauss-Seidel smoothing outperforms red-black Gauss-Seidel smoothing; this occurs for certain model convection-diffusion equations with high mesh Reynolds numbers. ii Preface This thesis describes the results of a refereed journal paper, L. Robert Hocking and Chen Greif. Closed-Form Multigrid Smoothing Factors for Lexicographic Gauss-Seidel. IMA Journal of Numerical Analysis, accepted for publication, August 2011. A body of unpublished work is described as well. The results of the paper are described in Chapters 3 and 5, while Chapters 2 and 4 represent new additions. I have taken a central role in all aspects of the paper. I am the first author, with collaboration by my supervisor Chen Greif. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Thesis Contributions . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Overview of Multigrid . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 3 4 2.1 Quantitative Convergence Analysis . . . . . . . . . . . . . . . . . 2.2 A Closer Look at the Interplay Between Errors on the Fine and 11 Coarse Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Smoothing Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1 Pointwise Relaxation . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 SOR Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Block Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Alternative Mesh-Ratios . . . . . . . . . . . . . . . . . . . . . . . . . 33 iv 5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.1 The d-Dimensional Poisson Problem . . . . . . . . . . . . . . . . 40 5.2 The Anisotropic Steady-State Diffusion Problem . . . . . . . . . 41 5.3 Alternative Coarsening for the 2D Poisson Problem . . . . . . . . 42 5.4 The Convection-Diffusion Problem . . . . . . . . . . . . . . . . . 44 5.5 The Time-dependent Diffusion Problem . . . . . . . . . . . . . . 47 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6 v List of Tables Table 4.1 Relative work of optimal coarsening and corresponding values of βopt , for various d-dimensional Poisson problems. . . . . . . Table 5.1 36 Smoothing factors of GS-LEX with point, line, and plane relaxation for the Poisson problem in 1D, 2D, and 3D. Moving from point to line or from line to plane relaxation is equivalent to reducing the dimension by one. . . . . . . . . . . . . . . . . Table 5.2 40 Comparison of predicted smoothing factors of GS-LEX with measured convergence rates of multigrid V (1, 1)-cycles, for the PDE −εuxx − uyy = f . The finest grid is 1023 × 1023, while the coarsest is 1 × 1. Various values of ε are considered. . . . . . . Table 5.3 42 Predicted smoothing factors smoothing vs. measured convergence rates of multigrid with H = kh coarsening and GS-LEX smoothing, applied to the 2D Poisson problem discretized using centered differences, for various values of k. V (1, 0)-cycles are used, with a finest grid of size kNk − 1 in each direction, and a coarsest grid consisting of a single point. Nk = logk (999) is chosen to make the total number of grid points approximately one million. Also included is the total work as measured by the running time of the program divided by the size of the problem - this is given both in microseconds per grid point and as a ratio relative to the reference case H = 2h. This latter form is compared against the theoretical relative work (4.2). . . . . . . vi 43 Table 5.4 Smoothing factors of pointwise GS-LEX applied to two common discretizations of the convection-diffusion equation with all mesh Reynolds numbers equal to γ, for d = 1, 2, 3. . . . . . Table 5.5 45 Comparison of the smoothing factor µ of GS-LEX with the measured asymptotic convergence rate of multigrid applied to the PDE −uxx − uyy + σ ux + σ uy = f , discretized uniformly. The problem is discretized using upwinding on a 1023 × 1023 grid. The Galerkin coarse grid operator is used. W (1, 0)-cycles are used, and various values of γ = σ h/2, where h is the grid spacing on the finest mesh, are considered. ‘NC’ stands for no convergence. ‘nL’ (with n = 2, 3, 4) signifies the number of levels. 45 Table 5.6 Comparison of the smoothing factor µ of GS-LEX with measured convergence rates of multigrid V (1, 0) and W (1, 0)-cycles, for the PDE −uxx − uyy + σ (ux + uy ) = f , discretized uniformly. The finest grid is 1023 × 1023, while the coarsest is 1 × 1. The coarse grid operator is obtained by direct discretization of the PDE on the coarse mesh. Various values of γ = σ h/2 are considered, where h is the grid spacing of the finest mesh. ‘2L’ stands for a two-level scheme, whereas V (1, 0) and W (1, 0) signify V and W -cycles, respectively, with one pre-smoothing sweep and no post-smoothing. . . . . . . . . . . . . . . . . . . Table 5.7 46 Comparison of the smoothing factor µ of GS-LEX with the measured asymptotic convergence rate of two-level and multilevel multigrid applied to the PDE −uxx − uyy + σ (ux + uy ) = f , discretized uniformly. The problem is discretized using centered differences on a 1023 × 1023 grid. The coarse grid operator is obtained by direct discretization of the PDE on the coarse mesh. W (1, 0)-cycles are used, and various values of γ = σ h/2, where h is the grid spacing on the finest mesh, are considered. ‘NC’ stands for no convergence. . . . . . . . . . . . . . . . . vii 47 List of Figures Figure 1.1 The computational molecule associated with the operator L h (d = 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.1 3 Top: a ‘smooth’ image is shrunk by a factor of 4 in each direction, and then blown back up. Bottom: the same is done for a ‘rough’ image. Reduction and subsequent enlargement were done using Photoshop’s bilinear image scaling option. . . . . Figure 2.2 Radii Petrovich Fedorenko (1930-2009) is credited with inventing the multigrid algorithm. . . . . . . . . . . . . . . . . Figure 2.3 6 7 A single iteration of two-level multigrid applied to the twodimensional Poisson equation on the square (0, 1)2 , discretized on a 23 × 23 uniform grid using centered differences. The initial error, containing both high and low frequencies, is shown on the left. In the middle we see the error after the smoother S (in this case three iterations of lexicographic Gauss-Seidel) is applied - the high frequency components of the error have been significantly dampened. On the right we see the new error after the coarse grid correction K has been applied - all frequencies are now highly reduced. . . . . . . . . . . . . . . . . . . . . Figure 2.4 8 The sequence of grids visited by one iteration of four level multigrid, for two different choices of γ. (a) γ = 1: a V-cycle (b) γ = 2: a W-cycle. . . . . . . . . . . . . . . . . . . . . . . viii 10 Figure 2.5 ||e(k) ||2 vs. k for two-level multigrid with lexicographic SOR smoothing applied to anisotropic diffusion problem 1 100 uxx + uyy + uzz = 0, discretized using centered differences. (a) ω = 1. (GS-LEX) (b) ω = 1.63 (SOR-LEX). Note that for ω = 1, the norm of the error goes down monotonically, whereas for ω = 1.63 there is a initial increase before the asymptotic regime is entered. Convergence is faster overall for ω = 1.63, despite the erratic behavior of the error initially. . . . . . . . . Figure 2.6 13 The distinct Fourier modes ϕ1 (k) = sin(kπ/20) and ϕ19 (k) = sin(19kπ/20) on the fine grid Gh = {k/20}20 k=0 restrict to the same mode ϕ1 (K) = sin(Kπ/10) on the coarse grid GH = {K/10}10 K=0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3.1 17 An illustration of Example 3.1.3: (a) If a double pendulum with arms of length 1 6 and 1 3 is allowed to swing freely, it re- mains inside a disk of radius 0.5. (b) If one arm is constrained to lie in the left half plane, the pendulum swings in the heartshaped region shown. In both cases the pendulum is unable to reach a set of points surrounding the origin. . . . . . . . . . . 23 Figure 3.2 The heart-shaped set D is the union of two disks and a half disk. 25 Figure 3.3 Left: Smoothing factor of SOR-LEX with relaxation parameter ω, as a function of ω, for the standard seven-point centered differences discretization of the 3D Poisson problem. Right: Smoothing factors of GS-LEX and optimal SOR-LEX applied to the anisotropic diffusion problem −εuxx − uyy − uzz = f , as a function of ε, assuming seven-point centered differences discretization. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.1 29 (a) An example pacman set with c = 0.5 and β = 1/4. (b) The corresponding set Dβ when cm = cr = 0.25. This example corresponds the standard centered differences discretization of the two dimensional Poisson problem, with the coarsening strategy H = 4h. . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 34 Figure 4.2 The relative work Wβ /W0.5 vs β for the d-dimensional Poisson problem with d = 2, 3, 4. . . . . . . . . . . . . . . . . . . . . Figure 4.3 36 A fine and coarse grid with mesh ratio h/H = 2/3. The stencil of the projection operator varies from point to point, depending on whether the coarse grid lines up with fine grid. . . . . . . . Figure 4.4 The relative work Wβ /W0.5 for the d-dimensional Poisson problem converges to the limit curve shown as d → ∞. . . . . . . . Figure 5.1 37 38 Theoretical smoothing factors of GS-LEX vs. experimental asymptotic convergence rates of two-level multigrid with GS-LEX and GS-RB smoothers, applied to the steady-state diffusion equation −uxx − uyy − εuzz = f . Pointwise relaxation appears on the left, while x-oriented line relaxation is depicted on the right. The discretization is on a 23×23×23 grid with centered differences. . . . . . . . . . . . . . . . . . . . . . . . . . . . x 41 Acknowledgments I would like to thank my supervisor Chen Greif for his ongoing guidance, support, and friendship. I would also like to acknowledge the group of rotations and scalings in two di√ mensional Euclidean space, for possessing a structure such that −1 and complex analysis may exist1 , thereby making the present work possible. 1 It is perhaps unfair to give sole credit to the aforementioned group of matrices. I should also give thanks to the structure of polynomials, in particular with respect to a certain equivalence relation. Perhaps more fundamentally I should thank R2 , for admitting a field, or maybe even reason itself, for making this inevitable. xi Chapter 1 Introduction Multigrid algorithms, a class of iterative methods for the numerical solution of elliptic partial differential equations, are based on the alternating application of two numerical procedures having complementary effects on the error. The first, called a smoother, is typically a simple relaxation scheme such as damped Jacobi or Gauss-Seidel. While ineffective at making the error small, smoothers are effective at making it smooth (in a geometric sense, relative to the underlying mesh). The second, called coarse grid correction, projects the error down to a coarser grid (with little loss in information, because the error is smooth), solves for it there, and interpolates the result onto the original grid. The result is an approximate error which may be used to update the current guess in such a way that the new error is small, but not necessarily smooth. The smoother may then be applied again, and the process continued. A detailed overview of the method is provided in Chapter 2. For certain elliptic problems (such as the Poisson problem) multigrid is the most efficient known method to date, the total work being linear in the number of unknowns with small hidden constants [19]. Unfortunately, the performance of the method is quite sensitive to the nature of the problem at hand. For example, if multigrid with lexicographic Gauss-Seidel (GS-LEX) smoothing is applied to the anisotropic diffusion problem εuxx + uyy = f 1 discretized using centered differences, convergence becomes arbitrarily slow as ε → 0 (see Chapter 5 Section 5.2). In this case, the issue is the smoother - although effective for ε ≈ 1, GS-LEX fails to provide adequate smoothing for ε 1. The problem may be resolved by, for example, replacing the point GS-LEX smoother with a block GS-LEX procedure that relaxes unknowns in the y-direction simultaneously. In general, one must tailor the smoother to the problem at hand, and finding an effective smoother is a major component in the design of an efficient multigrid algorithm. An important tool in this respect is Local Fourier Analysis (LFA), which provides a quantitative estimate (called the smoothing factor and denoted by µ) of the asymptotic error reduction factor per iteration of multigrid with a given smoother applied to a given problem [3]. LFA isolates the effect of the smoother by making idealizing assumptions regarding the coarse grid correction and doing away with the issue of boundary conditions by assuming an infinite grid - it therefore represents the efficiency a multigrid method can attain if the coarse grid correct works perfectly1 and if boundary conditions are handled suitably (in some cases by additional relaxations near the boundary) [13, p. 100]. Although smoothing factors may be computed numerically by straightforward techniques, closed form formulas are desirable if possible. This is particularly true when one wishes to understand how the smoothing factor varies as a function of one or more parameters in the PDE of interest (for example the mesh Reynolds numbers of a discrete convection-diffusion problem). An example of the success of the analytical approach in this regard is the 1995 paper of Irad Yavneh [17], where closed form expressions for pointwise and block Gauss Seidel with red-black ordering (GS-RB) are derived for the family of problems c1 ux1 x1 + c2 ux2 x2 + . . . cd uxd xd = f , discretized using centered differences, and normalized so that ∑dk=1 ck = 1. Here, among other observations, it is found that in the pointwise case µ depends only on cmin ≡ min(c1 , c2 , . . . , cd ).2 1 In Chapter 5 Section 5.4 we will see an example where the assumptions of LFA regarding the coarse grid correction are not realistic. 2 In Chapter 3, Theorem 3.1.5 I will show that the same is true of GS-LEX. 2 In contrast to GS-RB, the literature is relatively sparse with regards to closed form formulas for smoothing factors of GS-LEX, and seems to be limited to the case where d ≤ 2. In particular, while the smoothing factor of GS-LEX applied to the standard 5-point discretization of the 2D Poisson problem has been known to be µ2D = 0.5 since 1977 [3], an analogous expression for the 3D Poisson problem does not seem to have appeared in the literature (however, the numerical result µ3D ≈ 0.567 has appeared for example in [13, p. 74]). It is the goal of this thesis to fill a hole in the literature by providing closed form expressions for the smoothing factors of pointwise and block GS-LEX, valid in arbitrary dimensions, for a broad class of difference operators. I will show (Chapter 5) how the smoothing factors of a variety of d-dimensional problems, including Poisson’s equation, steady state diffusion, and a special case of convectiondiffusion, can all be obtained as special cases of a single expression. In particular, the family of problems considered for GS-RB in [17] is such a special case. My analysis relies on results from the theory of analytic functions of a complex variable, in particular the maximum modulus principle and properties of M¨obius Transformations [10]. 1.1 Thesis Contributions The purpose of this thesis is to provide smoothing factors for a family of discrete elliptic operators that includes a number of common model problems. To that end, suppose L u = f is a given elliptic differential equation (assumed linear), and let ∑ LI,hJ uhJ = fIh (1.1) J denote the linear system arising from its discretization on a uniform rectangular mesh, extending to infinity in all directions, so that I and J vary over Zd . We will assume the computational molecule forms a 2d + 1 point star, or equivalently that LI,hJ is of the form 3 Figure 1.1: The computational molecule associated with the operator L h (d = 2). LI,hJ = a −b+ k −b− k 0 if I − J = 0 if I − J = ek (1.2) if I − J = −ek otherwise (here ek denotes the kth basis vector of Rd ). The connection between LI,hJ and the computational molecule is illustrated in Figure 1.1. I assume throughout that LI,hJ satisfies d a≥ ∑ (b+k + b−k ) and − b+ k , bk > 0 for all k ∈ {1, . . . , d}, (1.3) k=1 so that [LI,hJ ] is a diagonally dominant M-matrix [15, p. 85]. In this thesis we will derive analytical smoothing factors for L h of the form (1.2) obeying the constraints (1.3), as well as one additional technical assumption (3.2) in Chapter 3. My analysis relies on results from the theory of analytic functions of a complex variable, in particular the maximum modulus principle and properties of M¨obius Transformations [10]. 1.2 Outline • Chapter 2 provides an overview of multigrid. Special attention is given to the theoretical justification for idealization of the coarse grid correction in Local Fourier Analysis, with a novel interpretation in terms of the NyquistShannon theorem [12] from sampling theory. The more conventional justification is also given, although in somewhat unorthodox terms using equivalence classes. • Chapter 3 is the main analytical portion of the thesis, in which the main results are derived. In particular, Theorem 3.1.5 provides an expression for the smoothing factor of point GS-LEX applied to operators of the form (1.2) 4 obeying the constraints (1.3) as well as the additional constraint (3.2). These results are then extended in Section 3.2 to incorporate the effect of a relaxation parameter. In Theorem 3.3.1, the results are extended in a different direction to block GS-LEX relaxations in which a subset I b ⊂ {1, . . . , d} of indices are relaxed simultaneously. However, in this case the additional complexity is such that our results only apply to the symmetric case where − b+ k = bk for all k, and a relaxation parameter is not included. • Chapter 4, explores the dependence of the smoothing factor on the ratio of the mesh spacing on the fine and coarse grids. An estimate of the total work is included, and the question of the optimal ratio minimizing this work is considered. • In Chapter 5, a series of example problems are considered, and the results of Chapter 3 is used to derive smoothing factors in each case. These smoothing factors are then compared against the actual performance of multigrid in numerical experiments. Differences between theoretical and experimental results are explored where they exist. • Finally, Chapter 6 draws some conclusions. 5 Chapter 2 Overview of Multigrid If one shrinks a sufficiently smooth image down to a reduced size and then blows it back up, it is essentially unaffected. However, the same is not true of a ‘rough’ image, in which the color changes rapidly from pixel to pixel - see Figure 2.1. Figure 2.1: Top: a ‘smooth’ image is shrunk by a factor of 4 in each direction, and then blown back up. Bottom: the same is done for a ‘rough’ image. Reduction and subsequent enlargement were done using Photoshop’s bilinear image scaling option. 6 A similar situation arises in the iterative solution of difference equations arising from the discretization of elliptic partial differential equations (PDEs). If the error between the current approximation and the exact solution is sufficiently smooth, one may project the problem onto a coarser mesh (multigrid researchers call this restriction), solve for the error there, and interpolate it back to the original mesh (prolongation), obtaining a good approximation to the error on the fine mesh. This process is called coarse grid correction. Let L u = f denote the PDE of interest (assumed linear) and let L h uh = f h (2.1) denote the linear system arising from its discretization; in this thesis it is always assumed that the discretization takes place on a uniform rectangular grid. If uˆh denotes the current approximate solution, then the error eh = uh − uˆh can be computed by solving the residual equation L h eh = r h , (2.2) where rh ≡ f h − L h uˆh . The coarse grid correction calculates an approximate error e˜h by computing rh on the fine grid, restricting it to a suitable coarse grid, solving the residual equation (2.2) on the coarse grid using the restricted residual as the right hand side, and prolonging the result back up Figure 2.2: Radii Petrovich Fedorenko (1930-2009) is updated via uˆh → uˆh + e˜h . Note that if e˜h was credited with inventing the multigrid algorithm. the true error, this update would yield the exact to the fine grid. The current guess uˆh is then answer uh . 7 Figure 2.3: A single iteration of two-level multigrid applied to the twodimensional Poisson equation on the square (0, 1)2 , discretized on a 23 × 23 uniform grid using centered differences. The initial error, containing both high and low frequencies, is shown on the left. In the middle we see the error after the smoother S (in this case three iterations of lexicographic Gauss-Seidel) is applied - the high frequency components of the error have been significantly dampened. On the right we see the new error after the coarse grid correction K has been applied all frequencies are now highly reduced. This process may be represented symbolically as e˜h uˆh → uˆh + IHh (L H )−1 IhH ( f h − L h uˆh ), (2.3) rh where IhH and IHh denote restriction and prolongation operators respectively between the fine grid (mesh spacing h) and the coarse grid (mesh spacing H), and L H denotes an appropriate representation of the operator L on the coarse grid. The quantity h/H is called the mesh ratio. Unless specified otherwise, in this thesis it is assumed that H = 2h - see Chapter 4 for a discussion of alternatives. We cannot expect (2.3) alone to be effective, since it requires the error to be smooth. However, if (2.3) is combined with a procedure for efficiently smoothing the error, a highly effective iterative method is born: one smooths the error, applies the coarse grid correction, smooths the new error, applies the coarse grid correction again, and so on (an illustration is provided in Figure 2.3). This idea is the basis of multigrid [11, 13, 16]. Credit for the inventing the algorithm is given to the late Soviet mathematician Radii Petrovich Fedorenko 8 (Figure 2.2), who wrote the first paper on the subject, “A Relaxation Method for Solving Elliptic Difference Equations”, in 1962 [8]. Achi Brandt popularized the method through his work in the 1970s, most notably his seminal 1977 paper “Multi-Level Adaptive Solutions to Boundary Value Problems” [3]; see also [1], [13, p. 23] for a historical account. We represent one iteration of this joint process symbolically as uˆh → Sν2 Sν1 uˆh + IHh (L H )−1 IhH ( f h − L h Sν1 uˆh ) , (2.4) where S is the smoothing operator of choice and ν1 , ν2 are the number of smoothing steps before and after the coarse grid correction. Typical values are ν1 = ν2 = 1 or ν1 = 1, ν2 = 0. The smoother S is typically a simple relaxation scheme such as damped Jacobi or Gauss-Seidel, which, while not necessarily effective at reducing the magnitude of the error, has a strong smoothing effect on it. Equation (2.4) represents what is called two-level multigrid, since it involves only two grids. It is not particularly efficient - discounting smoothing sweeps, two-level multigrid turns the solution of a problem with mesh spacing h into the solution of a small number of problems with mesh spacing H, leading to a speed up of at most a constant factor. Depending on what method is used to solve the coarse problem, there may be no benefit at all. For example, if the coarse grid problem is solved naively with Gaussian elimination, this is roughly 23d (d denotes the dimension of the problem) times cheaper on the coarse grid than it is on the fine grid. Clearly, convergence must be achieved in fewer than 23d iterations for there to be any gains in this case. Therefore, while two-level multigrid is of theoretical interest because of its relative simplicity, it is of little practical use. Full multigrid improves upon its two-level counterpart by (1) acknowledging that it is only necessary to the solve the coarse grid problem approximately and (2) doing so recursively with multigrid. The operator (L H )−1 (whose role is to solve for the error on the coarse grid) is replaced with a small number γ (called the cycle index) of iterations of multigrid. Full multigrid thus involves a whole hierarchy of grids - the coarsest of which is sufficiently small that the error may be cheaply computed by direct methods. See Figure 2.4 for a sketch of the sequence of grids visited in a single iteration of full multigrid (with four grids) for γ = 1 (called a 9 Figure 2.4: The sequence of grids visited by one iteration of four level multigrid, for two different choices of γ. (a) γ = 1: a V-cycle (b) γ = 2: a W-cycle. V -cycle) and γ = 2 (called a W -cycle). Unlike two-level multigrid, it can be shown that under certain conditions the work required for a single iteration of full multigrid (henceforth referred to as multigrid) is proportional to the number of grid points N. In particular, let us assume that the work required for one iteration of the smoother and intergrid operators is linear in N (as is the case for standard intergrid operators and pointwise smoothers) and that the cycle index obeys γ < 2d , where d denotes the dimension of the problem. If we let L denote the number of levels and assume H = 2h, then we find that the total work W obeys W = CN + γCN/2d + γ 2CN/22d + . . . + γ L−1CN/2(L−1)d ≤ CN + γCN/2d + γ 2CN/22d + . . . CN = . 1 − 2−d γ The cost of a single iteration is thus linear in N, just like alternate iterative methods such as Gauss-Seidel, SOR, or the Conjugate Gradient method, provided the aforementioned conditions are met. For most iterative methods, however, the total work, computed as the work per iteration times the number of iterations required for convergence, is not linear in N. This is because the number of iterations required for convergence is typically not a constant, but grows with N. Consider, for example, the three dimensional Poisson problem discretized on a uniform rectangular grid with mesh spacing h. 10 In this case N is proportional to h−3 , Gauss-Seidel and optimal SOR reduce the error by 1 − O(h2 ) and 1 − O(h) per iteration [2, p. 181], and it follows, using the Taylor expansion log(1 + ε) = ε + O(ε 2 ), that the number of iterations required for convergence grows like N 2/3 and N 1/3 respectively. The total work is then proportional to N 5/3 and N 4/3 respectively. By contrast, if multigrid with lexicographic Gauss-Seidel (GS-LEX) as a smoother is used (on the aforementioned 3D Poisson problem), it can be shown that the number of iterations required for convergence does not depend on h [9], and therefore the total work is linear in N. Conversely, one can easily show that no algorithm can do better than O(N)1 , so multigrid is optimal in terms of its scalability. However, optimal scaling is not enough to demonstrate the suitability of multigrid as a numerical method for real world problems. For example, a multigrid algo10 rithm that solves the 3D Poisson problem in a constant 1010 iterations is optimal in terms of its scaling, but will be beaten by simpler methods such as Gauss-Seidel unless the size of the problem is astronomical. With this in mind, we turn our attention to the quantitative convergence analysis of multigrid. 2.1 Quantitative Convergence Analysis The error in a multigrid algorithm with + 1 grids (labeled 0 through from coarsest to finest) obeys the recurrence relation e(k) = M e(k−1) , where e(k) denotes the error after k iterations and M is the iteration matrix of -level multigrid, defined recursively by M = Sν2 (I − I −1 (I −1 − (M −1 )γ )(L −1 −1 L −1 ) I )Sν1 ≡ Sν2 K Sν1 , where I is the identity operator on the th grid, I −1 and I −1 (2.5) are interpolation and projection operators respectively between grids and − 1, L is the discrete representation of L on grid , and M0 = 0 [13, p. 48]. The operator K is called the coarse grid operator. To gain quantitative predictions of convergence rates we are interested in the 1 Any algorithm must specify the value of the solution at each of the N grid points. Even if this can be done in O(1) operations per grid point, this is still O(N) overall work. 11 spectral radius of the iteration matrix, defined by ρ(M ) = max |λi |, i where {λi } denotes the spectrum of M . The spectral radius provides information on the asymptotic convergence rate, via the identity ρ(M ) = lim k→∞ e(k) e(k−1) . Although we are interested in ρ(M ), it is worth noting the limitations on what it can tell us. The spectral radius only provides us with asymptotic information, relevant after a sufficiently large number of iterations have occurred. However, it can in general say nothing regarding what happens before the asymptotic behavior takes over, or how long this will take. Let us digress briefly from multigrid in particular to iterative methods in general, and suppose Ax = b is solved iteratively by a process in which the error obeys e(k) = Me(k−1) , for a suitable matrix M. Even if ρ(M) < 1, the norm of the error may actually increase in the initial phases of the process (see Figure 2.5(b)). In extreme cases, overflow may occur before asymptotic behavior takes over. A quantity which has more to say regarding early behavior (and is unfortunately often much harder to compute) is the norm of the iteration matrix, defined by M ≡ max x=0 Mx , x where · is a suitable vector norm. This quantity provides a bound on the norm 12 Figure 2.5: ||e(k) ||2 vs. k for two-level multigrid with lexicographic SOR 1 smoothing applied to anisotropic diffusion problem 100 uxx + uyy + uzz = 0, discretized using centered differences. (a) ω = 1. (GS-LEX) (b) ω = 1.63 (SOR-LEX). Note that for ω = 1, the norm of the error goes down monotonically, whereas for ω = 1.63 there is a initial increase before the asymptotic regime is entered. Convergence is faster overall for ω = 1.63, despite the erratic behavior of the error initially. of the error after k iterations via the identity: e(k) ≤ M k e(0) ≤ M k e(0) . In particular, if M < 1, then monotonic reduction in the norm of the error is guaranteed. On the other hand, if M > 1, the quantity M k can sometimes provide useful information on the duration of the transient phase before asymptotic convergence takes over - see for example [7], where this is done in the context of Gauss-Seidel with various ordering strategies applied to the one-dimensional convection-diffusion equation. It is worth noting the inequality ρ(M) ≤ M , which becomes an equality when the matrix M is normal, that is M ∗ M = MM ∗ (here M ∗ denotes the conjugate transpose of M). In this special case ρ(M) gives a global upper bound on convergence rates. An example of this occurs when the matrix A is symmetric (for example, the discrete Poisson problem), and damped Jacobi is used as an iterative solver. In this case the iteration matrix is also symmetric, and therefore normal. 13 However, iteration matrices of GS-LEX are generally not normal, and neither are the iteration matrices of multigrid with GS-LEX smoothing. As a result, the convergence of multigrid is in general not guaranteed to be monotonic, and Figure 2.5 shows that it isn’t in general. Having seen that care must be taken interpreting the information provided by the spectral radius, we now move on to the method of its computation. The definition of iteration matrix given in (2.5) is typically too complex to derive analytic convergence rates from. Bounds may be obtained (see for example [13, p. 77]), but they typically lead to estimates of the work which are several orders of magnitude too large [3]. A practical alternative is Local Fourier Analysis, first described by Brandt [3]. Although non-rigorous, this approach can yield estimates of ρ(M ) that are much closer to the actual convergence rate of multigrid than the bounds provided by the rigorous approach [3]. As stated in the Introduction, the main idea behind Local Fourier Analysis is to idealize the effect of the coarse grid correction. Representing the error as a sum of Fourier modes, one partitions the modes into ‘low frequency’s (also called ‘smooth modes’) and ‘high frequency’s (‘rough modes’), and assumes the coarse grid correction leaves the high frequencies alone while wiping the low frequencies out completely. The smoothing factor µ is then defined in terms of this idealized coarse grid operator Kˆ by the formula µ≡ ν ρ(Sν2 Kˆ h Sν1 ), where ν = ν1 + ν2 . µ represents the average asymptotic error reduction factor per smoothing sweep. If bounds independent of h are desired, an alternative definition is µ ≡ sup ν ρ(Sν2 Kˆ h Sν1 ). h Smoothing factors are often analytically computable in situations where ρ(M ) is not. In particular, since Kˆ is by definition diagonal when written with respect to the Fourier basis, a highly tractable situation arises when the same is true of 14 the smoother S. In this case the smoothing factor reduces to the largest magnitude eigenvalue of the high frequency eigenfunctions of S. It is intuitively clear that a Fourier mode whose wavelength is as long as the grid should be classified as low frequency and conversely, a mode with a wavelength on the order of the mesh spacing h should be called high frequency. It is also clear that the former can be well represented on the coarse grid while the latter cannot (see Figure 2.1) - therefore, it is reasonable to idealize the coarse grid correction as totally effective and totally ineffective respectively in these cases. What isn’t so clear is where the line should be drawn between high and low frequencies, and whether it is realistic to draw such a line at all. 2.2 A Closer Look at the Interplay Between Errors on the Fine and Coarse Grids I conclude this chapter with a theoretical justification for the partitioning of modes into high and low frequencies, and their subsequent treatment by the idealized coarse grid correction. First, an (essentially) standard argument is provided based on the observation that distinct Fourier modes alias when restricted to the coarse mesh. It is followed by a novel argument based on the Nyquist-Shannon Theorem from sampling theory. For simplicity, let us assume our error lives on a one dimensional grid Gh = {h, 2h, . . . , 1 − h} ⊂ (0, 1). On the continuous interval (0, 1), any f ∈ L2 ((0, 1)) may be represented as an infinite sum of Fourier modes belonging to the set Ψ = {ϕn ≡ sin(nπx) : n ∈ N}. Although distinct on (0, 1), modes ϕn and ϕm look the same2 when restricted to Gh provided n ± m is a integer multiple of 2h−1 . As a result, one may show that the infinity of possible modes on (0, 1) collapses down to just N on Gh 3 , where N ≡ |Gh | is assumed odd. These modes span the function space { f : Gh → R}, viewed as a vector space of dimension N. 2 i.e differ at most in sign. 3 Technically there are N +1 possible modes, but one is the trivial mode which is zero everywhere, so I discount it. 15 We formalize this observation by introducing an equivalence relation ∼h on Ψ, where Definition 2.2.1 We say two modes ϕn , ϕm ∈ Ψ are h-equivalent if they cannot be distinguished on Gh . In this case we write ϕn ∼h ϕm . It is easy to verify that ∼h is an equivalence relation, and that ϕn ∼h ϕm ⇔ ϕn Gh = ±ϕm Gh ⇔ n ≡ ±m mod 2h−1 . In light of definition 2.2.1 we may identify the Fourier modes on Gh with equivalence classes of modes on (0, 1). That is, we identify the fourier mode ϕn Gh on Gh with the equivalence class {ϕn :∼h } of modes on (0, 1) modulo the equivalence relation ∼h . The set of distinct Fourier modes on Gh , denoted by Ψh , may then be represented as Ψh = {{ϕn :∼h } : n = 1, . . . , N}. Let us now introduce a second grid GH = {H, 2H, . . . , 1 − H}, where for simplicity H = 2h (note that this ensures GH ⊂ Gh ). It is natural to ask how the equivalence classes ΨH on GH are related to those on Gh - in particular, one might wonder if there are modes which can be distinguished on Gh but not on GH . To answer this, one may check that ϕn ∼H ϕm and ϕn h−1 h ϕm if and only if n + m is a multiple of = N + 1, from which it follows that {ϕn :∼H } = {ϕn :∼h } ∪ {ϕN+1−n :∼h }. In other words, when the mesh spacing is doubled the equivalence classes of Fourier modes merge pairwise; conversely, when the mesh spacing is halved each equivalence class splits in two. Let us now imagine a restriction operator IhH : { f : Gh → R} → { f : GH → R} defined by IhH ( f ) = f 16 GH . Figure 2.6: The distinct Fourier modes ϕ1 (k) = sin(kπ/20) and ϕ19 (k) = sin(19kπ/20) on the fine grid Gh = {k/20}20 k=0 restrict to the same mode ϕ1 (K) = sin(Kπ/10) on the coarse grid GH = {K/10}10 K=0 . Then IhH is an operation that maps Fourier modes on Gh to Fourier modes on GH in a 2 − 1 fashion since IhH ({ϕn :∼h }) = IhH ({ϕN+1−n :∼h }) = {ϕn :∼H }, see Figure 2.6. This fact creates a problem in the phase of the coarse grid correction in which the error on the coarse grid is sent to the fine grid. For if an error of {ϕn :∼H } has just been solved for on GH , we can only conclude that the error on Gh is some linear combination of {ϕn :∼h } and {ϕN+1−n :∼H } - in other words, we do not have enough information to reconstruct the error on Gh . However, if we somehow know that the error on GH contains no modes {ϕn :∼H } with n ≥ N/2 (for example if they had all been removed by a process of smoothing), we can say with confidence that an error of {ϕn :∼H } on GH implies an error of {ϕn :∼h } on Gh . Under these conditions we can construct a prolongation 17 operator IHh : { f : GH → R} → { f : Gh → R} by defining IHh ({ϕn :∼H }) = {ϕn :∼h }, (2.6) and transfer errors between GH and Gh without loss of information. Let us define a Fourier mode to be ‘high frequency’ if n ≥ N/2, and low frequency otherwise. It follows from (2.6) that the range of the prolongation operator IHh is contained in the span of the low frequency Fourier modes, and therefore all high frequency error components are unaffect by the coarse grid correction. On the other hand, restriction combines the low frequency mode {ϕn :∼h } and its high frequency ‘partner’ {ϕN+1−n :∼h } into a single mode on the coarse grid. The coarse grid correction will thus overcorrect {ϕn :∼h } by an amount equal to the amplitude of {ϕN+1−n :∼h }. However, if we assume the amplitude of the high frequency mode {ϕN+1−n :∼h } is sufficiently small, then this overcorrection is negligible and we may assume that the coarse grid correction simply eliminates all smooth modes. For another way of thinking about things, we turn to sampling theory, where the coarse grid operator can be thought of as a process that attempts to construct the error on Gh from a sampling on the subgrid GH . By the Nyquist-Shannon theorem, this is feasible if and only if the sampling frequency is at least twice the highest frequency present in the error (expanded in a Fourier series) [12]. Let us write the error on the grid Gh as a sum of complex Fourier modes of the form eıθ ·x/h ,4 where x = Ih ∈ Gh , and θ = (θ1 , θ2 , . . . , θd ) ∈ [−π, π]d . Then the sampling frequency on GH is 1 H in each direction, the frequency of mode eıθ x/h in direction ek is given by fk = θk /2πh, and we conclude that the coarse grid correction will be successful only if all modes in the error obey θk ≤ h π for all k = 1, . . . , n. H (2.7) Let us assume that the components of the error violating (2.7) (called ‘rough’ modes) have been dampened in the smoothing stage to the point that they are small in comparison with the components that obey (2.7) (called ‘smooth’ modes). Then we idealize the coarse grid operator by assuming it succeeds at reconstructing the 4 This is more convenient than the real Fourier modes used previously, and is preferred in LFA. 18 error in the fullest sense possible given (2.7); the ‘smooth’ components of the error are reconstructed perfectly while the ‘rough’ components are ignored. As a result, after the reconstructed error e˜h is added to the current guess, the ‘smooth’ part of the error is gone (because it was perfectly constructed) and the ‘rough’ part is as it was at the end of the smoothing step. 19 Chapter 3 Smoothing Analysis Let L h denote the elliptic difference operator defined in (1.2) and illustrated in Figure 1.1. Recall that L h is a d-dimensional operator with a stencil that forms a 2d + 1-point star. The central element of L h is denoted by a, while the element one unit in the ±ek direction from a is denoted by −b± k. I will consider an ordering in which grid points are ordered according to the following rule: along dimension k, unknowns are ordered from −ek to +ek if b+ k ≤ b− k , and from +ek to −ek otherwise. For convenience, let us define: − ck ≡ max(b+ k , bk )/a − dk ≡ min(b+ k , bk )/a. It is straightforward to show (following the technique of [13, p. 103] for example) that the smoothing factor of pointwise lexicographic Gauss-Seidel applied to L h is given by µ pt = max θ ∈Θd ∑dk=1 dk eıθk , 1 − ∑dk=1 ck e−ıθk (3.1) where Θd = [−π, π]d \(−π/2, π/2)d is the set of rough modes in d dimensions and √ ı ≡ −1. The argument of the max in (3.1) represents the eigenvalue of the point GS-LEX iteration matrix corresponding to the eigenfunction ϕθ : Gh → C defined by ϕθ (Ih) = eıI·θ . Although not an eigenfunction of GS-LEX on any finite grid, ϕθ becomes an eigenfunction when the grid Gh is assumed infinite, as in LFA. 20 We restrict our analysis to the case dk = αck , (3.2) where α ∈ (0, 1] is a constant independent of k. This restriction includes all cases where [LI,hJ ] is symmetric, and other cases, for example the case of the convectiondiffusion equation with all mesh Reynolds numbers equal. It is well-known that pointwise smoothing is ineffective for highly anisotropic problems [13, p. 131]; one way to resolve this is the use of block smoothers [13, p. 134]. Such smoothers are defined by a partitioning of {1, . . . , d} into disjoint sets Ib ⊂ {1, . . . , d} and I p = {1, . . . , d}\Ib , where coordinates belonging to Ib are relaxed simultaneously. To avoid degenerate cases, we assume both I p and Ib are non-empty. In this case the smoothing factor is given by µ block = max θ ∈Θd ∑k∈I p dk eıθk . 1 − ∑k∈I p ck e−ıθk − ∑k∈Ib (ck e−ıθk + dk eıθk ) (3.3) It is worth noting that the computational cost per iteration of a block smoother is higher than that of a point smoother. However, for highly anisotropic cases this overhead is typically small compared to the gains in convergence rates. Our analysis shows that µ pt depends only on α and any two of the quantities c ≡ c1 + c2 + . . . + cd cm ≡ min(c1 , . . . , cd ) cr ≡ c − cm , (3.4) regardless of the number of dimensions, d. Similarly, in the case that L h is symmetric, our analysis shows that the block smoothing factor µ block depends on any two of the quantities Cp = ∑ k∈I p ck , Cmp = min ck , k∈I p Crp = C p −Cmp as well as any two of Cb , Cmb and Crb (which are defined analoguously). 21 (3.5) Expressions for µ pt and µ block are given in Theorems 3.1.5 and 3.3.1 respectively. 3.1 Pointwise Relaxation Starting from (3.1), we complex-conjugate the denominator and apply (3.2) to obtain µ pt = max θ ∈Θd α ∑dk=1 ck eıθk αg(θ ) ∑dk=1 dk eıθk = max ≡ max , d d −ıθ ıθ 1 − ∑k=1 ck e k θ ∈Θd 1 − ∑k=1 ck e k θ ∈Θd 1 − g(θ ) where g(θ ) ≡ ∑dk=1 ck eıθk . The key to our analysis lies in the observation that this may be rewritten as µ pt = max | f (z)|, z∈g(Θd ) where f (z) ≡ αz 1−z is analytic in the punctured plane C\{1}. If the set g(Θd ) is sufficiently well-behaved (a connected open set or the closure of one) and does not contain the point z = 1, the maximum principle [10, p. 88] applies and we can write µ pt = max | f (z)|, z∈∂ g(Θd ) where ∂ g(Θd ) denotes the boundary of g(Θd ). However, computing g(Θd ) explicitly is challenging because for certain values of {ck } it contains a hole in the vicinity of the origin - see Example 3.1.3. Keeping track of this hole is difficult - we avoid the issue by proving we can replace g(Θd ) with a simply-connected set D equivalent to it in the following sense: Definition 3.1.1 Let A, B ⊂ C. We say A ∼ B if sup | f | = sup | f | B A 22 Figure 3.1: An illustration of Example 3.1.3: (a) If a double pendulum with arms of length 16 and 13 is allowed to swing freely, it remains inside a disk of radius 0.5. (b) If one arm is constrained to lie in the left half plane, the pendulum swings in the heart-shaped region shown. In both cases the pendulum is unable to reach a set of points surrounding the origin. for all functions f complex differentiable on A ∪ B. The following observation establishes an important sufficient condition for two sets to be equivalent. Observation 3.1.2 Let A, B ⊂ C, and suppose B is a bounded, connected open set (domain). If ∂ B ⊆ A ⊆ B, then A ∼ B ∼ B. To see this, let f be any function complex differentiable on A ∪ B = B – since B is a domain, f is analytic on B. We clearly have sup∂ B | f | ≤ supA | f | ≤ supB | f |, but the maximum modulus principle implies the equality of the leftmost and rightmost terms. In light of Observation 3.1.2, our goal is to find a simply-connected domain D such that (i) D = E, where E is a domain (ii) ∂ D ⊆ g(Θd ) ⊆ D. This is done in Lemma 3.1.4, but first we build some intuition with a simple example. Example 3.1.3 Suppose d = 2 and L h arises from a centered differences discretization of the Laplacian-like operator L = 0.5∂xx + ∂yy , giving a = 3, b+ 1 = 23 − + b− 1 = 0.5 and b2 = b2 = 1. Then c1 = 1/6, c2 = 1/3, and 1 1 g(θ1 , θ2 ) = eıθ1 + eıθ2 . 6 3 Geometrically, g(θ1 , θ2 ) may be viewed as a double pendulum with arms of length 1/6 and 1/3, making angles θk (k=1,2) with the x-axis. As (θ1 , θ2 ) varies over [−π, π]2 , this pendulum (and therefore the range of g) remains confined to the disk of radius 0.5 shown in Figure 3.1(a). Furthermore, the boundary of the disk is swept out as the double pendulum completes a full revolution with a constant angle of 180◦ between the arms. By Observation 3.1.2, g([−π, π]2 ) is equivalent to this disk in the sense of definition 3.1.1. However, if (θ1 , θ2 ) is constrained to lie in Θ2 ≡ [−π, π]2 \[− π2 , π2 ]2 , then |θk | ≥ π/2 for at least one k ∈ {1, 2}. In other words, at least one arm is constrained to lie in the left half plane. With this constraint, g(θ1 , θ2 ) now lives in the set shown in Figure 3.1(b) - the union of two disks and a half disk, as demonstrated in Figure 3.2. It is simple to show the boundary of this set is once again in the range of g, so this heart-shaped set provides us with the D we are looking for. Note that while g(Θ2 ) ∼ D it is not true that g(Θ2 ) = D, since 0 ∈ D but |g| ≥ 1 3 − 61 > 0. Lemma 3.1.4 Let B(r, z) denote the closed ball of radius r ≥ 0 centered at z ∈ C. Let C− denote the set of all complex numbers with negative real part, and let ∼ be the relation in definition 3.1.1. Then g(Θd ) ∼ [B(c, 0) ∩ C− ] ∪ B(cr , ıcm ) ∪ B(cr , −ıcm ) ≡ D, g([−π, π]d ) ∼ B(c, 0), where c, cm , and cr are defined in (3.4). Proof. By Observation 3.1.2 we clearly have g([−π, π]d ) ∼ B(c, 0), since |g(θ1 , . . . , θd )| ≤ c and g(θ , . . . , θ ) = ceıθ . On the other hand, if θ ∈ Θd then there is a j such that |θ j | ≥ π/2. Assuming we have made some fixed choice of j and set θ j to a fixed value α, the remaining variables {θk }k= j range over [−π, π]d−1 . It 24 follows that g is confined to B(c − c j , c j eıα ), and therefore d g(Θd ) ⊆ B(c − c j , c j eıα ). |α|≥ π2 j=1 But B(c − c j , c j eıα ) ⊆ B(cr , cm eıα ) for all j since |z − c j eıα | ≤ c − c j implies |z − cm eıα | ≤ |z − c j eıα | + |c j eıα − cm eıα | ≤ (c − c j ) + (c j − cm ) = cr . Therefore g(Θd ) ⊆ B(cr , cm eıα ). |α|≥ π2 It can be shown that any z ∈ B(cr , cm eıα ) obeying Re(z) ≥ 0 lies in B(cr , ıcm ) if Im(z) ≥ 0 and B(cr , −ıcm ) if Im(z) ≤ 0. Similarly, since |z| ≤ c we have z ∈ B(c, 0) ∩ C− if Re(z) < 0. It follows that g(Θd ) ⊆ D. It is trivial to show ∂ D ⊆ g(Θd ), which completes the proof. Lemma 3.1.4 and Definition 3.1.1 allow us to conclude µ pt = max | f | (3.6) ∂D provided f is complex differentiable in D, that is, provided we can show 1 ∈ / D. To that end, we note that the diagonal dominance of L h (1.3) implies c = ≤ − ∑dk=1 max(b+ k , bk ) a − d ∑k=1 max(b+ k , bk ) − ∑dk=1 (b+ k + bk ) Figure 3.2: The heart-shaped set D is the union of two disks and a half disk. < 1. The desired result follows from the obser25 vation D ∩ R = [−c, c2r − c2m ] ⊆ [−c, c). Theorem 3.1.5 The smoothing factor of pointwise GS-LEX applied to discrete operators of the form (1.2) obeying the constraints (1.3) and (3.2) is given by µ pt (cm , cr , α) = αcr + α c2m + (c2m − c2r )2 , 1 + c2m − c2r (3.7) where cm and cr are defined in (3.4) and α is defined in (3.2). Proof. The set ∂ D is the union of the three semicircular arcs π 3π , 2 2 S1 = ceıθ : θ ∈ S3 = z ∈ C : z¯ ∈ S2 . , S2 = ıcm + cr eıθ : θ ∈ − sin−1 cm cr , π 2 , By (3.6), it suffices to compute the maximum of the maxima of | f | over S1 , S2 , and S3 . However, since | f (¯z)| = | f (z)|, the maxima over S2 and S3 are the same and hence we omit S3 from consideration. Furthermore, it is trivial to show that the maximum of | f | over S1 is achieved at ıc, the point of intersection of S1 and S2 . Therefore, the maximum of | f | over S2 is at least as large as the maximum over S1 - hence we may omit S1 as well. It follows that the maximum of | f | over D is attained on the semicircle S2 . However, it is also true that the maximum is attained on the full circle |z−ıcm | = cr , as this circle contains S2 and is contained in D. It is more convenient to work with the full circle, so we conclude µ pt = Now, f (z) = αz 1−z max | f (z)|. |z−ıcm |=cr (3.8) is a M¨obius transform, and hence the image of a circle is either a circle or a line [10, p. 65]. Since the circle |z − ıcm | = cr does not contain any poles of f , it follows that f ({|z − ıcm | = cr }) is a circle. If we let zc ∈ C and r ≥ 0 denote the center and radius respectively of this 26 circle, then from (3.8) we have µ pt = max z∈ f ({|z−ıcm |=cr }) |z| = |zc | + r. (3.9) To find the parameters zc and r, we first write f as a sequence of elementary M¨obius transformations: f = f4 ◦ f3 ◦ f2 ◦ f1 , where f1 (z) = z − 1, f2 (z) = z−1 , f3 (z) = −αz, and f4 (z) = z − α. Then, starting with the circle |z − ıcm | = cr , we track the changes in its center and radius as we compute its image under f1 , then compute the image of the result under f2 , and so on. Since f1 , f3 , and f4 are all either translations or dilations, the steps involving them are straightforward. For f2 , we make the observation that the image of a circle with radius R and center z0 under inversion is the circle with radius R and center z0 given by z0 = |z0 z¯0 2 | − R2 R and R = ; |z0 |2 − R2 this fact can easily be derived from the argument found in [10, p. 66]. In the end we find r= and zc = α αcr 1 + c2m − c2r cm c2r − c2m + ıα , 2 2 1 + cm − cr 1 + c2m − c2r which completes the proof. Notice that the smoothing factor µ pt is proportional to α. Since the latter is small when the operator L h is strongly asymmetric, this suggests that GS-LEX may be particularly effective in this regime. In particular, as α → 0 the part of L h above its diagonal is converging to the zero matrix, and for α = 0 we have that L h is lower triangular and GS-LEX becomes a (direct) solver. See also Chapter 5 for further illustrations of this behavior. 27 3.2 SOR Relaxation It is well-known that the smoothing factors of GS-RB, Jacobi-RB1 , and Jacobi can be greatly improved by incorporating a relaxation parameter [13, 18, 20]. For example, for the standard seven-point centered difference discretization of the threedimensional Poisson problem, overrelaxation with ω = 1.15 improves the smoothing factor of GS-RB (equivalent to Jacobi-RB in this specific case) from µ ≈ 0.44 to µ ≈ 0.23 [18]. For the same problem, the smoothing factor of Jacobi is improved from µ = 1 (no convergence) to µ = 5/7 by underrelaxation with ω = 6/7 [13, p. 73]. We can analyze the effect of a relaxation parameter ω in GS-LEX smoothing by repeating the steps of Theorem 3.1.5 with the M¨obius transformation fω (z) = αz + ω1 − 1 1 ω −z 1 c2r −c2m used in place of f . Note that we must assume ω < √ in order to ensure that fω has no poles in D. Making use of the factorization fω = f4 ◦ f3 ◦ f2 ◦ f1 , where f1 (z) = z − ω1 , f2 (z) = z−1 , f3 (z) = 1 − 1+α z, f4 (z) = z − α, we obtain ω pt (ω) µSOR = cr 1 − 1+α ω + Figure 3.3(a) shows a plot of c2m 1 − 1+α ω 2 1 ω 1 + c2m − c2r ω2 pt µSOR + 1 − ω1 + α(c2m − c2r ) 2 . (3.10) as a function of ω ∈ [0, 2] for the 3D Poisson problem mentioned above (cm = 1/6, cr = 1/3, α = 1). The optimal smoothing factor µ ≈ 0.551 is attained at ω ≈ 1.1, a modest improvement over the smoothing factor µ ≈ 0.567 obtained when ω = 1. This finding supports the conclusion [13, p. 105] that for GS-LEX, the inclusion of a relaxation parameter is not necessarily worth the extra work per iteration (two operations per point per relaxation sweep). In Figure 3.3(b), the anisotropic problem −εuxx − uyy − uzz = f is discretized using standard seven-point centered differences, and the theoretical smoothing fac1 Jacobi-RB consists of a Jacobi sweep over the red points, followed by a Jacobi sweep over the black points using the updated values at red points; see for example [13, p 173]. 28 Figure 3.3: Left: Smoothing factor of SOR-LEX with relaxation parameter ω, as a function of ω, for the standard seven-point centered differences discretization of the 3D Poisson problem. Right: Smoothing factors of GS-LEX and optimal SOR-LEX applied to the anisotropic diffusion problem −εuxx − uyy − uzz = f , as a function of ε, assuming seven-point centered differences discretization. tors of both GS-LEX and optimal SOR-LEX (obtained by numerical minimization of (3.10)) are plotted against ε. 3.3 Block Relaxation In this section we analyze the smoothing properties of block GS-LEX relaxation. To keep the algebra manageable, we restrict the scope of our analysis to the case that [LI,hJ ] is symmetric (α = 1). Our main result is Theorem 3.3.1, which establishs a connection between block and pointwise smoothing factors. Since α = 1, we have ck = dk for all k - substituting this into equation (3.3), complex-conjugating the denominator and then simplifying, we obtain µ block = max θ ∈Θd = ∑k∈I p ck eıθk 1 − 2 ∑k∈Ib ck cos θk − ∑k∈I p ck eıθk max (z,x)∈G(Θd ) z , 1−x−z 29 where G : Θd → C × R is given by G(θ ) = (G1 (θ ), G2 (θ )) with G1 (θ ) = ∑ ck eıθk and G2 (θ ) = 2 k∈I p ∑ ck cos θk . k∈Ib Defining z = z/(1 − x), this becomes µ block = max z ∈Λ z = max | f (z )|, 1−z z ∈Λ (3.11) where f is defined in section 3.1 and Λ ≡ {z/(1 − x) : (z, x) ∈ G(Θd )}. We will show that Λ is ∼ (in the sense of Definition 3.1.1) to the union of a scaled copy of the set D from Lemma 3.1.4 and a particular ball centered at the origin. This fact will allow us to express µ block in terms of µ pt . If θ ∈ Θd , then there is a j such that |θ j | ≥ π/2. If j ∈ I p , it follows that {θk }k∈I p ∈ Θ|I p | and {θk }k∈Ib ∈ [−π, π]|Ib | , while if j ∈ Ib , then {θk }k∈I p ∈ [−π, π]|I p | and {θk }k∈Ib ∈ Θ|Ib | . The function G1 has the same form as g of Section 3.1 with the variables {θk }dk=1 replaced by {θk }k∈I p . Thus, Lemma 3.1.4 applies with C p , Cmp and Crp playing the roles of c, cm , cr : D(Cmp ,Crp ) G1 (Θd ) ∼ B(C p , 0) if k ∈ I p , otherwise. At the same time we have [−2Cb , 2Cb ] G2 (Θd ) = [−2Cb , 2Cb ] r 30 if k ∈ I p , otherwise. Assuming j ∈ I p , since G1 and G2 depend on disjoint variables we have Λ 1 G1 (Θd ) 1−x = j∈I p ∼ j∈I p x∈G2 (Θd ) x∈[−2Cb ,2Cb ] 1 D(Cmp ,Crp ). 1−x j∈I p Note that if 2Cb < 1, then for all x ∈ [−2Cb , 2Cb ] we have 1 1 D(Cmp ,Crp ). D(Cmp ,Crp ) ⊆ 1−x 1 − 2Cb It follows from the diagonal dominance of L h and the condition α = 1 that c ≤ 0.5. Since Cb < c, we have 2Cb < 1 as desired, and therefore Λ ∼ j∈I p Cmp Crp 1 p p ,C ) = D D(C , . m r 1 − 2Cb 1 − 2Cb 1 − 2Cb (3.12) A similar analysis with j ∈ Ib shows Λ ∼B j∈Ib Cp ,0 . 1 − 2Crb Theorem 3.3.1 The smoothing factor of block GS-LEX applied to symmetric discrete operators of the form (1.2) obeying the constraints (1.3) and (3.2) is given by µ block = max µ pt (c˜m , c˜r , α = 1), where Cp 1 −C p − 2Crb , (3.13) Cmp Crp c˜m = , c˜r = , 1 − 2Cb 1 − 2Cb and C p ,Cmp ,Crp ,Cb ,Crb are defined in (3.5). Proof. Clearly the maximum of | f | over Λ is the maximum of the separate maxima over Λ j∈I p and Λ j∈Ib . By (3.12) and Theorem 3.1.5 we have max | f (z)| = µ pt (c˜m , c˜r , 1) , z∈Λ| j∈I p provided f is analytic in D(c˜m , c˜r ). From the inequalities C p ,Cb < c ≤ 1/2 it is 31 easy to show c˜ ≡ c˜m + c˜r ≤ 1, and the desired result follows in the same way as in Section 3.1. Similarly, since Crb ,C p < 1/2 we have 1∈ /B Cp ,0 ∼ Λ 1 − 2Crb . j∈Ib It follows from the maximum modulus principle that f (z) attains its maximum modulus over Λ j∈Ib at z = C p /(1 − 2Crb ), and therefore max | f (z)| = z∈Λ| j∈Ib Cp , 1 −C p − 2Crb which completes the proof. Theorem 3.3.1 shows that applying a block smoother can be equivalent to applying a point-smoother to a related problem. 32 Chapter 4 Alternative Mesh-Ratios In his comprehensive 1977 paper [3], Brandt had the following to say regarding the choice of mesh ratio h/H. “...it is evident that the mesh-size ratio (H = 2h) is close to optimal...this ratio is more convenient and more economical in the interpolation process than any other efficient ratio. In practice, therefore, the ratio (H = 2h) should always be used, giving also a very desirable standardization.” Indeed, the choice H = 2h (called double coarsening) is very common amongst multigrid practitioners. However, under certain conditions alternatives may be more appropriate. For example, it has been argued that for cell-based discretization schemes triple coarsening (H = 3h) is a better choice because it better preserves cell structure [6]. See also [20] for an example involving the choice H = 4h. Therefore, to keep my results as general as possible, this chapter considers general mesh ratios h/H = β ∈ (0, 1). The smoothing factor of point GS-LEX is derived as a function of β , along with an estimate of the total work of multigrid. One finds that the standard choice β = 0.5 is nearly optimal when the dimension is small. In Theorem 4.0.2 I will show that a judicious choice of β can reduce the overall work by up to 50% in high dimensions, relative to the standard choice H = 2h; I call this optimal coarsening. High dimensional Poisson problems are of interest in applications such as mathematical finance; their efficient numerical solution by multigrid has been studied in [20]. It is possible that optimization of the mesh ratio could lead to reductions 33 Figure 4.1: (a) An example pacman set with c = 0.5 and β = 1/4. (b) The corresponding set Dβ when cm = cr = 0.25. This example corresponds the standard centered differences discretization of the two dimensional Poisson problem, with the coarsening strategy H = 4h. in the overall work for problems of this type. However, since the optimal β (denoted by βopt ) is typically such that h < Hopt ≡ h βopt < 2h, a potentially significant price is that the coarse grid nodes do not line up with those of the fine grid, leading to complications in the intergrid transfer operators (see Figure 4.3). The modified smoothing factor µ pt may be computed as a function of β from (3.1), replacing the rough modes Θd with their β -dependent generalization derived from (2.7), namely Θdβ = [−π, π]d \(−β π, β π)d . It is not hard to generalize the result of Theorem 3.1.5 to include β = 1/2. The set D is replaced with a set Dβ composed of two disks, centered at cm e±ıβ π , and a ‘pacman’ shape with mouth subtending an angle of 2β π (see Figure 4.1)1 . In symbols, we have Dβ = P(c, β ) ∪ B(cr , cm eıβ π ) ∪ B(cr , cm e−ıβ π ), 1 To see this, simply follow the proof of Lemma 3.1.4, replacing occurrences of π/2 with β π. 34 where P(c, β ) is the pacman set P(c, β ) ≡ {reıθ : r ≤ c and θ ∈ [−π, π]\[−β π, β π]}. Mimicking the approach taken in Chapter 3, the smoothing factor is found by taking the image of the circle |z − cm eıβ π | = cr under the M¨obius transform f (z). In the end one finds µ pt (β ) = αcr + α c2m + 2cm (c2r − c2m ) cos(β π) + (c2r − c2m )2 . 1 + c2m − c2r − 2cm cos(β π) (4.1) The qualitative behavior of µ pt as a function of β may not be obvious from (4.1). However, noting that β > β implies Θdβ ⊂ Θdβ , it follows that µ pt is a decreasing function of β . Therefore, the larger β is, the fewer iterations are required for convergence. On the other hand, the larger β is, the more grids are required between the coarsest and finest. As an example, suppose we are in two dimensions and the finest grid has 1002 points. With the mesh ratio H = 2h, five grids are required to obtain a coarse grid with fewer than 52 unknowns; for H = 1.5h, however, eight grids are required. Thus, the work per iteration increases with β . A natural question now arises: which β ∈ (0, 1) minimizes the total work, computed as the product of the work per iteration and the number of iterations required for convergence? Let N denote the number of grid points on the finest grid. Let us assume for simplicity that V -cycles are employed, and consider the limit where the number of grids tends to infinity. Then the work per iteration is wi = C(N + β d N + β 2d N + . . .) = CN , 1−βd where C is a constant of proportionality determined by the smoother and intergrid operators. At the same time, the number of iterations required to bring the error down to a predetermined tolerance is proportional to log(µ pt )−1 . It is natural to express the total work as a fraction of the work for the standard 35 Figure 4.2: The relative work Wβ /W0.5 vs β for the d-dimensional Poisson problem with d = 2, 3, 4. βopt Wopt /W0.5 d=2 0.57 0.97 d=3 0.65 0.86 d=4 0.7 0.77 d=8 0.76 0.63 d = 16 0.82 0.56 d = 32 0.88 0.53 d = 64 0.92 0.51 Table 4.1: Relative work of optimal coarsening and corresponding values of βopt , for various d-dimensional Poisson problems. choice β = 0.5. One obtains Wβ (1 − 0.5d ) log(µ pt (0.5)) = . W0.5 (1 − β d ) log(µ pt (β )) Figure 4.2 shows a plot of Wβ W0.5 (4.2) for a few d-dimensional Poisson problems. Table 4.1 displays the optimal β and corresponding savings for various values of d. From Figure 4.2 and Table 4.1, two things are apparent. Firstly, the choice β < 0.5 (that is H > 2h) seems to lead to an increase in the overall work. Thus, mesh ratios such as H = 3h or H = 4h will tend to increase the overall work (but may be desirable for other reasons). Second, while values of β > 0.5 (h < H = h β < 2h) can lead to decreased work, the gains appear to be small in low dimensions - in particular, for the two-dimensional Poisson problem optimal coarsening yields 36 a theoretical reduction of only 3% in the overall work. Figure 4.3: A fine and coarse grid with mesh ratio h/H = 2/3. The stencil of the projection operator varies from point to point, depending on whether the coarse grid lines up with fine grid. For three-dimensional Poisson, however, the gain is almost 15%, and may be worthwhile. To give an idea of how this method might be implemented in practice, Figure 4.3 provides a sketch of a one-dimensional fine and coarse grid with h/H = 2/3 (sufficiently close to βopt to give good results, but more convenient in terms of alignment of gridpoints), and gives a possible restriction operator. It is clear from Table 4.1 that the advantages of optimal coarsening become more significant as the dimension grows. In fact, we have the following Theorem: Theorem 4.0.2 Consider the standard centered differences discretization of the Poisson problem in d dimensions. As d → ∞, βopt → 1 and Wopt → 0.5. W0.5 Proof. As d → ∞, 1 − 0.5d → 1, while 1 − β d converges to the step function which is 1 for β ∈ [0, 1) and 0 for β = 1. As a result, (4.2) becomes Wβ d→∞ W0.5 lim and one finds pt log(µ (0.5)) log(µ pt (β )) = ∞ if β < 1 if β = 1 Wopt log(µ pt (0.5)) = , d→∞ W0.5 log(µ pt (1)) lim 37 while βopt tends to 1 (see Figure 4.4). One may expand µ pt (β ) as a Taylor series in 1/d, obtaining 1 1 µ pt = 1 − 2(1 − cos(β π)) + O 2 . d d Then, using the identity log(1 + x) = x + O(x2 ) (valid for x 1), one obtains Wopt = 0.5 d→∞ W0.5 lim for the d-dimensional poisson problem. Figure 4.4: The relative work Wβ /W0.5 for the d-dimensional Poisson problem converges to the limit curve shown as d → ∞. 38 Chapter 5 Examples In this section we apply Theorems 3.1.5 and 3.3.1 to a few examples, obtaining formulas for the smoothing factors of standard discretizations of several common PDEs; these are then compared with the measured asymptotic convergence rate of multigrid with GS-LEX smoothing. We also include measured convergence rates with GS-RB smoothing, and show that while the latter has better smoothing properties for the symmetric case [17], lexicographic smoothing is more effective in the strongly asymmetric setting of the convection-dominated convection-diffusion equation. All experiments are done on a rectangular domain with homogeneous Dirichlet boundary conditions. Unless stated otherwise, we use the Galerkin coarse h in our experiments, where I h is the prolongation opgrid operator LH = Ih2h Lh I2h 2h h )T is the restriction operator with erator with linear interpolation, and Ih2h = 2d (I2h full weighting. The asymptotic convergence rate ρ is estimated from the sequence of residuals {r(k) }m k=0 using the identity ρ = lim m m→∞ r(m) r(0) 2 2 and approximating the limit with a sufficiently large value of m (see [13, p. 54] for justification). All experiments were done in M ATLAB 7.9.0 on a 1.80 GHz Intel(R) Core(TM)2 Duo CPU. 39 d=1 0.447 relaxation/dimension point relaxation line relaxation plane relaxation d=2 0.5 0.447 d=3 0.567 0.5 0.447 Table 5.1: Smoothing factors of GS-LEX with point, line, and plane relaxation for the Poisson problem in 1D, 2D, and 3D. Moving from point to line or from line to plane relaxation is equivalent to reducing the dimension by one. 5.1 The d-Dimensional Poisson Problem Suppose the d dimensional Poisson problem −∆u = f is discretized using centered differences on a uniform grid. Theorem 3.1.5 gives a general formula for the smoothing factor of pointwise GS-LEX, namely √ 2(d − 1) + d 2 − 4d + 8 = . 3d + 2 √ It can be verified that this reduces to the known values 1/ 5 and 1/2 for d = 1 and pt (d) µPoisson d = 2 respectively [3], but also yields the exact expression pt µPoisson (3) √ 4+ 5 = ≈ 0.5669, 11 which as far as we know has not appeared in the literature in closed form. Similarly, Theorem 3.3.1 provides a general formula for block GS-LEX applied to the d-dimensional Poisson problem. If each block is a k-dimensional subproblem, then we have block µPoisson (d, k) = µ pt Poisson (d − k) d−k 2+d−k see Table 5.1. 40 if d − k < 3, otherwise; Figure 5.1: Theoretical smoothing factors of GS-LEX vs. experimental asymptotic convergence rates of two-level multigrid with GS-LEX and GS-RB smoothers, applied to the steady-state diffusion equation −uxx − uyy − εuzz = f . Pointwise relaxation appears on the left, while x-oriented line relaxation is depicted on the right. The discretization is on a 23 × 23 × 23 grid with centered differences. 5.2 The Anisotropic Steady-State Diffusion Problem Suppose the anisotropic steady-state diffusion problem −ρ1 ux1 x1 − ρ2 ux2 x2 − . . . − ρd uxd xd = f (5.1) is discretized over Rd with standard second-order centered differences. We then − obtain a linear system with L h of the form (1.2) where b+ k = bk = ρk , a = 2 ∑dk=1 ρk . In particular, if d = 3 and ρ1 = ρ2 = 1 while ρ3 = ε ∈ (0, 1], then µ pt √ 4 + 5ε 2 − 4ε + 4 = . 6 + 5ε Note that µ pt → 1 as ε → 0. If line relaxation is performed in which unknowns along the x or y directions are relaxed simultaneously, one obtains µ line √ 2 + 5ε 2 − 2ε + 1 = . 3 + 5ε 41 ε µ2 V (1, 1) 1.00 0.25 0.12 0.50 0.32 0.27 0.10 0.70 0.68 0.05 0.83 0.81 Table 5.2: Comparison of predicted smoothing factors of GS-LEX with measured convergence rates of multigrid V (1, 1)-cycles, for the PDE −εuxx − uyy = f . The finest grid is 1023 × 1023, while the coarsest is 1 × 1. Various values of ε are considered. It is interesting to note that µ line above is identical to µ pt in the lower dimensional case d = 2 and ρ1 = 1, ρ2 = ε. Figure 5.1 provides a comparison between the formulas for µ pt and µ line and measured asymptotic converge rates of two-level multigrid on a 23 × 23 × 23 grid, for both GS-LEX and GS-RB. The graphs illustrate that GS-RB is superior to GS-LEX in this symmetric case. We also consider multigrid V -cycles with one pre and one post GS-LEX smoothing step, applied to model problem (5.1) with d = 2, ρ1 = 1, ρ2 = ε on a 1023 × 1023 grid. These are compared against theoretical smoothing factors in Table 5.2. 5.3 Alternative Coarsening for the 2D Poisson Problem Suppose we wish to solve the two-dimensional Poisson problem uxx + uyy = f , (discretized using centered differences) by way of multigrid with a coarsening strategy H = kh, where k = 2, 3, 4 . . . It follows from Equation (4.1) that the smoothing factor of point GS-LEX in this case is given by µ pt = 1 2 − cos π k . Increasing k decreases the work per iteration by making the coarse grids diminish in size more rapidly, however, since µ pt → 1 as k → ∞, the number of iterations is at the same time increasing without bound. The analysis of Chapter 4 predicts that the standard choice k = 2 balances these 42 Smoothing factor V (1, 0)-cycle convergence rate Clock time per grid point (microseconds) Relative work (measured) Relative work (theoretical) H = 2h 0.5 0.43 74 1 1 H = 3h 0.66 0.64 125 1.7 1.4 H = 4h 0.77 0.76 188 2.5 2.2 H = 5h 0.84 0.84 295 3.9 3.1 Table 5.3: Predicted smoothing factors smoothing vs. measured convergence rates of multigrid with H = kh coarsening and GS-LEX smoothing, applied to the 2D Poisson problem discretized using centered differences, for various values of k. V (1, 0)-cycles are used, with a finest grid of size kNk − 1 in each direction, and a coarsest grid consisting of a single point. Nk = logk (999) is chosen to make the total number of grid points approximately one million. Also included is the total work as measured by the running time of the program divided by the size of the problem this is given both in microseconds per grid point and as a ratio relative to the reference case H = 2h. This latter form is compared against the theoretical relative work (4.2). opposing forces in such a way as to minimize their product, the total work1 . To test this prediction, numerical experiments with H = kh for k = 2, 3, 4, 5, 6 have been performed. The overall work is measured in terms of clock time on the computer running the experiment, however, this value is not used directly. Standard multigrid prolongation and restriction operators require that the resolution of the finest grid be (kn − 1) × (kn − 1) for some n ∈ N. As a result it is not possible run all experiments on the same grid; for this reason the total work is calculated by dividing the running time of the program by the number of grid points on the finest grid. For each k, a grid of size (kNk − 1) × (kNk − 1) is used, where Nk is chosen to be the largest integer such that the number of points on the finest grid is less than 106 . Thus N2 = 9, N3 = 6, N4 = N5 = 4, N6 = 3. The results are shown in Table 5.3. The first two rows compare predicted smoothing factors with measured V (1, 0)-cycle convergence rates, showing close agreement. The next three rows compare the total work of V (1, 0)-cycles with the theoretical work (4.2). 1 To be precise, H = 2h is optimal amongst the choices H = kh where k ∈ N. In fact, H = 1.75h is in theory slightly more efficient. 43 H = 6h 0.88 0.88 424 5.7 4.3 It is found that k = 2 is the most efficient of these possibilities. However, the agreement between the measured and theoretical work is not as close as the agreement between smoothing factors and measured convergence rates. Experimental results show that k > 2 is somewhat less efficient than expected. One possible explanation for this discrepancy lies in the observation that M ATLAB (the platform on which the experiments were run) has highly optimized built in functions and vectorization, making estimates of the form (4.2) less meaningful. 5.4 The Convection-Diffusion Problem Suppose the constant coefficient convection-diffusion problem −∆u + w · ∇u = f is discretized on a uniform rectangular grid with mesh spacing h. We define the mesh Reynolds numbers by γk = wk h/2, and assume they are all equal, that is γ1 = γ2 = . . . = γd = γ. If we discretize this PDE using centered differences, we obtain a = 2d, b+ k = 1 − γ, b− k = 1 + γ for all k. Provided |γ| < 1, (1.3) is satisfied and Theorem 3.1.5 gives µ pt = (1 − |γ|) (d − 1) + 1 + (1 + |γ|)2 (1 − d2 )2 2d + (1 + |γ|)2 (1 − d2 ) . Note that the condition |γ| < 1 means we are not in the convection-dominated regime. On the other hand, if we use first order upwinding, then (assuming for − simplicity that γ ≥ 0) we obtain a = 2d(1 + γ), b+ k = 1, bk = 1 + 2γ for all k. The constraints (1.3) are satisfied for all γ ≥ 0 (in other words, the convectiondominated case is included) and µ pt 2 γ (d − 1) + 1 + 1 + 1+γ 1 − d2 1 = 2 1+γ γ 2d + 1 + 1+γ 1 − d2 44 2 . centered differences, γ ∈ (−1, 1) √ d=1 upwinding, γ ≥ 0 1−|γ| 1 8γ 2 +12γ+5 1 2(1+γ) √ 4+(1+|γ|)2 1−|γ| 2 d=2 (1−|γ|) 2+ 1+ d=3 1+|γ| 2 2 γ 2+ 1+ 41 1+ 1+γ 6− 21 (1+|γ|)2 γ (1+γ) 6− 21 1+ 1+γ 2 2 Table 5.4: Smoothing factors of pointwise GS-LEX applied to two common discretizations of the convection-diffusion equation with all mesh Reynolds numbers equal to γ, for d = 1, 2, 3. Table 5.4 lists the above smoothing factors for the cases of interest d = 1, 2, 3. For asymmetric problems such as convection-diffusion, Galerkin coarse grid operators obtained by linear interpolation and full weighting may cause unstable coarse grid discretizations (even in diffusion dominated settings), leading to multigrid convergence problems, especially as the number of levels increases [14]. Indeed, in Table 5.5 we observe divergence for multigrid with three or more levels for upwinding in the convection-dominated regime; for two levels, however, we see a good agreement between measured convergence rates and the predicted smoothing factor. We note that there are ways of dealing with these instabilities within the Galerkin framework; see for example [5]. γ µ 2L 3L 4L 0.00 0.50 0.39 0.41 0.41 0.33 0.38 0.36 0.36 NC 0.66 0.30 0.30 0.31 NC 1.00 0.25 0.25 NC NC 3.00 0.13 0.13 NC NC 5.00 0.08 0.08 NC NC Table 5.5: Comparison of the smoothing factor µ of GS-LEX with the measured asymptotic convergence rate of multigrid applied to the PDE −uxx − uyy + σ ux + σ uy = f , discretized uniformly. The problem is discretized using upwinding on a 1023 × 1023 grid. The Galerkin coarse grid operator is used. W (1, 0)-cycles are used, and various values of γ = σ h/2, where h is the grid spacing on the finest mesh, are considered. ‘NC’ stands for no convergence. ‘nL’ (with n = 2, 3, 4) signifies the number of levels. 45 The coarse grid operator L H may also be defined by direct discretization of the underlying PDE on the coarse mesh. For upwinding with constant coefficients, this yields a stable discretization on all grid levels, and we can expect multigrid methods to converge. However, it has been shown that in this setting, for strong convection the determining factor of the multigrid convergence rate may not be the smoother, but rather the coarse grid correction [4]. In particular, in the high convection limit, there are certain smooth modes which are reduced by a constant factor that grows as the number of levels increases. In the convection-dominated regime where the smoothing factor µ becomes smaller than this factor, it is these smooth modes rather than the rough modes with reduction factor µ, that dominate the error asymptotically. Therefore, we expect that the asymptotic convergence of multigrid will follow the trend indicated by our smoothing analysis in the diffusiondominated regime, while tending to a constant in the convection-dominated regime; see Table 5.6. We also include GS-RB W -cycles in the table, to illustrate the superiority of GS-LEX in this case. γ µ 2L V (1, 0) W (1, 0)LEX W (1, 0)RB 0.0 0.50 0.37 0.42 0.37 0.23 1.0 0.25 0.25 0.56 0.29 0.37 2.0 0.17 0.24 0.54 0.34 0.46 3.0 0.13 0.27 0.52 0.35 0.59 4.0 0.10 0.29 0.50 0.36 0.71 5.0 0.08 0.30 0.48 0.36 0.79 Table 5.6: Comparison of the smoothing factor µ of GS-LEX with measured convergence rates of multigrid V (1, 0) and W (1, 0)-cycles, for the PDE −uxx − uyy + σ (ux + uy ) = f , discretized uniformly. The finest grid is 1023 × 1023, while the coarsest is 1 × 1. The coarse grid operator is obtained by direct discretization of the PDE on the coarse mesh. Various values of γ = σ h/2 are considered, where h is the grid spacing of the finest mesh. ‘2L’ stands for a two-level scheme, whereas V (1, 0) and W (1, 0) signify V and W -cycles, respectively, with one pre-smoothing sweep and no post-smoothing. For the centered difference discretization of the convection-diffusion problem, our smoothing analysis applies to the case where the mesh Reynolds number γ is less than one in magnitude. But in a multigrid setup, we have different mesh 46 γ µ 2L 3L 4L 5L 0.00 0.50 0.37 0.37 0.37 0.37 0.20 0.40 0.36 0.36 0.36 NC 0.40 0.30 0.30 0.30 NC NC 0.60 0.20 0.21 NC NC NC 0.80 0.10 0.10 NC NC NC 0.99 0.005 0.005 NC NC NC Table 5.7: Comparison of the smoothing factor µ of GS-LEX with the measured asymptotic convergence rate of two-level and multi-level multigrid applied to the PDE −uxx − uyy + σ (ux + uy ) = f , discretized uniformly. The problem is discretized using centered differences on a 1023 × 1023 grid. The coarse grid operator is obtained by direct discretization of the PDE on the coarse mesh. W (1, 0)-cycles are used, and various values of γ = σ h/2, where h is the grid spacing on the finest mesh, are considered. ‘NC’ stands for no convergence. Reynolds numbers on different grids - for our analysis to apply, they must all obey this constraint. Therefore, if L denotes the number of levels in our algorithm and γ denotes the mesh Reynolds number on the finest grid, what we really require is |γ| < 22−L . In Table 5.7, we compare our predictions with measured convergence rates for various values of L; W (1, 0)-cycles are applied. In practice, we can get away with a mesh Reynolds number moderately larger than 1 on the coarsest grid and still maintain expected convergence rates. 5.5 The Time-dependent Diffusion Problem Suppose the time-dependent diffusion problem ut = ρ1 ux1 x1 + ρ2 ux2 x2 + . . . + ρd uxd xd is discretized on a uniform mesh using backward Euler in time and centered differences in space, with multigrid used to solve the resulting linear system at each time step. If τ denotes the timestep size, then µ pt follows from Theorem 3.1.5 by setting cr = ρ − ρm , 2ρ + h2 /τ cm = 47 ρm , 2ρ + h2 /τ α = 1. where ρ = ∑dk=1 ρk and ρm = min({ρk }). In particular, in the isotropic case ρ1 = ρ2 = . . . = ρd = 1, we obtain µ pt = (d − 1)(2d + h2 /τ) + (2d + h2 /τ)2 + d 2 (2 − d)2 . (2d + h2 /τ)2 + d(2 − d) When d = 2 this reduces to the particularly simple expression µ pt = 1 2 h 2 + 2τ . Note the (expected) behavior of the smoothing factor as a function of τ: when h is fixed and τ gets smaller, the smoothing factor becomes smaller too. 48 Chapter 6 Conclusions Using results from complex analysis, primarily the maximum modulus principle and properties of M¨obius transformations, we have derived closed-form expressions for the smoothing factors of lexicographic pointwise and block Gauss-Seidel (Theorems 3.1.5 and 3.3.1). An extension of Theorem 3.1.5 that incorporates a relaxation parameter is provided in (3.10). In the pointwise case our results are valid for general operators of the form (1.2) satisfying the constraints (1.3) and (3.2), whereas in the block case we require the additional assumption of symmetry. In some cases, block smoothing on a high dimensional problem has the same smoothing factor as pointwise smoothing on a related lower dimensional problem. Our analysis provides smoothing factors for, among other equations, the following: • Pointwise GS-LEX applied to the d-dimensional anisotropic steady-state diffusion equation − ∑dk=1 ρk uxk xk = f , discretized with centered differences. • Pointwise and block GS-LEX applied to the d-dimensional Poisson equation −∆u = f , discretized with centered differences. • Pointwise GS-LEX applied to the d-dimensional constant coefficient convectiondiffusion equation −∆u + w · ∇u = f , discretized with centered differences or upwinding, and all mesh Reynolds numbers equal. • Pointwise GS-LEX applied to the linear systems arising in each timestep of 49 the solution of the d-dimensional time-dependent diffusion equation ut = ∑dk=1 ρk uxk xk , discretized with centered differences in space and backward Euler in time. For pointwise GS-LEX, we have generalized our results to allow mesh ratios besides H = 2h. We have shown that for high dimensional Poisson problems, a judicious choice of 2h > H > h can in theory lead to a reduction of up to 50% in the overall work. We have also observed that lexicographic Gauss-Seidel smoothing is effective for equations with strong asymmetry. In particular, for the constant coefficient convection-diffusion equation with equal mesh Reynolds numbers we have shown for upwind discretizations that GS-LEX has a smaller smoothing factor than GSRB in the convection-dominated regime. 50 Bibliography [1] E. L. Akim, O. M. Belotserkovskii, K. V. Brushlinskii, S. K. Godunov, V. F. Dyachenko, V. T. Zhukov, M. K. Kerimov, A. I. Lobanov, M. V. Maslennikov, I. B. Petrov, Y. P. Popov, G. P. Prokopov, V. S. Ryabenkii, L. G. Strakhovskaya, A. S. Kholodov, B. N. Chetverushkin, and T. M. Eneev. In Memory of Radii Petrovich Fedorenko. Computational Mathematics and Mathematical Physics, 50:1459–1463, 2010. [2] U. M. Ascher and C. Greif. A First Course in Numerical Methods. Computational Science and Engineering. SIAM, Philadelphia, 2011. [3] A. Brandt. Multi-Level Adaptive Solutions to Boundary-Value Problems. Mathematics of Computation, 31(138):333–390, 1977. [4] A. Brandt and I. Yavneh. Accelerated Multigrid Convergence and High-Reynolds Recirculating Flows. SIAM J. Sci. Comput., 14(3):607–626, 1993. [5] P. de Zeeuw. Matrix-dependent Prolongations and Restrictions in a Blackbox Multigrid Solver. Journal of Computational and Applied Mathematics, 33(1):1–27, 1990. [6] J. Dendy. Black Box Multigrid with Coarsening by a Factor of Three. Numerical Linear Algebra with Applications, 17:577–598, 2010. [7] H. C. Elman and M. P. Chernesky. Ordering Effects on Relaxation Methods Applied to the Discrete One- Dimensional Convection-Diffusion Equation. SIAM Journal on Numerical Analysis, 30(5):1268–1290, 1993. [8] R. P. Fedorenko. A relaxation method for solving elliptic difference equations. USSR Computational Mathematics and Mathematical Physics, 1 (4):1092–1096, 1962. 51 [9] R. P. Fedorenko. The speed of convergence of one iterative process. USSR Computational Mathematics and Mathematical Physics, 4(3):227–235, 1964. [10] T. W. Gamelin. Complex Analysis. Springer, 2001. [11] W. Hackbusch. Multigrid methods and applications, volume 4 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1985. [12] H. Nyquist. Certain Topics in Telegraph Transmission Theory. Trans. Amer. Inst. Electr. Eng., pages 617–644, 1928. [13] U. Trottenberg, C. W. Oosterlee, and A. Sch¨uller. Multigrid. Academic Press Inc., San Diego, CA, 2001. [14] E. J. van Asselt and P. de Zeeuw. The Convergence Rate of Multi-Level Algorithms Applied to the Convection-Diffusion Equation. SIAM Journal on Scientific and Statistical Computing, 6(2):492–503, 1985. [15] R. S. Varga. Matrix Iterative Analysis. Prentice-Hall Inc., Englewood Cliffs, N.J., 1962. [16] P. Wesseling. An Introduction to Multigrid Methods. Pure and Applied Mathematics (New York). John Wiley & Sons Ltd., Chichester, 1992. [17] I. Yavneh. Multigrid Smoothing Factors of Red-Black Gauss-Seidel Applied to a Class of Elliptic Operators. SIAM J. Numer. Anal., 32(6):1126–1138, 1995. [18] I. Yavneh. On Red-Black SOR Smoothing in Multigrid. SIAM J. Sci. Comput., 17(1):180–192, 1996. [19] I. Yavneh. Why Multigrid Methods are so Efficient. Computing in Science and Engineering, 8(6):12–22, 2006. [20] H. B. Zubair, C. Oosterlee, and R. Wienands. Multigrid for High-Dimensional Elliptic Partial Differential Equations on Non-equidistant Grids. SIAM J. Sci. Comput., 29(4):1613–1636, 2007. 52
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A complex analysis based derivation of multigrid smoothing...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A complex analysis based derivation of multigrid smoothing factors of lexicographic Gauss-Seidel Hocking, Laird Robert 2011
pdf
Page Metadata
Item Metadata
Title | A complex analysis based derivation of multigrid smoothing factors of lexicographic Gauss-Seidel |
Creator |
Hocking, Laird Robert |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | This thesis aims to present a unified framework for deriving analytical formulas for multigrid smoothing factors in arbitrary dimensions, under certain simplifying assumptions. To derive these expressions we rely on complex analysis and geometric considerations, using the maximum modulus principle and Mobius transformations. We restrict our attention to pointwise and block lexicographic Gauss-Seidel smoothers on a d-dimensional uniform mesh, where the computational molecule of the associated discrete operator forms a 2d+1 point star. In the pointwise case the effect of a relaxation parameter, as well as different choices of mesh ratio, are analyzed. The results apply to any number of spatial dimensions, and are applicable to high-dimensional versions of a few common model problems with constant coefficients, including the Poisson and anisotropic diffusion equations as well as the convection-diffusion equation. We show that in most cases our formulas, exact under the simplifying assumptions of Local Fourier Analysis, form tight upper bounds for the asymptotic convergence of geometric multigrid in practice. We also show that there are asymmetric cases where lexicographic Gauss-Seidel smoothing outperforms red-black Gauss-Seidel smoothing; this occurs for certain model convection-diffusion equations with high mesh Reynolds numbers. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-08-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052104 |
URI | http://hdl.handle.net/2429/37012 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2011_fall_hocking_laird.pdf [ 702.48kB ]
- Metadata
- JSON: 24-1.0052104.json
- JSON-LD: 24-1.0052104-ld.json
- RDF/XML (Pretty): 24-1.0052104-rdf.xml
- RDF/JSON: 24-1.0052104-rdf.json
- Turtle: 24-1.0052104-turtle.txt
- N-Triples: 24-1.0052104-rdf-ntriples.txt
- Original Record: 24-1.0052104-source.json
- Full Text
- 24-1.0052104-fulltext.txt
- Citation
- 24-1.0052104.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0052104/manifest