Preconditioning for saddle point problems

[Jed Brown]

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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 197 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20080429/c75d7f29/attachment.sig>