Residual value goes to NaN while solving a system of nonlinear equations

Question

I am solving a system of coupled nonlinear equations using Newton's method, similar to $$\begin{split} c_A(A, B)\partial_tA&=\nabla\left(k_A(A, B)\nabla A\right) + f_A(A, B, t)\\ c_B(A, B)\partial_tB&=\nabla\left(k_B(A, B)\nabla B\right) + f_B(A, B, t) \end{split}$$ with $f_x(), k_x()\text{ and }c_x()$ nonlinear functions depending on $A$ and $B$.

For calculating the Jacobian matrix I am using an AD-tool (Sacado/ADOL-C).
When calculating just one equation at a time (and setting the other equation to a constant value), the program works fine. When calculating both equations at the same time, sooner or later the solver will fail in iteration 0 (GMRES-solver) with "residual is equal to NaN". This also will happen if I reduce one of the equations to (with $B$ as an example): $$c_B(A, B)\partial_tB=0$$

What would be the best approach for debugging this problem? I assume it is not related to a race condition, due to always coming up at the same time, regardless of the number of threads used.

One idea I am testing is to use AD for generating the Jacobian for the first equation and add the Jacobian of the second equation by hand.
Are there other reasons for this behaviour?

Answer 1

Several possible other leads:

• check your RHS (size, presence of NaNs or infs) for the case when you solve equations simultaneously. There is a chance that something goes wrong on that end.
• compile your code in a Debug mode and run with 1 thread for a small example and see how it goes.
• Valgrind to detect potential access to unallocated memory (and memory leaks as well) to see potential issues mostly on your side (unless you can mentally power through valgrinding an external library)
• another advanced option would be to try seeing at which step those NaNs occur in the first place. Take a look at this post describing catching NaNs with gcc/gdb in C and C++ using a debugger as well as a classic runtime signaling NaNs one.