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