The discretization scheme results in an algebraic system with 4nc + nb/2 unknowns, where nc is the number of cells in the problems and nb is the number of boundary faces.2 This is roughly four times the number of unknowns that a method like the standard seven-point orthogonal operator produces, but the new method cannot be expressed in terms of only the cell center variables without obtaining a dense diffusion matrix. There are no methods which use only cell center unknowns that are second-order accurate on skewed meshes. The matrix system for this method is unsymmetric, which necessitates an unsymmetric solver. There is a maximum of 11 non-zero elements in any row, due to the cell face equation stencil.
When a Krylov subspace solver is used, a specialized preconditioning system has been developed to speed convergence. The preconditioner consists of a low-order version of the new method itself, which is derived by setting the minor direction terms for each face to zero (see section entitled ``Flux Terms''). This results in a matrix which can be reduced to a system involving the cell center unknowns only. In addition to having four times fewer unknowns, this system has a maximum of seven non-zero elements per row and is symmetric. After solution with a symmetric solver (which tends to be faster than an unsymmetric solver), the face unknowns are obtained via back-substitution.
The Krylov subspace solvers used for the solution of the main system are GMRES and TFQMR, among others. The low-order preconditioning system is solved using the Conjugate Gradients method, which is in turn preconditioned using SSOR. Alternately, an incomplete direct method called the unstructured multi-frontal method may be used. This solver gives very exact answers rapidly for small systems, but ultimately loses with the respect to the Krylov solvers when large time-dependent problems with loose tolerances are solved.