Algebraic Solution

The discretization scheme results in an algebraic system with $\left(\vphantom{ 4 n_c + n_b/2 }\right.$4n_c + n_b/2$\left.\vphantom{ 4 n_c + n_b/2 }\right)$ unknowns, where n_c is the number of cells in the problems and n_b is the number of boundary faces.² 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.