Multigrid method

Last updated
Multigrid method
Class Differential equation

In numerical analysis, a multigrid method (MG method) is an algorithm for solving differential equations using a hierarchy of discretizations. They are an example of a class of techniques called multiresolution methods, very useful in problems exhibiting multiple scales of behavior. For example, many basic relaxation methods exhibit different rates of convergence for short- and long-wavelength components, suggesting these different scales be treated differently, as in a Fourier analysis approach to multigrid. [1] MG methods can be used as solvers as well as preconditioners.

Contents

The main idea of multigrid is to accelerate the convergence of a basic iterative method (known as relaxation, which generally reduces short-wavelength error) by a global correction of the fine grid solution approximation from time to time, accomplished by solving a coarse problem. The coarse problem, while cheaper to solve, is similar to the fine grid problem in that it also has short- and long-wavelength errors. It can also be solved by a combination of relaxation and appeal to still coarser grids. This recursive process is repeated until a grid is reached where the cost of direct solution there is negligible compared to the cost of one relaxation sweep on the fine grid. This multigrid cycle typically reduces all error components by a fixed amount bounded well below one, independent of the fine grid mesh size. The typical application for multigrid is in the numerical solution of elliptic partial differential equations in two or more dimensions. [2]

Multigrid methods can be applied in combination with any of the common discretization techniques. For example, the finite element method may be recast as a multigrid method. [3] In these cases, multigrid methods are among the fastest solution techniques known today. In contrast to other methods, multigrid methods are general in that they can treat arbitrary regions and boundary conditions. They do not depend on the separability of the equations or other special properties of the equation. They have also been widely used for more-complicated non-symmetric and nonlinear systems of equations, like the Lamé equations of elasticity or the Navier-Stokes equations. [4]

Algorithm

Visualization of iterative Multigrid algorithm for fast O(n) convergence. Multigrid Visualization.png
Visualization of iterative Multigrid algorithm for fast O(n) convergence.

There are many variations of multigrid algorithms, but the common features are that a hierarchy of discretizations (grids) is considered. The important steps are: [5] [6]

There are many choices of multigrid methods with varying trade-offs between speed of solving a single iteration and the rate of convergence with said iteration. The 3 main types are V-Cycle, F-Cycle, and W-Cycle. These differ in which and how many coarse-grain cycles are performed per fine iteration. The V-Cycle algorithm executes one coarse-grain V-Cycle. F-Cycle does a coarse-grain V-Cycle followed by a coarse-grain F-Cycle, while each W-Cycle performs two coarse-grain W-Cycles per iteration. For a discrete 2D problem, F-Cycle takes 83% more time to compute than a V-Cycle iteration while a W-Cycle iteration takes 125% more. If the problem is set up in a 3D domain, then a F-Cycle iteration and a W-Cycle iteration take about 64% and 75% more time respectively than a V-Cycle iteration ignoring overheads. Typically, W-Cycle produces similar convergence to F-Cycle. However, in cases of convection-diffusion problems with high Péclet numbers, W-Cycle can show superiority in its rate of convergence per iteration over F-Cycle. The choice of smoothing operators are extremely diverse as they include Krylov subspace methods and can be preconditioned.

Any geometric multigrid cycle iteration is performed on a hierarchy of grids and hence it can be coded using recursion. Since the function calls itself with smaller sized (coarser) parameters, the coarsest grid is where the recursion stops. In cases where the system has a high condition number, the correction procedure is modified such that only a fraction of the prolongated coarser grid solution is added onto the finer grid.

These steps can be used as shown in the MATLAB style pseudo code for 1 iteration of V-Cycle Multigrid:

functionphi=V_Cycle(phi,f,h)% Recursive V-Cycle Multigrid for solving the Poisson equation (\nabla^2 phi = f) on a uniform grid of spacing h% Pre-Smoothingphi=smoothing(phi,f,h);% Compute Residual Errorsr=residual(phi,f,h);% Restrictionrhs=restriction(r);eps=zeros(size(rhs));% stop recursion at smallest grid size, otherwise continue recursionifsmallest_grid_size_is_achievedeps=coarse_level_solve(eps,rhs,2*h);elseeps=V_Cycle(eps,rhs,2*h);end% Prolongation and Correctionphi=phi+prolongation(eps);% Post-Smoothingphi=smoothing(phi,f,h);end

