Self-consistent NSE floating-point overflow
Opened this issue · 1 comments
yut23 commented
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 commented
Note: updating X2 to 0.9999999983 to fix the normalization (re: AMReX-Astro/Castro#2806) gives the same behavior.