This paper presents the application of a preconditioned conjugate-gradient-like method to a non-self-adjoint problem of interest in underground flow simulation. The method furnishes a reliable iterative solution scheme for the non-symmetric matrices arising at each iteration of the non-linear time-stepping scheme. The method employs a generalized conjugate residual scheme with nested factorization as a preconditioner. Model runs demonstrate significant computational savings over direct sparse matrix solvers. KEY WORDS Generalized Conjugate Residual Method Nested Factorization Buckley-Leverett Equation Preconditioned Conjugate Gradients.
1. Introduction
Recent years have seen increasing interest in algorithms of the conjugate-gradient type for solving algebraic analogues of differential equations. This trend has been especially strong in petroleum reservoir simulation, where the need to solve large, sparse matrix equations arising from coupled sets of flow equations has sparked intensive research into iterative solution techniques. The use of the conjugate gradient method to solve linear systems Ju = -r dates at least to Hestenes and Stiefel,' who examine the standard algorithms applicable to symmetric, positive-definite matrices. Such matrices arise in the approximate solution of many elliptic partial differential equations, Laplace's equation being a prototype. We owe to Reid2 the view of the conjugate-gradient method as an iterative technique appropriate for sparse matrix systems and to Meijerink and van der Vorst3 the development of a practical salver using preconditioning to speed convergence in systems involving large, symmetric, positive-definite matrices.
However, in fluid flow problems the differential operators are rarely self-adjoint, and as a result discrete approximations typically give rise to non-symmetric matrices. The importance of such applications has motivated the development of a variety of techniques, related to the conjugategradient method, that accommodate non-symmetric matrices. Among these are Manteuffel's Chebyshev iteration m e t h ~d , ~ Kershaw's application of the conjugate-gradient procedure to the normal equations JTJu = -JTr,' Saad's incomplete orthogonalization methods,6 the biconjugate gradient (BCG) method presented by Fletcher' and a class of preconditioned iterative methods based on Elman's generalized conjugate residual (GCR) m e t h ~d . ~, ~
The GCR-based algorithms are attractive because none of them relies on computations involving JTJ, whose condition number may be much larger than that of the non-symmetric matrix J. Also, the GCR algorithm always converges provided the symmetric part of J is positive-definite, whereas with the BCG method