The following represents F-cycle multigrid. This multigrid cycle is slower than V-Cycle per iteration but does result in faster convergence.

functionphi=F_Cycle(phi,f,h)% Recursive F-cycle multigrid for solving the Poisson equation (\nabla^2 phi = f) on a uniform grid of spacing h% Pre-smoothingphi=smoothing(phi,f,h);% Compute Residual Errorsr=residual(phi,f,h);% Restrictionrhs=restriction(r);eps=zeros(size(rhs));% stop recursion at smallest grid size, otherwise continue recursionifsmallest_grid_size_is_achievedeps=coarse_level_solve(eps,rhs,2*h);elseeps=F_Cycle(eps,rhs,2*h);end% Prolongation and Correctionphi=phi+prolongation(eps);% Re-smoothingphi=smoothing(phi,f,h);% Compute residual errorsr=residual(phi,f,h);% Restrictionrhs=restriction(r);% stop recursion at smallest grid size, otherwise continue recursionifsmallest_grid_size_is_achievedeps=coarse_level_solve(eps,rhs,2*h);elseeps=V_Cycle(eps,rhs,2*h);end% Prolongation and Correctionphi=phi+prolongation(eps);% Post-smoothingphi=smoothing(phi,f,h);end

Similarly the procedures can modified as shown in the MATLAB style pseudo code for 1 iteration of W-cycle multigrid for an even superior rate of convergence in certain cases:

functionphi=W_cycle(phi,f,h)% Recursive W-cycle multigrid for solving the Poisson equation (\nabla^2 phi = f) on a uniform grid of spacing h% Pre-smoothingphi=smoothing(phi,f,h);% Compute Residual Errorsr=residual(phi,f,h);% Restrictionrhs=restriction(r);eps=zeros(size(rhs));% stop recursion at smallest grid size, otherwise continue recursionifsmallest_grid_size_is_achievedeps=coarse_level_solve(eps,rhs,2*h);elseeps=W_cycle(eps,rhs,2*h);end% Prolongation and correctionphi=phi+prolongation(eps);% Re-smoothingphi=smoothing(phi,f,h);% Compute residual errorsr=residual(phi,f,h);% Restrictionrhs=restriction(r);% stop recursion at smallest grid size, otherwise continue recursionifsmallest_grid_size_is_achievedeps=coarse_level_solve(eps,rhs,2*h);elseeps=W_cycle(eps,rhs,2*h);end% Prolongation and correctionphi=phi+prolongation(eps);% Post-smoothingphi=smoothing(phi,f,h);end

Computational cost [ citation needed ]

Assuming a 2-dimensional problem setup, the computation moves across grid hierarchy differently for various multigrid cycles. MultigridWork.svg
Assuming a 2-dimensional problem setup, the computation moves across grid hierarchy differently for various multigrid cycles.

This approach has the advantage over other methods that it often scales linearly with the number of discrete nodes used. In other words, it can solve these problems to a given accuracy in a number of operations that is proportional to the number of unknowns.

Assume that one has a differential equation which can be solved approximately (with a given accuracy) on a grid with a given grid point density . Assume furthermore that a solution on any grid may be obtained with a given effort from a solution on a coarser grid . Here, is the ratio of grid points on "neighboring" grids and is assumed to be constant throughout the grid hierarchy, and is some constant modeling the effort of computing the result for one grid point.

The following recurrence relation is then obtained for the effort of obtaining the solution on grid :

Example of Convergence Rates of Multigrid Cycles in comparison to other smoothing operators. Convergence Rate of Multigrid Cycles.svg
Example of Convergence Rates of Multigrid Cycles in comparison to other smoothing operators.

And in particular, we find for the finest grid that

Combining these two expressions (and using ) gives

Using the geometric series, we then find (for finite )

that is, a solution may be obtained in time. It should be mentioned that there is one exception to the i.e. W-cycle multigrid used on a 1D problem; it would result in complexity. [ citation needed ]

Multigrid preconditioning

