Re: [Libmesh-devel] Preconditioning of saddle point problems

[Jed B. <je...@59...>]

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