Cubic NSE interpolation can turn neutrino cooling into heating

[zingale]

Reviewed commit: 3ac0bda4f00da3eed0446beadfc1cd5c51b99a40

Suggested priority: High

Affected code: nse_tabular/nse_table.H:468-472; integration/nse_update_sdc.H:89-90 and integration/nse_update_sdc.H:155-156.

Description

The default cubic interpolant does not preserve positivity of e_nu. Positive neutrino-loss data can produce a negative interpolated loss inside the table. The SDC energy source subtracts this value, so the weak-neutrino contribution becomes heating. This does not imply that the total energy source, including other contributions, is positive.

Reproduction

Use the shipped aprox19 table and network.nse_table_interp_linear = 0. Call nse_interp() with:

T   = 3667943011.533362 K
rho = 9434630294.888599 g/cm^3
Ye  = 0.433159096449449

The result is e_nu = -4.353326581794286e12 erg/g/s. All 64 stencil values are positive (minimum 1.3703405e11). Linear interpolation at the same state returns 9.761675725144252e12.

Expected behavior / suggested fix

A nonnegative tabulated loss should remain nonnegative. Use positivity-preserving interpolation or a local fallback to linear interpolation when cubic interpolation violates this constraint. Add this in-bounds state as a regression test.

Validation: Executed against production interpolation code; see reproduce.cpp and reproduction-output.txt

reproduce.cpp
reproduction-output.txt

[zingale]

The easiest fix is to replace the 1D cubic interpolation routine we currently use with the monotone PCHIP method:

F. N. Fritsch and R. E. Carlson, “Monotone Piecewise Cubic Interpolation,” SIAM Journal on Numerical Analysis 17, 238–246 (1980).
doi:10.1137/0717021

This is what we've switched the partition functions over to.

We would then do this in each direction independently. So it will not be formally monotone with the cross-terms, but it should be much better. We can still add a check.

[zingale]

I think I mostly have this fixed now. I want to do a bit of cleaning and have the option for species to use trilinear.