Multi-scale linear solvers for very large systems derived from PDEs

TL;DR

The paper addresses solving very large linear systems $A x = b$ arising from PDE discretizations by introducing a recursive, multi-scale preconditioned conjugate gradient method. It constructs a high-quality preconditioner $M^{-1}$ from a coarse-grid inverse and smoothing, implemented recursively so that coarse problems are solved with the same algorithm. The approach achieves near-linear scaling and robust convergence on porous-flow tests, handling large coefficient contrasts and anisotropy; it also outperforms standard multigrid in several scenarios and remains applicable to broader PDEs. This work provides a scalable, flexible framework that blends multiscale decomposition with Krylov methods and is amenable to extensions such as parallelization and nonsymmetric systems.

Abstract

We present a novel linear solver that works well for large systems obtained from discretizing PDEs. It is robust and, for the examples we studied, the computational effort scales linearly with the number of equations. The algorithm is based on a wavelength decomposition that combines conjugate gradient, multi-scaling and iterative splitting methods into a single approach. On the surface, the algorithm is a simple preconditioned conjugate gradient with all the sophistication of the algorithm in the choice of the preconditioning matrix. The preconditioner is a very good approximate inverse of the linear operator. It is constructed from the inverse of the coarse grained linear operator and from smoothing operators that are based on an operator splitting on the fine grid. The coarse graining captures the long wavelength behavior of the inverse operator while the smoothing operator captures the short wavelength behavior. The conjugate gradient iteration accounts for the coupling between long and short wavelengths. The coarse grained operator corresponds to a lowerresolution approximation to the PDEs. While the coarse grained inverse is not known explicitly, the algorithm only requires that the preconditioner can be a applied to a vector. The coarse inverse applied to a vector can be obtained as the solution of another preconditioned conjugate gradient solver that applies the same algorithm to the smaller problem. Thus, the method is naturally recursive. The recursion ends when the matrix is sufficiently small for a solution to be obtained efficiently with a standard solver. We have tested our solver on the porous flow equation. On a workstation we have solved problems on grids ranging in dimension from $10^3$ to $10^6$, and found that the linear scaling holds.

Figures (8)

• Figure 4.1: Flow lines (in black) superimposed on the log of permeability field for the base case. The permeability field is random log-normal with zero mean, a variance of 2 and a power law auto-correlation. The field is discretized on a $1000\times 1000$ grid and its log ranges from a low of -7.3 (blue) to a high of 7.9 (red). We note for this "fractal" permeability field the flow lines tend to follow in narrow channels. This is in contrast to our experience with Gaussian auto-correlations. (The raggedness of the flow lines is a consequence of the compression algorithm used to reduce the size of the plot file and is not an indication of either the resolution or the accuracy of the solution.)
• Figure 4.2: Level tree showing the dynamic transition between levels for the base case. Level 0 is the finest grid and level 4 is the coarsest grid.
• Figure 4.3: The effect of the scale factor between grids on the base case. Shown are the number of levels, iterations on the finest grid, and the time per point.
• Figure 4.4: Number of iterations on the finest grid and time per point versus number of grid points for the recursive multi-scale conjugate gradient algorithm.
• Figure 4.5: Number of iterations on the finest grid and time per point versus number of grid points for the multi-grid preconditioned conjugate gradient (Tatebe's) algorithm.