Self-consistent NSE floating-point overflow

[yut23]

test_nse inputs:

unit_test.small_temp = 1.e5
unit_test.small_dens = 1.e5

unit_test.density = 2.e8
unit_test.temperature = 3.9999999999999981e9

unit_test.X1 = 1.e-10
unit_test.X2  = 0.99999999839999998
unit_test.X3  = 1.e-10
unit_test.X4  = 1.e-10
unit_test.X5  = 1.e-10
unit_test.X6  = 1.e-10
unit_test.X7  = 1.e-10
unit_test.X8  = 1.e-10
unit_test.X9  = 1.e-10
unit_test.X10 = 1.e-10
unit_test.X11 = 1.e-10
unit_test.X12 = 1.e-10
unit_test.X13 = 1.e-10
unit_test.X14 = 1.e-10
unit_test.X15 = 1.e-10
unit_test.X16 = 1.e-10
unit_test.X17 = 1.e-10
unit_test.X18 = 1.e-10

unit_test.mu_p = -5.0
unit_test.mu_n = -10.0

unit_test.y_e = 0.5

Run as:

make NETWORK_DIR=subch_base SCREEN_METHOD=chabrier1998 DEBUG=TRUE
gdb --args ./main3d.gnu.DEBUG.ex inputs_subch_base.fpe amrex.fpe_trap_{zero,overflow,invalid}=1

Inspection in gdb:

Program received signal SIGFPE, Arithmetic exception.
0x000000000041b8b9 in dogleg<2> (r=..., diag=..., qtb=..., delta=7.8803954902403048e+159, x=..., wa1=..., wa2=...) at ../../util/hybrj/hybrj_dogleg.H:106
106                 wa1(i) += r(l) * temp;
(gdb) bt
#0  0x000000000041b8b9 in dogleg<2> (r=..., diag=..., qtb=..., delta=7.8803954902403048e+159, x=..., wa1=..., wa2=...) at ../../util/hybrj/hybrj_dogleg.H:106
#1  0x0000000000415b50 in hybrj<2, nse_solver_data<burn_t> > (hj=..., data=...) at ../../util/hybrj/hybrj.H:187
#2  0x00000000004114db in nse_hybrid_solver<burn_t> (state_data=..., eps=1e-10) at ../../nse_solver/nse_solver.H:314
#3  0x000000000040fe21 in get_actual_nse_state<burn_t> (state=..., eps=1e-10, input_ye_is_valid=true) at ../../nse_solver/nse_solver.H:530
#4  0x00000000004099d6 in nse_example_c () at ./nse_example.H:69
#5  0x0000000000409c18 in main (argc=6, argv=0x7fffffffd398) at main.cpp:31
(gdb) p l
$1 = 1
(gdb) p r
$2 = (amrex::Array1D<double, 1, 3> &) @0x7fffffffbf90: {arr = {-1.2593613952543642e+161, -1.2593613952543639e+161, -0}}
(gdb) p r.arr[0]
$3 = -1.2593613952543642e+161
(gdb) p temp
$4 = -1.5503337516169112e+159

Increasing the temperature by 1 ULP to 3.9999999999999986e9 or decreasing by 5 ULP to 3.9999999999999957e9 avoids the overflow.

These inputs come from AMReX-Astro/Castro#2768, in the get_actual_nse_state() call in problem_initialize_state_data.H.

[yut23]

Note: updating X2 to 0.9999999983 to fix the normalization (re: AMReX-Astro/Castro#2806) gives the same behavior.