Self-Gravity error: "Bicgstab:: beta is nan in step 2"

[Anto6453]

Hi all,
I'm trying to start a simulation from a previous dump file with a higher resolution. But when I increase the resolution too much, a problem occurs with the calculation of the self-gravity field.
I'm using PBICGSTAB to calculate the self-gravity field, and it gives me the error "Bicgstab:: beta is nan in step 2". I tried to check the reason by defining DEBUG_BICGSTAB in the definition.hpp file, and it appears that at some point the "rho" parameter becomes zero, which means that "beta" becomes nan in the next step due to division by zero.
So I wanted to know if there was a way of avoiding this (it also happens with BICGSTAB without a preconditioner, and MINRES).

Also, when I use PMINRES, CG or PCG, I reach the maximum iteration, and as I'm restarting from a dump file, I was wondering if it was possible to access the old gravity field in the dump file, so that I could use it to define the new field close to the correct solution on the new grid (which would make convergence faster, I suppose).

Thanks in advance for the help.

[glesur]

Hi @Anto6453, unfortunately there is no way to guarantee that Bicgstab will always to converge without any trouble, especially when your guess is far from the the solution. With Jacobi, you're sure you will get convergence, but then it is very slow. It might be a good idea however to do the first guess with Jacobi, and then revert to Bicgstab.

On the number of iterations, this can be changed in your input file using the parameter maxIter in the [SelfGravity] block.

Finally, it is possible to add the potential in the restart dump to get faster convergence. I've just opened a PR (#245 ) with that specific modification which I might merge in the upcoming version if you think it improves the situation.

[Anto6453]

Hi @glesur.
Thank you for your answer!

I tried to use Jacobi but it take too much time to converge, and even by increasing a lot the maxIter parameter, it still not converge.
With the self-gravity potential in the restart dump (#245) , it seems that it converge a lot faster when doing a normal restart with the -restart option, so this could probably solve my problem I think.

The new problem I have now is that if i want to change the grid to increase the resolution, I need to set manually the potential in the setup.cpp file from the image of the dump file. As the self-Gravity potential is not contained inside the dataBlockHost class, I tried to add it by myself using replace_idefix_source in the CMakeLists.txt file to modify the dataBlocklockHost.cpp/hpp and selfGravity.cpp/hpp files. But it seems that the modification in the .hpp files I make are not take into account during the compilation. So I was wondering, does the function replace_idefix_source
is supposed to work with .hpp files, and if yes should I open an issue?

[glesur]

You can't replace a .hpp file with replace_idefix_source, that's a limitation of the way c++ compilers work, there is no proper workaround, but to modify directly the idefix source tree. I've removed replace_idefix_source from the documentation, as I think we should progressively move away from this.

[glesur]

@Anto6453 would be useful to know a bit more about the fix you told me about to avoid singularities when bicgstab is iterating. If you have a reference or a link, let me know!

[Anto6453]

Hi @glesur,
The main problem I encountered with BICGSTAB is that during its iteration, the variable rho = ComputeDotProduct(res0, res) can become equal to zero, so that in the next step the variable beta = rho/rhoOld becomes nan.

When this happened, the solver would restart from a flat potential that was too far from the solution to converge correctly. I therefore followed the recommendation of this Lecture Note (ln07.pdf) which explains on page 6 that when rho becomes equal to (or very close to) zero, it is possible to set only res0 and dir to res (as when we initialize the solver) and continue the iteration without restarting the potential.

This works well enough in my case to avoid the BICGSTAB crashes I had in my simulation.