I use the parameter
TwoPunctures::TP_epsilon = 1e-6
to ensure that the NaN never appears, no matter where I place the BHs. This replaces a denominator
sqrt (pow (r2_plus, 2))
with
sqrt (pow (r2_plus, 2) + pow (TP_epsilon, 4))
See
so adds a fixed source of error to your evolution, but the effect of this is negligible.