Efficient Use of Mass Storage During Elimination for Sparse Sets of Simulation Equations
Authors: McDonald A.E. , Trimble R.H. 
Published: December 2007

Abstract
Matrices too large to invert in central storage must be scratched onto mass storage. The scratched information must be ordered carefully to avoid overloading the computer's input/output system. Efficient user software is preferable to system software because matrix sparsity can be taken into account in a natural way. Storage management can be simplified by combining mass storage scratch with an appropriate sparse matrix factorization. This paper presents a factorization for the alternating-column ordering, and shows that it is faster than the alternating diagonal ordering (D-4) of Price and Coats when the grid dimensions are large. The method is applicable to both five- and nine-point operators, defined on a rectangular grid

Introduction
Recent publications have presented or evaluated methods for direct solution of linearized reservoir simulation equations. Practical aspects when systems are too large for central storage were not treated. For five-point operators, this paper presents an algorithm that is competitive with the presents an algorithm that is competitive with the alternate diagonal (D-4) method of Price and Coats, and that permits easy management of both central and mass storage. For both five- and nine-point operators, the method competes with nested dissection methods.
Matrices too large to invert in central storage must be scratched onto mass storage. The scratched information must be ordered carefully to avoid overloading the computer's input/output system. Efficient user software is preferable to system software since the special structure of the matrix can be taken into account. To some extent, a Virtual storage environment can reduce user concern over 1/0. Nevertheless, care must be taken to avoid excessive paging.
Woo et al. state that when their (Virtual) row algorithm is running alone, the elapsed time is about 1.5 times central processor (CPU) time. For their GNSO method, the factor increases to 2.0 times CPU time. Scratch penalties of 50 to 100 percent such as these are not essential when percent such as these are not essential when inverting large matrices. For the methods of this paper scratch penalties do not exceed 20 percent, paper scratch penalties do not exceed 20 percent, while central storage scratch requirements are very small.
For a five-point operator on a grid with dimensions I and J (where J less than I), central storage scratch requires 2J locations. lJ2/4 - J(J - 4)/12 mass storage locations are needed. These are requirements for the factorization and back substitution. Arrays that define the linear system may reside in central or mass storage. When in mass storage, their percentage contribution to scratch overhead is small for percentage contribution to scratch overhead is small for large I and J. Central processor requirements are made small by using a two-step procedure. In the first step, the "full grid" is subjected to odd-even reduction. This eliminates half the nodes, leaving a nine-point operator on a "reduced grid." Details of the reduction are given in Appendix A. The second step passes the reduced grid and its nine-point operator to the nine-point Zebra algorithm.
For a nine-point operator, the grid is given an ordering that permits efficient solution by elimination. The ordering alone does not prescribe a solution to the system of equations. A solution algorithm is also provided. It factors the matrix corresponding to the grid into the product of lower and upper triangular matrices, where the factors are sparse. It then uses sparse matrix arithmetic, in parallel with mass storage scratching, to effect an parallel with mass storage scratching, to effect an economical solution. George discussed a similar ordering, but gave no computational algorithm.

THE ZEBRA ORDERING FOR A NINE-POINT OPERATOR
The method is easy to visualize. In two dimensions, let the nodes be ordered by alternate columns, as indicated in Fig. 1. The columns with nodes denoted by () may be referred to as black stripes, and those with (O) as white stripes. This pattern of alternating black and white is called the zebra ordering.
SPEJ
P. 300