NaN reaction rates for some reversible reactions

[roystgnr]

This can be triggered in our kinetics_unit tests, when air_5sp.xml is evaluated at T=500K, using StateType=float, but I fear the problem is more general; I encountered it while trying to debug a double precision code.

The equilibrium constant for reaction 0001 in this unit test is so close to zero (exppower is -213.401) that it gets rounded down to 0. kbkwd_times_products = kfwd/Keq is then NaN, which propagates through the evaluations.

I'm not sure what the proper fix is. The correct interpretation of these results is "N should recombine to N2 ludicrously fast", right? We could test for Keq=0, but then what? Return numeric_limits::max or numeric_limits::infinity, and hope the user's code is good at handling really stiff integration?

Guess how those tests were passing before I started throwing in assert(!isnan)? omega_dot was NaN, which means sum was NaN, which means abs(sum) was NaN, which means "abs(sum) > sum_tol" returns false, because all comparisons with NaN return false. We needed to be doing "if (!(error < tol)) fail" instead of "if (error > tol) fail" to be safe.

[SylvainPlessis]

The correct interpretation of these results is "N should recombine to N2 ludicrously fast", right?

This, and "N2 never dissociates". Actually, my view on this is that the reaction

2 N + M -> N2 + M

is irreversible and instantaneous (hence the arrow), which means writing it as they do in the file (N2 + M [=] 2 N + M) has no sense...

My personal view on how to treat it is to actually reverse the reaction and make it instantaneous and irreversible, but that'd be much more (quickly seen from a 2am perspective) than just a fix.