A multigrid method with an intentionally reduced tolerance can be used as an efficient preconditioner for an external iterative solver, e.g., [7] The solution may still be obtained in time as well as in the case where the multigrid method is used as a solver. Multigrid preconditioning is used in practice even for linear systems, typically with one cycle per iteration, e.g., in Hypre. Its main advantage versus a purely multigrid solver is particularly clear for nonlinear problems, e.g., eigenvalue problems.

If the matrix of the original equation or an eigenvalue problem is symmetric positive definite (SPD), the preconditioner is commonly constructed to be SPD as well, so that the standard conjugate gradient (CG) iterative methods can still be used. Such imposed SPD constraints may complicate the construction of the preconditioner, e.g., requiring coordinated pre- and post-smoothing. However, preconditioned steepest descent and flexible CG methods for SPD linear systems and LOBPCG for symmetric eigenvalue problems are all shown [8] to be robust if the preconditioner is not SPD.

Bramble–Pasciak–Xu preconditioner

Originally described in Xu's Ph.D. thesis [9] and later published in Bramble-Pasciak-Xu, [10] the BPX-preconditioner is one of the two major multigrid approaches (the other is the classic multigrid algorithm such as V-cycle) for solving large-scale algebraic systems that arise from the discretization of models in science and engineering described by partial differential equations. In view of the subspace correction framework, [11] BPX preconditioner is a parallel subspace correction method where as the classic V-cycle is a successive subspace correction method. The BPX-preconditioner is known to be naturally more parallel and in some applications more robust than the classic V-cycle multigrid method. The method has been widely used by researchers and practitioners since 1990.

Generalized multigrid methods

Multigrid methods can be generalized in many different ways. They can be applied naturally in a time-stepping solution of parabolic partial differential equations, or they can be applied directly to time-dependent partial differential equations. [12] Research on multilevel techniques for hyperbolic partial differential equations is underway. [13] Multigrid methods can also be applied to integral equations, or for problems in statistical physics. [14]

Another set of multiresolution methods is based upon wavelets. These wavelet methods can be combined with multigrid methods. [15] [16] For example, one use of wavelets is to reformulate the finite element approach in terms of a multilevel method. [17]

Adaptive multigrid exhibits adaptive mesh refinement, that is, it adjusts the grid as the computation proceeds, in a manner dependent upon the computation itself. [18] The idea is to increase resolution of the grid only in regions of the solution where it is needed.

Algebraic multigrid (AMG)

Practically important extensions of multigrid methods include techniques where no partial differential equation nor geometrical problem background is used to construct the multilevel hierarchy. [19] Such algebraic multigrid methods (AMG) construct their hierarchy of operators directly from the system matrix. In classical AMG, the levels of the hierarchy are simply subsets of unknowns without any geometric interpretation. (More generally, coarse grid unknowns can be particular linear combinations of fine grid unknowns.) Thus, AMG methods become black-box solvers for certain classes of sparse matrices. AMG is regarded as advantageous mainly where geometric multigrid is too difficult to apply, [20] but is often used simply because it avoids the coding necessary for a true multigrid implementation. While classical AMG was developed first, a related algebraic method is known as smoothed aggregation (SA).

In an overview paper [21] by Jinchao Xu and Ludmil Zikatanov, the "algebraic multigrid" methods are understood from an abstract point of view. They developed a unified framework and existing algebraic multigrid methods can be derived coherently. Abstract theory about how to construct optimal coarse space as well as quasi-optimal spaces was derived. Also, they proved that, under appropriate assumptions, the abstract two-level AMG method converges uniformly with respect to the size of the linear system, the coefficient variation, and the anisotropy. Their abstract framework covers most existing AMG methods, such as classical AMG, energy-minimization AMG, unsmoothed and smoothed aggregation AMG, and spectral AMGe.

Multigrid in time methods

Multigrid methods have also been adopted for the solution of initial value problems. [22] Of particular interest here are parallel-in-time multigrid methods: [23] in contrast to classical Runge–Kutta or linear multistep methods, they can offer concurrency in temporal direction. The well known Parareal parallel-in-time integration method can also be reformulated as a two-level multigrid in time.

Multigrid for nearly singular problems

