The increasing use of geostatistical techniques to infer physically meaningful realizations of porous media (e.g., petroleum reservoirs and subsurface aquifers) has served to highlight the weaknesses of early flow simulators. In particular, porous media typically exhibit fine-scale anisotropies that may not be aligned with the coordinate axes. Consequently, in the case of single phase saturated flow, it is well known from homogenization theory, and clearly demonstrated in practice, that the coarse-scale permeability will be, in general, a full tensor. In addition, many natural geological structures, such as fractures, significantly influence the flow field and should be preserved in the grid generation process. In many cases, a logically rectangular grid of irregular quadrilaterals (in two dimensions, irregular hexahedrals in three dimensions) is employed in the simulator. These two facts have placed new demands on the discretization methods that are used in single-phase saturated flow simulators.
In response to these demands a renewed interest in various discretization methods that preserve certain underlying mathematical properties of the model, such as mass conservation, has resulted. In fact, many of these sophisticated discretizations, (e.g., Mixed-Hybrid FEM, Local Support Operator and Nodal Methods) satisfy mass conservation locally over each cell in the grid. Typically, these discretizations treat the first order form (i.e., Darcy's Law and Mass Conservation) in an effort to obtain comparable accuracy in both the pressure and velocity. This approach leads to a saddle point problem for which the discrete system is indefinite. Localization (e.g., hybridization of Mixed FEMs), through the addition of scalar unknowns on the cell edges, facilitates the local elimination of the velocity without compromising the overall sparsity of the system. In fact, the resulting system for the pressure is symmetric positive definite, although, the sparsity structure is unusual, coupling unknowns on the edges and in the cells. Moreover, different stencils are associated with cells and edges of each orientation, suggesting that the discrete equations possess similarities with staggered grid discretizations of PDE systems. Thus, a difficult and important aspect of developing multigrid methods for these discretizations is the development of efficient relaxation algorithms. To this end we analyze the smoothing properties and computational cost of several relaxation schemes including, Gauss-Seidel, alternating lines, alternating block lines, distributed Gauss-Seidel and alternating distributed block lines. We limit the scope of our study to the lowest-order discretization of each family, noting that such a multigrid algorithm may be employed as a preconditioner for Conjugate Gradient iterations of higher order discretizations or other transport models.