Preconditioning for saddle point problems

[Barry Smith]

   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