The Multigrid Solvers
The different Multigrid solvers types — geometric multigrid (GMG) and algebraic multigrid (AMG) solvers— are discussed in this section as well as the multigrid algorithm.
The Geometric Multigrid Solver/Preconditioner
The geometric multigrid solver uses a hierarchy of multigrid levels where each level corresponds to a mesh and a choice of shape functions. Thus, in addition to coarsening the mesh, it is possible to construct a new “coarser” level by lowering the order of the shape functions. The number of degrees of freedom decreases when you go to a coarser multigrid level. The meshes for the different levels can be constructed manually or automatically. The automatic options use a coarsening algorithm to the fine mesh, which leads to meshes that are not aligned to each other. There is also an option to generate the finer meshes from the coarsest mesh by successive mesh refinements, which leads to aligned (nested) meshes. The manual option can be useful when you have a quadrilateral, hexahedral, or prism mesh, or when you for some other reason want to control details in the meshes.
The geometric multigrid (GMG) solver or preconditioner is a fast and memory-efficient iterative method for elliptic and parabolic models. It performs one or several cycles of the geometric multigrid method. The classical multigrid algorithm uses one or several auxiliary meshes that are coarser than the original (fine) mesh. The idea is to perform just a fraction of the computations on the fine mesh. Instead, it performs computations on the coarser meshes when possible, which leads to fewer operations. The size of the extra memory used for the coarser meshes and associated matrices is comparable to the size of the original data. This leads to an iterative algorithm that is both fast and memory efficient. See Ref. 18 for more information.
The Algebraic Multigrid Solvers/Preconditioners
The algebraic multigrid (AMG) solvers or preconditioners perform one or several cycles of the algebraic multigrid method. This is similar to the geometric multigrid algorithm, the difference being that it constructs the multigrid levels directly from the finest-level system matrix A0. That is, it constructs the prolongations Pi from A0 without using auxiliary meshes. It constructs the coarse level matrices Ai from A0 with the Galerkin projection method. The advantage is that you need not bother about the coarse multigrid levels.
The Multigrid Algorithm
To describe the multigrid algorithm, assume that you have N + 1 multigrid levels numbered from 0 to N, where 0 is the finest level (the level for which you seek the solution). To solve the linear system A0x = b (corresponding to level 0), the algorithm must reform the system matrices A1, …, AN for the coarse multigrid levels. It must also compute the prolongation matrices Pi that map a solution x vector on level  i to the corresponding solution vector Pi x on the next finer level i − 1.
The prolongation matrices are constructed using plain interpolation from one multigrid level to the other. The system matrices for the coarse levels can be constructed in two ways:
•
By assembling Ai on the mesh of level i (the default method).
•
By projection from the finer level: Ai = PiTAi−1Pi. This is also called the Galerkin method. It typically leads to more nonzero elements in the system matrix Ai, but the convergence can be faster than in the default method.
The following algorithm describes one multigrid cycle:
1
The input to the algorithm is some initial guess x0 to the solution of the system A0x = b.
2
Starting with x0, apply a few iterations of a presmoother to the linear system A0x = b, yielding a more accurate iterate x0s. Typically the presmoother is some simple iterative algorithm such as SOR, but you can also choose any iterative solver.
3
Compute the residual r0 = b − A0 x0s. The presmoother “smooths” the residual so the oscillations in r have such a long wavelength that they are well resolved on the next coarser level (1). Therefore, project the residual onto level 1 by applying the transpose of the prolongation: r1 = P1Tr0.
4
If N = 1 use the coarse solver to solve the system A1x1 = r1. The coarse solver is typically a direct solver such as MUMPS. The number of degrees of freedom on level 1 is less than for level 0, which means that solving A1x1 = r1 is less expensive. If instead N > 1, solve the system A1x1 = r1 (approximately) by recursively applying one cycle of the multigrid algorithm for levels 1, 2, …, N. In both cases the obtained solution x1 is called the coarse grid correction.
5
Map the coarse grid correction to level 0 using the prolongation matrix: x0c = x0s + P1x1.
6
Starting with x0c, apply a few iterations of a postsmoother to the linear system A0x = b, yielding a more accurate iterate x0mg. The default postsmoother is SORU (the version of SOR using the upper triangle of the matrix). The iterate x0mg is the output of the multigrid cycle.
The cycle just described is called the V-cycle. The recursive call in step 4 (when N > 1) is a also a V-cycle. For the W-cycle and the F-cycle, steps 1–6 above are the same but with the twist that the recursive call in step 4 is substituted with two multigrid calls for the coarser levels. For the W-cycle these two calls are recursive calls (W-cycle calls). For the F-cycle the first call is a W-cycle and the second a V-cycle.
For only two multigrid levels (N = 1) these cycles are the same because the algorithm uses the coarse solver in step 4. Also the amount of work on the finest level is the same for the different cycles. Normally the V-cycle is sufficient, but the W-cycle and the F-cycle can be useful for more difficult problems.
When using multigrid as a preconditioner, the action of this preconditioner is obtained by applying a fixed number of multigrid cycles. When using multigrid as a solver, the multigrid cycle repeats until it reaches convergence.
When using multigrid as a preconditioner for the conjugate gradients method for a symmetric matrix A, the preconditioning matrix M should also be symmetric. This requirement is fulfilled if the matrices M associated with the presmoother and the postsmoother are transposes of each other. For instance, this is the case if the presmoother is SOR and the postsmoother is SORU, and if the same number of smoothing steps is used. This combination with two smoothing steps is the default.
Notes on the Efficiency of Smoothers
COMSOL Multiphysics performs smoothing on all but the coarsest multigrid level. A smoother should be computationally cheap and effective at reducing the part of the error that has a high spatial frequency on the mesh to which it is applied. Therefore, applying a smoother on several meshes with a hierarchy of mesh sizes results in a more efficient solver than if the smoother were applied only on the finest mesh.
The efficiency of the multigrid method with simple iterations as a smoother (such as the Jacobi and SOR iteration) hinges on the ellipticity of the underlying mathematical problem. For Helmholtz problems originating from an equation
or
the obtained linear problem is indefinite for large frequencies ω. For these problems, a simple iteration amplifies smooth eigenmodes if the mesh is too coarse and makes these methods unsuitable as smoothers. To determine when to use a simple iteration, apply the Nyquist criterion:
which says that the mesh must have at least two mesh elements per wavelength. Thus, when using the geometric multigrid solver for these types of problems, ensure that this criterion is fulfilled on the coarsest mesh if simple iterations are used as a smoother. In situations where the criterion is not fulfilled on coarse meshes, GMRES can be a suitable smoother (Ref. 21). However, this setting makes smoothers on all levels more expensive and might not always pay off compared to choosing a coarse grid that satisfies the Nyquist criterion. Note also that a smoother based on a Krylov preconditioner like GMRES requires the (outer) iterative solver to be set to FGMRES.
The Complex Shifted Laplace Operator
The convergence of the Multigrid method can be improved by adding a complex shifted Laplace (CSL) contribution — that is, a purely complex term that has the effect of damping oscillations in the solution. CSL is typically applied to wave problems at high frequency. Its effect is limited to the preconditioner and the convergence rate, and it can be applied independently both on the fine and on the multigrid levels. CSL does not change the original system matrix and the final solution. The CSL contribution is a function of the wave number of the problem and is generally multiplied by a relaxation factor O(1). The choice of the relaxation factor is a trade-off between no damping (0) and large damping but deterioration of the preconditioner performance (large values).