[Dune] ILU Preconditioner fails on zero diagonal entries

[Christian Engwer]

Dear Kilian,

a) improved preconditioners are always welcome

b) perhaps you don't even need this...

> I have implemented a fully implicit incompressible Stokes solver in Dumux,
> yielding a global matrix with some zero entries on its
> diagonal. Currently, I use SuperLU or UMFPACK to solve it, which gets of
> course very slow for larger systems.

what kind of discretization do you use exactly. If you can do entity
wise blocking and store everything in a
BCRSMatrix<FieldMatrix<double,N,N>> you should be finde. This works
for example for stabilized Q1/Q1 with N=4 (3 velocity
+ 1 pressure component) or for a stable DG-P2/DG-P1 pair. The block
will be inverted directly, which works, and the ILU-0 runs trough
smoothly.

> Apparently, many people use an ILU-preconditioned GMRES solver for this kind
> of application. However, the dune ILU preconditioner
> runs into problems, since it divides by zero (the diagonal entries) at some
> points.
> 
> I have then tried to store my matrix in a file and solve it using python's
> ILU-preconditioned GMRES which works rather well.
> Their ILU implementation seems to rely on SuperLU ( https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.spilu.html
> )
> 
> Furthermore, there is a paper addressing the issue of zero pivots (
> http://crd-legacy.lbl.gov/~xiaoye/slu_ilu.pdf , "3.4 Handling breakdown due
> to zero pivots" ). I have tried to implement their solution strategy in
> ilu.hh. However, I am no expert on that and this only seems to work for
> small systems in my case.

I think the approach sketched by Simon should be fairly easy to
implement. The downside is, tat you have to store the matrix twice
(once again in the SuperLU format) and that you loose any information
about the block structure of your problem.

An other apporach (which is the classical one) is to use a
Schur-Complement solver. You can look for example at the Rannacher
Numerics 3 lecture notes.

Ciao
Christian