From: Jed B. <je...@59...> - 2008-04-29 11:18:20
|
This issue arose while discussing improved preconditioners for a colleague's libmesh-based 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 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. Perhaps some of the libmesh developers have ideas for the right way to abstract this method since it does require a lot of low-level PETSc-specific 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={1--137}, year={2005}, publisher={Cambridge Univ Press} } Jed |
From: John P. <jwp...@gm...> - 2008-04-29 14:10:49
|
Hi Jed, A few years ago I also wrote a (non-general) Schur complement preconditioner in LibMesh for the Navier-Stokes 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 Navier-Stokes problem? -J On Tue, Apr 29, 2008 at 6:16 AM, Jed Brown <je...@59...> wrote: > This issue arose while discussing improved preconditioners for a colleague's > libmesh-based 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 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. > > Perhaps some of the libmesh developers have ideas for the right way to abstract > this method since it does require a lot of low-level PETSc-specific > 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={1--137}, > 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 > _______________________________________________ > Libmesh-devel mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-devel > > |
From: Jed B. <je...@59...> - 2008-04-29 14:52:04
|
On Tue 2008-04-29 09:10, John Peterson wrote: > Hi Jed, > > A few years ago I also wrote a (non-general) Schur complement > preconditioner in LibMesh for the Navier-Stokes 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. Were you solving the Schur complement system completely? In my experiments, that is quite expensive and you might be better off just using a sparse approximate inverse (-pc_type spai or -pc_type hypre -pc_hypre_type parasails) instead. However, if the Schur complement is only solved approximately, I see large performance gains. > 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 Navier-Stokes problem? My current experiment was with a spectral collocation method because I was interested in preconditioning the high order nodal approximation and dealing with the nonlinear rheology in an unforgiving environment. I expect the results to carry over to p-refinement in an hp-element. My interest is primarily in low Reynolds number flow with non-Newtonian rheology. For high Reynold's number, the Schur complement needs a better preconditioner than the mass matrix. There have been lots of papers about this issue, for instance @article{elman2002paa, title={{Performance and analysis of saddle point preconditioners for the discrete steady-state Navier-Stokes equations}}, author={Elman, H.C. and Silvester, D.J. and Wathen, A.J.}, journal={Numerische Mathematik}, volume={90}, number={4}, pages={665--688}, year={2002}, publisher={Springer} } My experiments with the Stokes problem shows that it takes about four times as long to solve the indefinite Stokes system as it takes to solve a poisson problem with the same number of degrees of freedom. For instance, for a 3D problem with half a million degrees of freedom, the Stokes problem takes 2 minutes on my laptop while the poisson problem takes 30 seconds (both are using algebraic multigrid as the preconditioner). Note that these tests are for a Chebyshev spectral method where the (unformed because it is dense) system matrix is applied via DCT, but a low-order finite difference or finite element approximation on the collocation nodes is used to obtain a sparse matrix with equivalent spectral properties, to which AMG is applied. With a finite difference discretization ($PETSC_DIR/src/ksp/ksp/examples/tutorials/ex22.c) the same sized 3D poisson problem takes 13 seconds with AMG and 8 with geometric multigrid. This is not a surprise since the conditioning of the spectral system is much worse, O(p^4) versus O(n^2), since the collocation nodes are quadratically clustered. I expect that the rough factor of 4 will carry over to the low-order finite element world. My next experiment is to use a nodal Galerkin formulation and compare it to the Chebyshev case. I intend to incorporate it into an hp-element method where the high order elements are applied matrix-free. It is likely that I'll use libmesh for adaptivity, but I am using Sieve (from petsc-dev) for the data structures. I'll let you know how it goes. Jed |