Cubic NSE interpolation can turn neutrino cooling into heating

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

Start with nse_tabular/nse_table.H:468-472 and the SDC source-use sites in integration/nse_update_sdc.H:89-90 and :155-156. Reproduce the supplied aprox19 case with network.nse_table_interp_linear = 0, then verify the chosen interpolation or fallback keeps nonnegative tabulated losses nonnegative. Add the stated in-bounds state as a regression test and confirm the linear result remains valid.