Introduction and need
Almost any method used for the discretization of the Navier-Stokes equations (and most other physical phenomena) yield a linear system of equations that needs to be solved. The need for robust, efficient, and relatively fast linear solvers becomes a crucial issue. While some linear systems exhibit specific behavior (e.g. symmetric, positive definite, tridiagonal...), the general rule excludes this specialty. For example, matrices arising from the finite volume method when ussed on unstructured grids are unsymmetric, sparse, and sometimes singular. Therefore, a huge body of research dedicated for designing robust algorithms has been accumulated over the last two decades.
Without loss of generality, most linear systems arising form the FV method are sparse (8 or 9 nonzero elements per row at most for a 3D solver). It is thus essential to plan ahead and choose a general solver that is able to handle the sparsity and skewness of the coefficient matrix.
The choice of implementors is an iterative solver. Direct solvers are computationally expensive, and inefficient for large sparse matrices. The most compelling reason for choosing iterative solvers is the inherent non-linearity of the coefficient matrix. As the coefficient matrix is updated at the end of each outer iteration, it is unjustifiable to spend the extra cost for a direct solution that will be iteratively driven to the final solution.