Dear Gael,

I'm working on the development of an incompressible Navier Stokes solver and I am using
Eigen for the solution of the resulting pressure-Poisson equation. I am mostly interested on iterative solvers
and more particularly BiCGSTAB since my LHS is non-symmetric.

Until now I was using an in-house code based on the implementation of the BiCGSTAB algorithm that can be found
on the "Numerical Recipes in C" book and I was not facing any problems with respect to convergence of my problem etc.

Now, since I want to have an optimized multi-threaded version of my code, I chose Eigen's BiCGSTAB version in order to solve
the Poisson equation but the code does not converge, not even for tolerance of 1e-4.

I tried both diagonal and ILU preconditioner, also tuned them, but no improvement was observed.

Regarding the LHS construction I am following the same strategy on both codes.
In fact, before solving the pressure-Poisson equation I am solving another linear system, of the same size, that is though very diagonally dominant
and it converges in a couple of iterations. When I compare the output of the two BiCGSTAB solvers I get almost identical results(difference 1e-6)

Do you have any idea why the solver does not converge in the case that I am using Eigen's BiCGSTAB?
Any hints on the method I could follow in order to find what's wrong?

Kind Regards,

Mike

Dear Gael,

It seems that we were using the L1 norm in our code while the squaredNorm() is being used in the Eigen's implementation.
I tried to use the lpNorm<1>() in the BiCGSTAB of EIGEN but I am getting the following error

/usr/local/include/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h:104:28: error: invalid operands of types ‘<unresolved overloaded function type>’ and ‘int’ to binary ‘operator<’
tol_error = sqrt(r.lpNorm<1>()/rhs_sqnorm);


Is it possible to change the norm in the BiCGSTAB implementation of EIGEN? If yes can you please explain me how?

Kind Regards,

Mike

You should add the template keyword because you are in a template code:

r.template lpNorm<1>()

I'll compare our implementation of Numerical Recipes's one later.

btw, it would be nice if you could share your problematic matrix and right-hand-side so that we can include it in our regression tests. See the lines of code to copy-pasted: http://eigen.tuxfamily.org/dox-devel/gr ... tml#title7

Dear Gael,

Thanks for your reply,

I did use the L1 norm to the Eigen's BICGSTAB implementation but the results I am getting are even worse.
The linear system that I was using before the pressure-poisson, the one that the diagonal term is much stronger,
does not ever converge now while before it was performing 3-4 iterations and then converge. In addition, the pressure-poisson returns a NaN as an error.

I was thinking that the reinitialization of the solvWithGuess technique could play a role in the different behaviour of the Eigen's BiCGSTAB solver.
We are not performing this action to our implementation. What do you think? Do you have any other propositions with respect to the way that I could examine the causes of this problem?

Btw I sent the LHS and RHS of my problem to your email address.

Kind Regards,
Mike

Could you give the L2 error( A*x-b).norm()/b.norm() you get with the initial solver ? Looking at the different iteration, using ILUT, it quickly reach a relative error of 0.0005, to then start diverging....

ok, I've tried the NR's linbcg routine with Jacobi preconditioner (diagonal part), and itol==4 (L-inf norm on x stopping criteria), and 1e-4 for the tolerance. It stops after 734 iterations, with a reported error of 7.62-05, but if I compute the L2 residual, I get a huge error: 1.9. Do you observe the same? If so, then the "solution" you've obtained might not be as good as you thought.

Dear Gael,

I observe exactly the same behaviour as the one you said with the ILU preconditioner. It converges for a tolerance of 1e-5 for a couple of time
steps but then the solution of the linear system diverges.

With respect to the linbcg solver, we are only using the L1 norm because it was the one having the
best ratio with respect to computational time and accuracy. Honestly, I didn't try any other convergence criteria recently, I ll give them a try as soon as I can and I ll let
you know about the results.

Based on your experience, do you think that a our convergence issue is related to the nature of our problem, incompressible Navier Stokes pressure poisson equation, or you think that we have a buggy LHS or RHS implementation?

Thanks for your help,
Mike

I'm confused about linbcg: the initial implementation provides 4 criterion, but none of them is based on L1. Moreover, using the L1 or L2 norm for the stopping criterion should not make much difference. Also, I did not managed to get any good solution with linbcg for your problem.

Then, the BICGSTAB algorithm in general is known to exhibit breakdown for some kind of problems, typically for problems having eigenvalues with large imaginary parts. Detecting them is rather difficult and kind of pointless without a remedy... maybe we could keep track of the best solution and residual norm found during the iterations, and stop if the current residual becomes much larger (say x10^6) and return the best solution which might still be ok in some cases.

Dear Gael,

The convergence criterion based on the L1 norm was introduced by a former member of our group. I contacted him to ask him why he
was using that but did not receive a reply yet. The previous version of the code though, the one using the linbcg version of the algorithm was validated with several benchmark turbulent fluid flow problems and based on this fact I can say that the "solution" of the pressure-poisson problem is at least as good as it should be for my needs.

I am afraid that my problem is not set-up correctly and more precisely I have a suspicion that the Matrix I created with the Triplets technique is wrong.
In the serial version of my code I am following the PCGPACK sparse storage convention which creates two vectors s_A and ij_A. To constuct these vectors I followed the guidelines given on the 2nd edition of the Numerical Recipes book. Is there a way to create an EIGEN sparse matrix object by providing these two vectors? If not, can you propose a way to compare the two LHS created by the two methods in order to see if they resemble the same LHS matrix?

In any case, I would like to sincerely thank you for your fast replies and for your efforts you make to help me.

Kind Regards,

Mike