[Dune] [#928] breakdown of BiCGSTAB

[Christian Engwer]

Hi Bernd,

did you try reducing the breakdown criterion? In other implementations
it is '== 0' instead of '< epsilon'. Does it work then, read "it is
just the choice of the epsilon", or does it still fail, read "there
is a real bug in the BiCGStab"?

Christian

On Thu, Jun 09, 2011 at 02:14:05PM +0200, Bernd Flemisch wrote:
> Hey Markus,
> 
> thank you for your answer. Although you apparently think that my
> answer is not interesting since you already closed this task, I
> provide it here anyway.
> 
> This matrix and right hand side are coming from a to my knowledge
> correct CVFE discretization of a simple diffusion equation on the
> unit square discretized with 8x4 rectangular elements. Dirichlet BCs
> on top and bottom, Neumann zero left and right. The solution of the
> equation system is the physically correct one, meaning constant in x
> and linear in y direction.
> 
> You are somehow right about suggesting a different right hand side.
> The unusual form is coming from the fact that our Dirichlet BCs are
> implemented a bit differently leading to 1's and 2's on the
> diagonals in the original example. I changed that in the newly
> attached example with only 1's on the Dirichlet diagonals and your
> suggested rhs. BiCGSTAB still fails.
> 
> Matlab's bicgstab works.
> 
> This is a boiled down simple example, and of course I could use CG
> here. I might not for more complex problems where I have to rely on
> BiCGSTAB or GMRes.
> 
> Thanks for maybe taking one more look. Kind regards
> Bernd
> 
> 
> On Thu, 09 Jun 2011 11:39:54 +0200
>  Dune <flyspray at dune-project.org> wrote:
> >THIS IS AN AUTOMATED MESSAGE, DO NOT REPLY.
> >
> >The following task has a new comment added:
> >
> >FS#928 - breakdown of BiCGSTAB
> >User who did this - Markus Blatt (mblatt)
> >
> >----------
> >Does your right hand side make sense from the physics of the
> >problem.
> >
> >Whenever this happened to us, the discretisation was wrong. Most
> >often the left and right hand side are not compatible in terms of
> >boundary conditions if breakdowns occur.
> >
> >try e.g. X=0; x[0:8] = -1e+05; which makes more sense and
> >satisfies BiCGSTAB!
> >
> >With BiCGStab breakdowns may occur. But this is rather a problem
> >of the algorithm than our implementation.
> >If CG works, you should therefore use it.
> >
> >Did you also check bicgstab in matlab?
> >----------
> >
> >More information can be found at the following URL:
> >http://www.dune-project.org/flyspray/index.php?do=details&task_id=928#comment2595
> >
> >You are receiving this message because you have requested it from
> >the Flyspray bugtracking system.  If you did not expect this
> >message or don't want to receive mails in future, you can change
> >your notification settings at the URL shown above.
> >
> >_______________________________________________
> >Dune mailing list
> >Dune at dune-project.org
> >http://lists.dune-project.org/mailman/listinfo/dune
> 
> ___________________________________________________________
> 
> Bernd Flemisch                     phone: +49 711 685 69162
> IWS, Universitaet Stuttgart          fax: +49 711 685 67020
> Pfaffenwaldring 61        email: bernd at iws.uni-stuttgart.de
> D-70569 Stuttgart        url: www.hydrosys.uni-stuttgart.de
> ___________________________________________________________

> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
> 
> #include "config.h"
> #include <dune/istl/bcrsmatrix.hh>
> #include <dune/istl/io.hh>
> #include <dune/istl/operators.hh>
> #include <dune/istl/solvers.hh>
> 
> /*!
>  * \brief fill a BCRS matrix from a Matlab input file.
>  *
>  * An input file of the form
>  *
>  * i0 j0 a_i0j0
>  * ...
>  * in jn a_injn
>  *
>  * is read, where
>  * i0, ..., in are row indices,
>  * j0, ..., jn are column indices,
>  * a_i0j0, ..., a_injn are corresponding entries,
>  * and stored into a BCRS matrix.
>  * The indices have to start with 1.
>  *
>  * \param fileName name of the file containing the matrix indices and entries
>  * \param A the matrix to be filled
> */
> template<class Matrix>
> void readBCRSMatrixFromMatlabFile(const char *fileName, Matrix& A)
> {
>     typedef typename Matrix::size_type size_type;
>     typedef typename Matrix::field_type Scalar;
> 
>     std::vector<size_type> row;
>     std::vector<size_type> col;
>     std::vector<Scalar> entry;
> 
>     std::ifstream input(fileName);
>     if (!input)
>         DUNE_THROW(Dune::Exception, "could not open file " << fileName);
> 
>     while(input)
>     {
>         size_type i;
>         if (!(input >> i))
>             break;
>         row.push_back(i-1);
> 
>         size_type j;
>         if (!(input >> j))
>             break;
>         col.push_back(j-1);
> 
>         Scalar aij;
>         if (!(input >> aij))
>             break;
>         entry.push_back(aij);
>     }
> 
>     size_type n = 1 + *(std::max_element(row.begin(), row.end()));
>     size_type nnz = row.size();
> 
>     A.setSize(n, n, nnz);
>     A.setBuildMode(Matrix::random);
> 
>     for (size_type i = 0; i < n; i++)
>     {
>         size_type cols = count(row.begin( ), row.end( ), i);
>         A.setrowsize(i, cols);
>     }
>     A.endrowsizes();
> 
>     for (size_type i = 0; i < nnz; i++)
>         A.addindex(row[i], col[i]);
>     A.endindices();
> 
>     for (size_type i = 0; i < nnz; i++)
>         A[row[i]][col[i]] = entry[i];
> }
> 
> int main(int argc, char** argv)
> {
>   try
>   {
>     typedef Dune::FieldMatrix<double,1,1> M;
>     typedef Dune::BCRSMatrix<M> Matrix;
>     Matrix A;
>     
>     readBCRSMatrixFromMatlabFile(argv[1], A);
> 
>     typedef Dune::FieldVector<double, 1> VB;
>     typedef Dune::BlockVector<VB> Vector;
>     Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);        // make linear operator from A
>     Dune::SeqILU0<Matrix,Vector,Vector> prec(A, 1.0);        // preconditioner object
> 
>     int n = A.N();
> 
>     Vector b(n);
>     b = 0;
>     b[0] = -1e+05;
>     b[1] = -1e+05;
>     b[2] = -1e+05;
>     b[3] = -1e+05;
>     b[4] = -1e+05;
>     b[5] = -1e+05;
>     b[6] = -1e+05;
>     b[7] = -1e+05;
>     b[8] = -1e+05;
> 
>     Vector x(n);
>     x = 0;
> 
>     Dune::InverseOperatorResult r;
> 
>     double reduction = 1e-14;
>     int maxit = 500;
>     int verbose = 2;
>     Dune::LoopSolver<Vector> loop(op, prec, reduction, maxit, verbose); // an inverse operator
>     Vector bTemp(b);
>     loop.apply(x, bTemp, r);
>     std::cout << "LoopSolver finished.\n";
> 
>     Dune::BiCGSTABSolver<Vector> bicgstab(op, prec, reduction, maxit, verbose); // an inverse operator
>     x = 0;
>     bicgstab.apply(x, b, r);
>     std::cout << "BiCGSTABSolver finished.\n";
> 
>     return 0;
>   }
>   catch (Dune::Exception &e){
>     std::cerr << "Dune reported error: " << e << std::endl;
>   }
>   catch (...){
>     std::cerr << "Unknown exception thrown!" << std::endl;
>   }
> } 