Nearly singular problems arise in a number of important physical and engineering applications. Simple, but important example of nearly singular problems can be found at the displacement formulation of linear elasticity for nearly incompressible materials. Typically, the major problem to solve such nearly singular systems boils down to treat the nearly singular operator given by robustly with respect to the positive, but small parameter . Here is symmetric semidefinite operator with large null space, while is a symmetric positive definite operator. There were many works to attempt to design a robust and fast multigrid method for such nearly singular problems. A general guide has been provided as a design principle to achieve parameters (e.g., mesh size and physical parameters such as Poisson's ratio that appear in the nearly singular operator) independent convergence rate of the multigrid method applied to such nearly singular systems, [24] i.e., in each grid, a space decomposition based on which the smoothing is applied, has to be constructed so that the null space of the singular part of the nearly singular operator has to be included in the sum of the local null spaces, the intersection of the null space and the local spaces resulting from the space decompositions.

Notes

  1. Roman Wienands; Wolfgang Joppich (2005). Practical Fourier analysis for multigrid methods. CRC Press. p. 17. ISBN   978-1-58488-492-7.
  2. U. Trottenberg; C. W. Oosterlee; A. Schüller (2001). Multigrid. Academic Press. ISBN   978-0-12-701070-0.
  3. Yu Zhu; Andreas C. Cangellaris (2006). Multigrid finite element methods for electromagnetic field modeling. Wiley. p. 132 ff. ISBN   978-0-471-74110-7.
  4. Shah, Tasneem Mohammad (1989). Analysis of the multigrid method (Thesis). Oxford University. Bibcode:1989STIN...9123418S.
  5. M. T. Heath (2002). "Section 11.5.7 Multigrid Methods". Scientific Computing: An Introductory Survey. McGraw-Hill Higher Education. p. 478 ff. ISBN   978-0-07-112229-0.
  6. P. Wesseling (1992). An Introduction to Multigrid Methods. Wiley. ISBN   978-0-471-93083-9.
  7. Andrew V Knyazev, Klaus Neymeyr. Efficient solution of symmetric eigenvalue problems using multigrid preconditioners in the locally optimal block conjugate gradient method. Electronic Transactions on Numerical Analysis, 15, 38–55, 2003.
  8. Bouwmeester, Henricus; Dougherty, Andrew; Knyazev, Andrew V. (2015). "Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods 1". Procedia Computer Science. 51: 276–285. arXiv: 1212.6680 . doi: 10.1016/j.procs.2015.05.241 . S2CID   51978658.
  9. Xu, Jinchao. Theory of multilevel methods. Vol. 8924558. Ithaca, NY: Cornell University, 1989.
  10. Bramble, James H., Joseph E. Pasciak, and Jinchao Xu. "Parallel multilevel preconditioners." Mathematics of Computation 55, no. 191 (1990): 1–22.
  11. Xu, Jinchao. "Iterative methods by space decomposition and subspace correction." SIAM review 34, no. 4 (1992): 581-613.
  12. F. Hülsemann; M. Kowarschik; M. Mohr; U. Rüde (2006). "Parallel geometric multigrid". In Are Magnus Bruaset; Aslak Tveito (eds.). Numerical solution of partial differential equations on parallel computers. Birkhäuser. p. 165. ISBN   978-3-540-29076-6.
  13. For example, J. Blaz̆ek (2001). Computational fluid dynamics: principles and applications. Elsevier. p. 305. ISBN   978-0-08-043009-6. and Achi Brandt and Rima Gandlin (2003). "Multigrid for Atmospheric Data Assimilation: Analysis". In Thomas Y. Hou; Eitan Tadmor (eds.). Hyperbolic problems: theory, numerics, applications: proceedings of the Ninth International Conference on Hyperbolic Problems of 2002. Springer. p. 369. ISBN   978-3-540-44333-9.
  14. Achi Brandt (2002). "Multiscale scientific computation: review". In Timothy J. Barth; Tony Chan; Robert Haimes (eds.). Multiscale and multiresolution methods: theory and applications. Springer. p. 53. ISBN   978-3-540-42420-8.
  15. Björn Engquist; Olof Runborg (2002). "Wavelet-based numerical homogenization with applications". In Timothy J. Barth; Tony Chan; Robert Haimes (eds.). Multiscale and Multiresolution Methods. Vol. 20 of Lecture Notes in Computational Science and Engineering. Springer. p. 140 ff. ISBN   978-3-540-42420-8.
  16. U. Trottenberg; C. W. Oosterlee; A. Schüller (2001). Multigrid. Academic Press. ISBN   978-0-12-701070-0.
  17. Albert Cohen (2003). Numerical Analysis of Wavelet Methods. Elsevier. p. 44. ISBN   978-0-444-51124-9.
  18. U. Trottenberg; C. W. Oosterlee; A. Schüller (2001). "Chapter 9: Adaptive Multigrid". Multigrid. Academic Press. p. 356. ISBN   978-0-12-701070-0.
  19. Yair Shapira (2003). "Algebraic multigrid". Matrix-based multigrid: theory and applications. Springer. p. 66. ISBN   978-1-4020-7485-1.
  20. U. Trottenberg; C. W. Oosterlee; A. Schüller (2001). Multigrid. Academic Press. p. 417. ISBN   978-0-12-701070-0.
  21. Xu, J. and Zikatanov, L., 2017. Algebraic multigrid methods. Acta Numerica, 26, pp.591-721.
  22. Hackbusch, Wolfgang (1985). "Parabolic multi-grid methods". Computing Methods in Applied Sciences and Engineering, VI: 189–197. ISBN   9780444875976 . Retrieved 1 August 2015.
  23. Horton, Graham (1992). "The time-parallel multigrid method". Communications in Applied Numerical Methods. 8 (9): 585–595. doi:10.1002/cnm.1630080906.
  24. Young-Ju Lee, Jinbiao Wu, Jinchao Xu and Ludmil Zikatanov, Robust Subspace Correction Methods for Nearly Singular Systems, Mathematical Models and Methods in Applied Sciences, Vol. 17, No 11, pp. 1937-1963 (2007)

Related Research Articles

In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the i-th approximation is derived from the previous ones.

Spectral methods are a class of techniques used in applied mathematics and scientific computing to numerically solve certain differential equations. The idea is to write the solution of the differential equation as a sum of certain "basis functions" and then to choose the coefficients in the sum in order to satisfy the differential equation as well as possible.

<span class="mw-page-title-main">Computational fluid dynamics</span> Analysis and solving of problems that involve fluid flows

Computational fluid dynamics (CFD) is a branch of fluid mechanics that uses numerical analysis and data structures to analyze and solve problems that involve fluid flows. Computers are used to perform the calculations required to simulate the free-stream flow of the fluid, and the interaction of the fluid with surfaces defined by boundary conditions. With high-speed supercomputers, better solutions can be achieved, and are often required to solve the largest and most complex problems. Ongoing research yields software that improves the accuracy and speed of complex simulation scenarios such as transonic or turbulent flows. Initial validation of such software is typically performed using experimental apparatus such as wind tunnels. In addition, previously performed analytical or empirical analysis of a particular problem can be used for comparison. A final validation is often performed using full-scale testing, such as flight tests.

<span class="mw-page-title-main">Smoothed-particle hydrodynamics</span> Method of hydrodynamics simulation

Smoothed-particle hydrodynamics (SPH) is a computational method used for simulating the mechanics of continuum media, such as solid mechanics and fluid flows. It was developed by Gingold and Monaghan and Lucy in 1977, initially for astrophysical problems. It has been used in many fields of research, including astrophysics, ballistics, volcanology, and oceanography. It is a meshfree Lagrangian method, and the resolution of the method can easily be adjusted with respect to variables such as density.

Numerical methods for partial differential equations is the branch of numerical analysis that studies the numerical solution of partial differential equations (PDEs).

Linear multistep methods are used for the numerical solution of ordinary differential equations. Conceptually, a numerical method starts from an initial point and then takes a short step forward in time to find the next solution point. The process continues with subsequent steps to map out the solution. Single-step methods refer to only one previous point and its derivative to determine the current value. Methods such as Runge–Kutta take some intermediate steps to obtain a higher order method, but then discard all previous information before taking a second step. Multistep methods attempt to gain efficiency by keeping and using the information from previous steps rather than discarding it. Consequently, multistep methods refer to several previous points and derivative values. In the case of linear multistep methods, a linear combination of the previous points and derivative values is used.

In numerical linear algebra, the method of successive over-relaxation (SOR) is a variant of the Gauss–Seidel method for solving a linear system of equations, resulting in faster convergence. A similar method can be used for any slowly converging iterative process.

In mathematics, preconditioning is the application of a transformation, called the preconditioner, that conditions a given problem into a form that is more suitable for numerical solving methods. Preconditioning is typically related to reducing a condition number of the problem. The preconditioned problem is then usually solved by an iterative method.

In differential geometry, a Kähler–Einstein metric on a complex manifold is a Riemannian metric that is both a Kähler metric and an Einstein metric. A manifold is said to be Kähler–Einstein if it admits a Kähler–Einstein metric. The most important special case of these are the Calabi–Yau manifolds, which are Kähler and Ricci-flat.

In numerical mathematics, relaxation methods are iterative methods for solving systems of equations, including nonlinear systems.

In numerical linear algebra, the alternating-direction implicit (ADI) method is an iterative method used to solve Sylvester matrix equations. It is a popular method for solving the large matrix equations that arise in systems theory and control, and can be formulated to construct solutions in a memory-efficient, factored form. It is also used to numerically solve parabolic and elliptic partial differential equations, and is a classic method used for modeling heat conduction and solving the diffusion equation in two or more dimensions. It is an example of an operator splitting method.

In computational fluid dynamics, shock-capturing methods are a class of techniques for computing inviscid flows with shock waves. The computation of flow containing shock waves is an extremely difficult task because such flows result in sharp, discontinuous changes in flow variables such as pressure, temperature, density, and velocity across the shock.

<span class="mw-page-title-main">Domain decomposition methods</span>

In mathematics, numerical analysis, and numerical partial differential equations, domain decomposition methods solve a boundary value problem by splitting it into smaller boundary value problems on subdomains and iterating to coordinate the solution between adjacent subdomains. A coarse problem with one or few unknowns per subdomain is used to further coordinate the solution between the subdomains globally. The problems on the subdomains are independent, which makes domain decomposition methods suitable for parallel computing. Domain decomposition methods are typically used as preconditioners for Krylov space iterative methods, such as the conjugate gradient method, GMRES, and LOBPCG.

In mathematics, Neumann–Neumann methods are domain decomposition preconditioners named so because they solve a Neumann problem on each subdomain on both sides of the interface between the subdomains. Just like all domain decomposition methods, so that the number of iterations does not grow with the number of subdomains, Neumann–Neumann methods require the solution of a coarse problem to provide global communication. The balancing domain decomposition is a Neumann–Neumann method with a special kind of coarse problem.

In numerical analysis, coarse problem is an auxiliary system of equations used in an iterative method for the solution of a given larger system of equations. A coarse problem is basically a version of the same problem at a lower resolution, retaining its essential characteristics, but with fewer variables. The purpose of the coarse problem is to propagate information throughout the whole problem globally.

Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) is a matrix-free method for finding the largest eigenvalues and the corresponding eigenvectors of a symmetric generalized eigenvalue problem

In the mathematical discipline of numerical linear algebra, a matrix splitting is an expression which represents a given matrix as a sum or difference of matrices. Many iterative methods depend upon the direct solution of matrix equations involving matrices more general than tridiagonal matrices. These matrix equations can often be solved directly and efficiently when written as a matrix splitting. The technique was devised by Richard S. Varga in 1960.

In linear algebra, a convergent matrix is a matrix that converges to the zero matrix under matrix exponentiation.

<span class="mw-page-title-main">Parareal</span> Parallel algorithm from numerical analysis

Parareal is a parallel algorithm from numerical analysis and used for the solution of initial value problems. It was introduced in 2001 by Lions, Maday and Turinici. Since then, it has become one of the most widely studied parallel-in-time integration methods.

<span class="mw-page-title-main">Quasilinearization</span> Quasilinearization replaces a nonlinear operator equation with a sequence of linear operator equatio

In mathematics, quasilinearization is a technique which replaces a nonlinear differential equation or operator equation with a sequence of linear problems, which are presumed to be easier, and whose solutions approximate the solution of the original nonlinear problem with increasing accuracy. It is a generalization of Newton's method; the word "quasilinearization" is commonly used when the differential equation is a boundary value problem.

References