We are starting to provide this functionality in petsc-dev and it will be in the next release. Barry On Apr 29, 2008, at 6:54 AM, Jed Brown wrote:

This issue arose while discussing improved preconditioners for a colleague's libmesh-based ice flow model. Since libmesh is one of several finite element libraries that use PETSc as the backend, perhaps the abstraction should be implemented in PETSc.

As a model, consider the Stokes problem with mixed discretization in matrix form

(1) [A, B1'; B2, 0] [u; p] = [f; 0]

where A is the viscous velocity system, B1' is a gradient and B2 is divergence. The system (1) is indefinite so most preconditioners and many direct methods will fail, however A is a uniformly elliptic operator so we can use exotic preconditioners. An effective preconditioner is to solve the Schur complement problem

(2) S p = -B2 A^{-1} f

where S = -B2 A^{-1} B1' is the Schur complement. Once p is determined, we can find u by solving

(3) A u = f - B1' p

If we solve (2,3) exactly, this is the ultimate preconditioner and we converge in one iteration. It is important to note that applying S in (2) requires solving a system with A. This can be implemented in PETSc using a MatShell object. The classical Uzawa algorithm is a Richardson iteration on this problem. A better method (in my experience) is to precondition a Krylov iteration (FGMRES) on (1) with an approximate solve of (2,3). For instance I have had success with 3 iterations of GMRES on (2) where the A^{-1} in S is replaced by one V-cycle of algebraic multigrid and the A^{-1} appearing explicitly in (2,3) are replaced by 4 iterations of GMRES preconditioned with AMG.

It is desirable to have all the details of this iteration controllable from the command line. For an example implementation (a spectral method) of this problem, check out StokesMatMultXXX() and StokesPCApply() at

http://github.com/jedbrown/spectral-petsc/tree/master/stokes.C

Note that typically, S is preconditioned with the mass matrix (or a ``pressure-laplacian'' at higher Reynolds number) but since this code uses a collocation method and is only intended for slow flow, there is no preconditioner for S.

My question for the PETSc developers: how much of this can/should be abstracted out? While it's not difficult to code, perhaps there should at least be an example. If it would be useful, I can strip my experiment above down to make it suitable.

PS: An excellent review paper: @article{benzi2005nss, title={{Numerical solution of saddle point problems}}, author={Benzi, M. and Golub, G.H. and Liesen, J.}, journal={Acta Numerica}, volume={14}, pages={1--137}, year={2005}, publisher={Cambridge Univ Press} }

Jed