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

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$$.

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?

• 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.
• 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.