> 1 1 1
> 1 2 0
> 1 10 0
> 1 11 0
> 2 1 0
> 2 2 1
> 2 3 0
> 2 10 0
> 2 11 0
> 2 12 0
> 3 2 0
> 3 3 1
> 3 4 0
> 3 11 0
> 3 12 0
> 3 13 0
> 4 3 0
> 4 4 1
> 4 5 0
> 4 12 0
> 4 13 0
> 4 14 0
> 5 4 0
> 5 5 1
> 5 6 0
> 5 13 0
> 5 14 0
> 5 15 0
> 6 5 0
> 6 6 1
> 6 7 0
> 6 14 0
> 6 15 0
> 6 16 0
> 7 6 0
> 7 7 1
> 7 8 0
> 7 15 0
> 7 16 0
> 7 17 0
> 8 7 0
> 8 8 1
> 8 9 0
> 8 16 0
> 8 17 0
> 8 18 0
> 9 8 0
> 9 9 1
> 9 17 0
> 9 18 0
> 10 1 6.25e-06
> 10 2 -3.125e-05
> 10 10 0.0001875
> 10 11 -0.0001375
> 10 19 6.25e-06
> 10 20 -3.125e-05
> 11 1 -3.125e-05
> 11 2 1.25e-05
> 11 3 -3.125e-05
> 11 10 -0.0001375
> 11 11 0.000375
> 11 12 -0.0001375
> 11 19 -3.125e-05
> 11 20 1.25e-05
> 11 21 -3.125e-05
> 12 2 -3.125e-05
> 12 3 1.25e-05
> 12 4 -3.125e-05
> 12 11 -0.0001375
> 12 12 0.000375
> 12 13 -0.0001375
> 12 20 -3.125e-05
> 12 21 1.25e-05
> 12 22 -3.125e-05
> 13 3 -3.125e-05
> 13 4 1.25e-05
> 13 5 -3.125e-05
> 13 12 -0.0001375
> 13 13 0.000375
> 13 14 -0.0001375
> 13 21 -3.125e-05
> 13 22 1.25e-05
> 13 23 -3.125e-05
> 14 4 -3.125e-05
> 14 5 1.25e-05
> 14 6 -3.125e-05
> 14 13 -0.0001375
> 14 14 0.000375
> 14 15 -0.0001375
> 14 22 -3.125e-05
> 14 23 1.25e-05
> 14 24 -3.125e-05
> 15 5 -3.125e-05
> 15 6 1.25e-05
> 15 7 -3.125e-05
> 15 14 -0.0001375
> 15 15 0.000375
> 15 16 -0.0001375
> 15 23 -3.125e-05
> 15 24 1.25e-05
> 15 25 -3.125e-05
> 16 6 -3.125e-05
> 16 7 1.25e-05
> 16 8 -3.125e-05
> 16 15 -0.0001375
> 16 16 0.000375
> 16 17 -0.0001375
> 16 24 -3.125e-05
> 16 25 1.25e-05
> 16 26 -3.125e-05
> 17 7 -3.125e-05
> 17 8 1.25e-05
> 17 9 -3.125e-05
> 17 16 -0.0001375
> 17 17 0.000375
> 17 18 -0.0001375
> 17 25 -3.125e-05
> 17 26 1.25e-05
> 17 27 -3.125e-05
> 18 8 -3.125e-05
> 18 9 6.25e-06
> 18 17 -0.0001375
> 18 18 0.0001875
> 18 26 -3.125e-05
> 18 27 6.25e-06
> 19 10 6.25e-06
> 19 11 -3.125e-05
> 19 19 0.0001875
> 19 20 -0.0001375
> 19 28 6.25e-06
> 19 29 -3.125e-05
> 20 10 -3.125e-05
> 20 11 1.25e-05
> 20 12 -3.125e-05
> 20 19 -0.0001375
> 20 20 0.000375
> 20 21 -0.0001375
> 20 28 -3.125e-05
> 20 29 1.25e-05
> 20 30 -3.125e-05
> 21 11 -3.125e-05
> 21 12 1.25e-05
> 21 13 -3.125e-05
> 21 20 -0.0001375
> 21 21 0.000375
> 21 22 -0.0001375
> 21 29 -3.125e-05
> 21 30 1.25e-05
> 21 31 -3.125e-05
> 22 12 -3.125e-05
> 22 13 1.25e-05
> 22 14 -3.125e-05
> 22 21 -0.0001375
> 22 22 0.000375
> 22 23 -0.0001375
> 22 30 -3.125e-05
> 22 31 1.25e-05
> 22 32 -3.125e-05
> 23 13 -3.125e-05
> 23 14 1.25e-05
> 23 15 -3.125e-05
> 23 22 -0.0001375
> 23 23 0.000375
> 23 24 -0.0001375
> 23 31 -3.125e-05
> 23 32 1.25e-05
> 23 33 -3.125e-05
> 24 14 -3.125e-05
> 24 15 1.25e-05
> 24 16 -3.125e-05
> 24 23 -0.0001375
> 24 24 0.000375
> 24 25 -0.0001375
> 24 32 -3.125e-05
> 24 33 1.25e-05
> 24 34 -3.125e-05
> 25 15 -3.125e-05
> 25 16 1.25e-05
> 25 17 -3.125e-05
> 25 24 -0.0001375
> 25 25 0.000375
> 25 26 -0.0001375
> 25 33 -3.125e-05
> 25 34 1.25e-05
> 25 35 -3.125e-05
> 26 16 -3.125e-05
> 26 17 1.25e-05
> 26 18 -3.125e-05
> 26 25 -0.0001375
> 26 26 0.000375
> 26 27 -0.0001375
> 26 34 -3.125e-05
> 26 35 1.25e-05
> 26 36 -3.125e-05
> 27 17 -3.125e-05
> 27 18 6.25e-06
> 27 26 -0.0001375
> 27 27 0.0001875
> 27 35 -3.125e-05
> 27 36 6.25e-06
> 28 19 6.25e-06
> 28 20 -3.125e-05
> 28 28 0.0001875
> 28 29 -0.0001375
> 28 37 6.25e-06
> 28 38 -3.125e-05
> 29 19 -3.125e-05
> 29 20 1.25e-05
> 29 21 -3.125e-05
> 29 28 -0.0001375
> 29 29 0.000375
> 29 30 -0.0001375
> 29 37 -3.125e-05
> 29 38 1.25e-05
> 29 39 -3.125e-05
> 30 20 -3.125e-05
> 30 21 1.25e-05
> 30 22 -3.125e-05
> 30 29 -0.0001375
> 30 30 0.000375
> 30 31 -0.0001375
> 30 38 -3.125e-05
> 30 39 1.25e-05
> 30 40 -3.125e-05
> 31 21 -3.125e-05
> 31 22 1.25e-05
> 31 23 -3.125e-05
> 31 30 -0.0001375
> 31 31 0.000375
> 31 32 -0.0001375
> 31 39 -3.125e-05
> 31 40 1.25e-05
> 31 41 -3.125e-05
> 32 22 -3.125e-05
> 32 23 1.25e-05
> 32 24 -3.125e-05
> 32 31 -0.0001375
> 32 32 0.000375
> 32 33 -0.0001375
> 32 40 -3.125e-05
> 32 41 1.25e-05
> 32 42 -3.125e-05
> 33 23 -3.125e-05
> 33 24 1.25e-05
> 33 25 -3.125e-05
> 33 32 -0.0001375
> 33 33 0.000375
> 33 34 -0.0001375
> 33 41 -3.125e-05
> 33 42 1.25e-05
> 33 43 -3.125e-05
> 34 24 -3.125e-05
> 34 25 1.25e-05
> 34 26 -3.125e-05
> 34 33 -0.0001375
> 34 34 0.000375
> 34 35 -0.0001375
> 34 42 -3.125e-05
> 34 43 1.25e-05
> 34 44 -3.125e-05
> 35 25 -3.125e-05
> 35 26 1.25e-05
> 35 27 -3.125e-05
> 35 34 -0.0001375
> 35 35 0.000375
> 35 36 -0.0001375
> 35 43 -3.125e-05
> 35 44 1.25e-05
> 35 45 -3.125e-05
> 36 26 -3.125e-05
> 36 27 6.25e-06
> 36 35 -0.0001375
> 36 36 0.0001875
> 36 44 -3.125e-05
> 36 45 6.25e-06
> 37 28 0
> 37 29 0
> 37 37 1
> 37 38 0
> 38 28 0
> 38 29 0
> 38 30 0
> 38 37 0
> 38 38 1
> 38 39 0
> 39 29 0
> 39 30 0
> 39 31 0
> 39 38 0
> 39 39 1
> 39 40 0
> 40 30 0
> 40 31 0
> 40 32 0
> 40 39 0
> 40 40 1
> 40 41 0
> 41 31 0
> 41 32 0
> 41 33 0
> 41 40 0
> 41 41 1
> 41 42 0
> 42 32 0
> 42 33 0
> 42 34 0
> 42 41 0
> 42 42 1
> 42 43 0
> 43 33 0
> 43 34 0
> 43 35 0
> 43 42 0
> 43 43 1
> 43 44 0
> 44 34 0
> 44 35 0
> 44 36 0
> 44 43 0
> 44 44 1
> 44 45 0
> 45 35 0
> 45 36 0
> 45 44 0
> 45 45 1

> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune