[petsc-users] Many issues solving matrix derived from a CFD coupled Solver

[Jed Brown]

Some aspects of your questions come up frequenty, especially the
unrealistic expectations from black-box iterative solvers.  Hopefully
this somewhat in-depth reply will be useful.

On Thu, 27 May 2010 09:20:55 +0200, Luca Mangani <luca.mangani at hslu.ch> wrote:
> Good Morning,
> 
> After many steps I'm not able to solve my coupled matrix with
> iterative solver. The matrix is defined as the discretization of the
> Navier Stokes equation in a pressure based form.

This form is indefinite.  While we often talk about there being no
scalable black-box solvers, the statement is especially true for
indefinite systems.

As an example of the spectacular failure of black-box methods for
elliptic systems when applied to other problems, I have observed
BoomerAMG being massively singular when applied to a mixed finite
element discretization of incompressible flow.  This resulted in
preconditioned residuals falling in essentially textbook manner, but the
unpreconditioned residuals did not converge at all (after satisfying the
boundary conditions).  It is very important to check this.

For *some* incompressible flow problems, it is possible to order
unknowns such that incomplete factorization causes useful fill into the
pressure-pressure block.  Unfortunately, incomplete factorization has
very unpredictable properties and is certainly not a robust solution.
When playing with orderings, usually you would like all velocity
unknowns in a subdomain to come before the pressure unknowns, a number
of automatic orderings are also available with
-pc_factor_mat_ordering_type.  You may also try changing the shift
behavior, see -pc_factor_shift_type, using an ILU with pivoting, or
using sparse approximate inverse -pc_hypre_type parasails or -pc_type
spai.

While some of these black-box approaches will be better (worse) than
others, and depending on your needs, perhaps even "good enough", none
will scale to large problem size or be robust to parameters.  They may
all fail completely.  There are essentially two scalable approaches:

1. Schur complement methods

  This class contains classical methods like Uzawa and SIMPLE, as well
  as newer approaches based on approximate commutators.  See Elman's
  2008 paper "A taxonomy and comparison of parallel block multi-level
  preconditioners for the incompressible Navier-Stokes equations" for a
  current overview.  These methods break the problem into definite
  operators for which conventional preconditioners are suitable.  It is
  relatively easy to implement using PCFieldSplit, but note that the
  approximation to the Schur complement is a crucial choice and involves
  problem-specific information.

2. Coupled multilevel methods

  This class contains multigrid methods based on Vanka smoothers as well
  as domain decomposition methods.  It is crucial to define subdomains,
  interpolants, and smoothers/subdomain solvers that are compatible with
  the inf-sup condition.  These stability requirements are highly
  dependent on your spatial discretization, therefore it really cannot
  be encapsulated in a black-box solver.  For the multigrid end of the
  spectrum, see papers by Brandt in the finite difference context or
  Rannacher and Turek for certain finite element discretizations.  At
  the DD end, you might like papers by Jing Li for the non-overlapping
  case or Xiao-Chuan Cai for overlapping cases.

  Note that 1-level Schwarz methods essentially only require a choice of
  stable subdomains which you can do with PCASMSetLocalSubdomains().
  Many other methods can be posed as a multilevel Schwarz method in
  which case you can use PCMG and define your own problem-specific
  interpolants.

There are currently no scalable algorithms with convergence rates that
are independent of Reynolds number, the best cases with the latter class
(coupled multigrid or DD) seem to be somewhat less sensitive than the
Schur-complement methods.  As usual, there are difficulties in the
presence of large variation in coefficients or high aspect ratio, though
some methods are indpendent or very weakly dependent on these
parameters.  The methods in class 1 almost always involve more frequent
communication, therefore their scalability is likely to be somewhat
worse than class 2, though some people have demonstrated quite good weak
scalability using class 1.

Benzi, Golub, and Liesen's 2005 Acta Numerica review paper may also be
of interest.


More things are likely to work with non-mixed discretizations, but those
discretizations are generally less robust.  A group at Sandia has had
some luck applying ML to stabilized Q1-Q1 finite element methods, it is
important to keep all fields coupled (PETSc will inform ML correctly if
you order unknowns [u0,v0,w0,p0,u1,...] and call MatSetBlockSize or use
BAIJ, note that this ordering is better for memory performance and
smoother strength anyway).  Their experience is that this approach is
usually faster than approximate commutators, but sometimes fails to
converge at very moderate parameter choices.  Note that this Q1-Q1
discretization does not exactly enforce incompressibility and this
divergence error can be quite large on under-resolved meshes or in the
presence of non-smooth boundary conditions.

Galerkin least squares methods offer a way to produce a positive
definite system from a velocity-pressure formulation, therefore standard
solvers can be applied.  But note that it generally commits very large
divergence errors (see 99% mass loss a few years ago, the newer methods
can keep it to a few percent for simple problems), this can be
circumvented by choosing inf-sup stable spaces, but then the solver
issues become somewhat more complicated again.

Once you have decided on a discretization and algorithm from the
literature that you would like to implement (or your own variant), we
can give suggestions.  Please recognize that scalable solvers for
incompressible flow is very much an active research topic, you cannot
expect black-box solvers to produce good results.  If you do not need
very large problem sizes, you might want to stick with a direct solver
like MUMPS; if do not need extreme parameter choices, you might get by
with naive 1-level domain decomposition or coupled algebraic multigrid
(especially with collocated discretization).

Jed