[petsc-users] [Libmesh-users] Preconditioning of a Stokes problem

Mon Feb 4 12:30:39 CST 2013

Francesco,

Are you using any pressure stabilization scheme? If not, the 11-block in
your Jacobian
A = [A00 A01; A10 A11] would typically be zero, and preconditioning it with
jacobi wouldn't really work.

If A11 = 0, you ought to use the fieldsplit of type Schur, for example
(assuming petsc-3.3; prior versions use slightly different option names):
-pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type full
The question is then how you want to solve the Schur complement system.
Since the A11 block for you is zero, the easiest choice is to use S as the
preconditioning matrix:
-pc_fieldsplit_schur_precondition self

Preconditioning options for the Schur solver have prefix -fieldsplit_1_.
Since for Stokes the Schur complement is spectrally equivalent to the
pressure mass matrix, you can
hope that an unpreconditioned GMRES (default when using
-pc_fieldsplit_schur_precondition self) will solve it.

A better choice might be to use the Least Squares Commutator preconditioner
-fieldsplit_1_pc_type lsc
See
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCLSC.html
and note the relevant references at the bottom of that manual page.
You might also find this useful:
http://www.icfd.rdg.ac.uk/ICFD25/Talks/AWathen.pdf

PCLSC preconditions S = A10 inv(A00) A01 with  inv(A10A01) A10 A00 A01
inv(A10A01)
The "Laplacian" (a true Laplacian in the case of incompressible Stokes) can
be preconditioned
using AMG:
-fieldsplit_1_lsc_pc_type ml
(assuming you have ML built, e.g., with --download-ml)
or
-fieldsplit_1_lsc_pc_type gamg
with some appropriate choice of gamg options:
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCGAMG.html

To start with I would use -pc_fieldsplit_type schur
-pc_fieldsplit_schur_fact_type full -pc_fieldsplit_schur_precondition self
-fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml
gamg may also be very good in place of ml, but will require more option
setting.

You might find src/ksp/ksp/examples/tutorials/ex42.c and ex43.c useful,
albeit, possibly, a bit heavyweight for your needs.

A few more comments below:

On Mon, Feb 4, 2013 at 10:20 AM, John Peterson <jwpeterson at gmail.com> wrote:

> On Thu, Dec 27, 2012 at 4:31 AM, Francesco Ballarin
> <francesco.ballarin at mail.polimi.it> wrote:
> > Dear libMesh users,
> > I am trying to solve a Stokes problem in a 3D domain, but I am currently
> > struggling with the choice of the preconditioner.
> > I have slightly modified systems_of_equations_ex1 to accommodate for the
> > third velocity component and boundary conditions of the problem at hand
> > (Dirichlet BCs on the inlet on the walls, enforced using a
> > DirichletBoundary object, and homogeneous Neumann BC on the outflow)
> >
> > I am trying to use PETSc block preconditioners as shown, e.g., in
> > http://acts.nersc.gov/events/Workshop2012/slides/PETSc.pdf
> > section 4. I have tried most of the preconditioners discussed during the
> > talk, but I am not satisfied because I get either PETSc errors or a
> > performance that I hope can be increased.
> >
> > For example, with the following options:
> > *******
> >   PetscOptionsSetValue("-ksp_type", "gmres");
> > PetscOptionsSetValue("-ksp_right_pc", PETSC_NULL);
> >  PetscOptionsSetValue("-pc_type", "fieldsplit");
> > PetscOptionsSetValue("-pc_fieldsplit_detect_saddle_point", PETSC_NULL);
> >  PetscOptionsSetValue("-pc_fieldsplit_type", "multiplicative");
> > PetscOptionsSetValue("-fieldsplit_0_ksp_type", "preonly");
> >  PetscOptionsSetValue("-fieldsplit_0_pc_type", "ml");
> > PetscOptionsSetValue("-fieldsplit_1_ksp_type", "preonly");
> >  PetscOptionsSetValue("-fieldsplit_1_pc_type", "jacobi");
>

You do not need to set options to PETSC_NULL, if you do not intend to use
them.
Also, how are you forming the splits?
-pc_fieldsplit_detect_saddle_point is actually rather useful here.

Further, there are programmatic ways to set most of these options, for
example:
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFieldSplitSetType.html

Even so, we typically recommend setting these on the command line whenever
possible:
this way you do not hardwire your solver choices and can quickly experiment
with different
solvers from the command-line without recompiling your code.

>  > *******
> > I obtain a reasonable solution but after thousands of gmres iterations.
> Is
> > there any way to improve that?
> >
> > Or, if I try to embed the previous block preconditioner in mg:
> > *******
> > PetscOptionsSetValue("-pc_type", "mg");
> >  PetscOptionsSetValue("-pc_mg_levels", "5");
> > PetscOptionsSetValue("-pc_mg_galerkin", PETSC_NULL);
> >  PetscOptionsSetValue("-ksp_type", "gmres");
> > PetscOptionsSetValue("-mg_levels_pc_type", "fieldsplit");
> >  PetscOptionsSetValue("-mg_levels_pc_fieldsplit_detect_saddle_point",
> > PETSC_NULL);
> > PetscOptionsSetValue("-mg_levels_pc_fieldsplit_type", "multiplicative");
> >  PetscOptionsSetValue("-mg_levels_fieldsplit_0_ksp_type", "preonly");
> > PetscOptionsSetValue("-mg_levels_fieldsplit_0_pc_type", "gamg");
> >  PetscOptionsSetValue("-mg_levels_fieldsplit_1_ksp_type", "preonly");
> > PetscOptionsSetValue("-mg_levels_fieldsplit_1_pc_type", "jacobi");
> > *******
> > I get the following error:
>
Always send the full stack trace, when reporting errors.

Finally, in the future you might want to post these questions to
petsc-users at mcs.anl.gov

Dmitry.

> > *******
> > PETSC ERROR: Null argument, when expecting valid pointer!
> > PETSC ERROR: Null Object: Parameter # 2!
> > PETSC ERROR: MatPtAP() line 8237 in
> petsc-3.3-p4/src/mat/interface/matrix.c
> > *******
> >
> > Do you have any suggestion? What are you currently using as
> preconditioner
> > for saddle point problems?
>
> Apologies, it appears that this email fell through the cracks
> somewhere back around the holidays.
>
> Did you ever find a satisfactory answer to your question?
>
> Petsc developers may occasionally read this list, but for questions
> that are this advanced you would probably have better luck mailing the
> petsc developers list directly.
>
> And of course if you do find a good block preconditioning strategy for
> systems_of_equations_ex1, please feel free to report it here!
>
> --
> John
>
>
> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_jan
> _______________________________________________
> Libmesh-users mailing list
> Libmesh-users at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>

-- 
Dmitry Karpeev, Ph.D.
Assistant Computational Mathematician
Mathematics and Computer Science
Argonne National Laboratory
Argonne, Illinois, USA
and
Fellow
Computation Institute
University of Chicago
5735 S. Ellis Avenue
Chicago, IL 60637
-----------------------
Phone: 630-252-1229
 Fax:   630-252-5986
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130204/b985b326/attachment.html>