Multiscale linear solvers for very large systems derived from PDEs

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, multiscaling, 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 lower resolution approximation to the PDEs. While the coarse grained inverse is not known explicitly, the algorithm requires only 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. The local feedback provided by the conjugate gradient step at every level makes the algorithm very robust. In spite of the effort required for the coarse inverse, the algorithm is efficient because the increased quality of the approximate inverse greatly reduces the number of times the preconditioner needs to be evaluated. A feature of the algorithm is that the transition between coarse grids is determined dynamically by the accuracy requirement of the conjugate gradient solver at each level. Typically, later iterations on the finer scales need fewer iterations on the coarser scales and the computational effort is proportional to N rather than N log N, where N is the number of equations. We have tested our solver on the porous flow equation. On a workstation we have solved problems on grids ranging in dimension over three orders of magnitude, from 10³ to 10⁶, and found that the linear scaling holds. The algorithm works well, even when the permeability tensor has spatial variations exceeding a factor of 10⁹.