From: John Peterson <jwpeterson@gm...>  20080429 14:10:49

Hi Jed, A few years ago I also wrote a (nongeneral) Schur complement preconditioner in LibMesh for the NavierStokes problem, based around ex13. I built the normal global stiffness matrix and made use of LibMesh's PetscMatrix::create_submatrix() and PetscVector::create_subvector() routines to logically work with the individual A, B, and B^t blocks, as well as a MatShell object you mentioned below. This code will most likely no longer compile as is, but you would be welcome to take a look if you are interested. It also doesn't appear to have been working (successfully solving problems) at the time I stopped using it ... I'm not sure if this was from a bug in the implementation, or if not having a good preconditioner for the "S" system was killing me at the time. I think our best shot at generalizing would be within the FEMSystem framework. If I were going to do it myself, I would probably use a trick similar to what's been done in the EigenTimeSolver class, which also assembles multiple matrices (using a flag to control the assembly) for solving generalized eigenproblems. Do you have any results for Schur complement preconditioning in FEM with a particular combination of preconditioners and solvers that significantly accelerates the solution of the Stokes or NavierStokes problem? J On Tue, Apr 29, 2008 at 6:16 AM, Jed Brown <jed@...> wrote: > 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 > >  > This SF.net email is sponsored by the 2008 JavaOne(SM) Conference > Don't miss this year's exciting event. There's still time to save $100. > Use priority code J8TL2D2. > http://ad.doubleclick.net/clk;198757673;13503038;p?http://java.sun.com/javaone > _______________________________________________ > Libmeshdevel mailing list > Libmeshdevel@... > https://lists.sourceforge.net/lists/listinfo/libmeshdevel > > 