From: Jed Brown <jed@59...>  20080429 11:18:20

This issue arose while discussing improved preconditioners for a colleague's libmeshbased ice flow model. It seems to me that libmesh ought to have an abstraction for this since normal preconditioners are not effective. 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, B2 is the weak 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 Vcycle 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/spectralpetsc/tree/master/stokes.C Note that typically, S is preconditioned with the mass matrix (or a ``pressurelaplacian'' 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. Perhaps some of the libmesh developers have ideas for the right way to abstract this method since it does require a lot of lowlevel PETScspecific manipulation that ought to be hidden from the typical libmesh user. 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={1137}, year={2005}, publisher={Cambridge Univ Press} } Jed 