Efficient Sparse-Matrix Implementation For Direct Methods In Reservoir Simulation
Authors: Graham Michael F. , Jennings Gray 
Published: December 2007

Abstract
In the past few years, some elimination orderings that represent improvements over the well-known grid-row ordering have appeared. Analyses using number of multiplicative operations as the theoretical measure of work demonstrate significant gains. However, implementation of these orderings for a computer requires programming techniques more sophisticated than those for a band solver.
In a recent paper, Alan George discusses experiments using dissection methods to solve square two-dimensional grid problems. In his experiments, the behavior of code-execution time reflected multiplication counts. But he did not determine - in practical sense - whether or not a problem-oriented implementation could be made superior to existing implementations for grid-row ordering. We have designed and produced an implementation of sparse-matrix techniques for direct solution of the pressure equation in a general-purpose three-phase reservoir simulator. Our implementation realizes in practice much of the gain promised in principle by the theoretical measure of work.
This paper discusses the design of our implementation and the impact of programming-language overhead. Performance data for nested-dissection and alternate-diagonal orderings in two and in three dimensions are given, and comparisons are made with existing band-algorithm and SIP implementations in the reservoir simulator.

Introduction
In the past few years, some elimination orderings that represent improvements over the well-known grid-row (or so-called standard) ordering have been published and analyzed. The authors' analyses, published and analyzed. The authors' analyses, using number of multiplicative operations as the theoretical measure of work, indeed demonstrate significant gains. A practical difficulty is that implementation of these orderings in a computer program requires programming techniques more program requires programming techniques more sophisticated than those for a band solver. If one is not careful in his design, gains expected from the improved orderings will be negated by overhead introduced in managing the ordering scheme.
In a recent paper, Alan George discusses experiments using dissection methods to solve two-dimensional grid problems. In his experiments, the behavior of his code-execution times reflected multiplication counts. But he did not determine - in practical sense - whether a problem-oriented implementation of dissection orderings could be made superior to existing ones for grid-row ordering. In this paper, we discuss an implementation of sparse-matrix techniques for direct solution of the pressure equation in a general-purpose three-phase reservoir simulator. Our implementation realizes in practice much of the gain promised in principle by the theoretical measure of work.

ORDERINGS OF INTEREST
The work expended in solving a linear system of equations is usually measured by counting the number of multiplicative operations required to effect the solution. A general expression for the work using Gaussian elimination is
(1)
where m is the number of equations in the system and ci is the number of nonzero coefficients in the i-th equation to the right of the main diagonal when the i-th unknown is eliminated. The ci(ci + 2) term enters from the so-called forward solution. There are ci equations from which the i-th unknown must be eliminated, and it takes ci + 2 operations to remove it from each of them. The ci + 1 term enters from the so-called backward solution.