Stabilization of Agglomerative Algebraic Multigrid Solver by Recursive Projection Method

Aleksandar Jemcov

Fluent Inc., 10 Cavendish Court, Lebanon, NH, 03766, USA

Joseph P. Maruszewski
Hrvoje Jasak


Abstract

The Classical approach to constructing an Algebraic Multigrid Solver (AMG) for the solution of linear system of equations is to select a smoother and devise an agglomeration strategy to follow the smooth error component on successive levels. Performance of an AMG solver is dependent on the ability of agglomeration to closely represent the distribution of the smooth error components on a given level and, to a lesser degree, on a choice of smoother. Successful agglomeration strategies result in better performance of AMG solvers.

Agglomeration strategies are hardly universal and depend on the physical properties of the problem and numerical properties of the discretization technique used to create the system of linear equations under consideration: agglomeration that leads to success in one case may lead to failure in another. Furthermore, complex agglomeration may incur a computational penalty associated with the construction of restriction and prolongation operators, with potentially unacceptable computational time overhead. This is particularly important in time-dependent solution of partial differential equations, where a linear system is solved many times per time step. Additional problems are associated with the need to tune the agglomeration, when the primary goal is the analysis of the underlying problem and not a search for new agglomeration strategies. It is therefore desirable to to devise a computationally efficient agglomeration robust enough to handle wide range of problems.

One simple and popular matrix coarsening algorithm is the Agglomerative Method. the basic assumption in agglomerative multigrid is the homogeneity of the smooth error distribution within the agglomeration stencil. The result of the application of this agglomeration strategy is that every fine level equation will have a representation on a coarse level and moreover, coarse level corrections will be constant within a given interpolation stencil. Homogeneous smooth error distribution assumption is correct only in a very special cases of strongly elliptic problems with constant coefficients discretized on uniform grids. In practical application this assumption is rarely met and an Algebraic Multigrid solver that uses this agglomeration strategy can easily fail to converge. The attraction of this agglomeration strategy is in its simplicity and computational efficiency since prolongation and restriction operators are known a priori thus reducing the overall number of operations per multigrid cycle resulting in a very fast algorithm. This agglomeration strategy is effective if the discretization method produces linear system of equation whose matrix has M-matrix properties. In that case, simple scaling of corrections based on variational considerations can substantially improve performance of the multigrid solver. However, if the underlying matrix violates M-matrix properties, the assumption of homogeneous smooth error distribution within a stencil is not applicable and agglomeration strategy that closely follows smooth error distribution is more effective.This, of course, comes at increased computational time and larger memory usage.

One solution to this problem is to use Agglomerative Algebraic Multigrid Solver as a preconditioner to Krylov subspace methods. This approach allows for flexibility of the selection of the most appropriate Krylov subspace solver that fits the needs of the problem at hand and significant performance gain results from the proper choice of the solver and preconditioner. However, often the spectrum of the iteration matrix may be influenced by the quality of the computational grid and peculiarities of the physical models used in the discretization process thus making the process of selection of the appropriate Krylov subspace method harder. When this is the case, a careful choice of the type of the cycle, smoother and agglomeration stencil for the multigrid preconditioner together with appropriate Krylov subspace solver will strike the right balance between computational cost and efficiency. In principle, these choices are not known a priori and every problem requires some amount of tuning of the solver and preconditioner parameters.

Another solution to this problem is in the enhancement of the Agglomerative Algebraic Multigrid Solver by Recursive Projection Method (RPM). The idea of the RPM consists of the projection of the solution vector and fixed-point function onto stable and unstable spaces thus separating unstable eigenvalues of the iteration matrix from the stable ones. The principal question in this approach is the one of the construction of stable and unstable orthogonal projectors. Projectors are constructed from the recursive property of the the algebraic multigrid algorithm which states that the error propagation through iterations is according to powers of the iteration matrix. If the error is defined as the difference between two consecutive solves, then these differences form the Krylov basis of the problem. The dominant eigenspace is extracted from this basis by performing the QR-decomposition of the rectangular matrix whose columns are difference vectors. Diagonal entries of the upper triangular matrix R in QR-decomposition are used to determine which columns of the matrix Q span the unstable space. This procedure is performed when the AMG starts experiencing stalling or divergence of L2 norm of the linear system residual. Once the unstable basis is determined, the unstable projector is constructed and stable projector is obtained from the fact that these projectors are orthogonal.With the knowledge of stable and unstable projectors, the fixed-point function of Agglomerative Algebraic Multigrid Solver together with the current solution is projected on stable and unstable spaces. The stable part of the algorithm proceeds to iterate according to AMG algorithm with stable part of the solution, whereas for the unstable part of the solution the special Newton algorithm is constructed. At the convergence, full solution is constructed from unstable and stable parts by the direct summation.

RPM solver is implemented as a wrapper around agglomerative AMG solver in such way that it becomes active only when divergence or stalled convergence is detected. In that case, dominant eigenspace is determined and stable and unstable projectors are constructed and solution is projected onto those spaces. The basis of unstable space is recursively enriched until all unstable modes are captured at which time rapid convergence is experienced. The advantage of this method is in its dynamic nature since it is activated only when needed. In contrast to preconditioned Krylov space solvers, RPM is universal in nature and very efficient from the computational point of view.