[Dune] ILU Preconditioner fails on zero diagonal entries

[Kilian Weishaupt]

Dear Christian and Simon,

thanks a lot for your fast response!

Concerning the details of my implementation:
I use a staggered-grid (MAC) approach.
The global jacobian matrix is a multitype-blockmatrix, including four 
large sub-blocks (pressure-pressure, pressure-velocity,
velocity-pressure and velocity-velocity). For UMFPACK and SuperLU, if 
have to copy everything into one single BCRSMatrix.
I currently do the same for GMRES. I'm not sure if entity wise blocking 
would work here, I need to have a look at that.
If it is of any use, I can also attach the matrix file so you might get 
an idea of what it looks like.

I haven't tried different preconditioners yet, thanks for the hint.

The idea of adapting the SuperLU backend in istl seems promising. I 
think I will give this a try first.

Thanks for your help!

Best regards
Kilian


On 04/20/2017 10:48 AM, Christian Engwer wrote:
> 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
>

-- 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Kilian Weishaupt M.Sc.
Institut für Wasser- und Umweltsystemmodellierung (IWS)
Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
Universität Stuttgart, Pfaffenwaldring 61, 70569 Stuttgart
Email: kilian.weishaupt at iws.uni-stuttgart.de
Telefon: 0049 711 685-60461 ** fax: 0049-711-685-60430
http://www.hydrosys.uni-stuttgart.de
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%