diff --git a/.github/workflows/test_nse_interp.yml b/.github/workflows/nse_interp_check.yml similarity index 74% rename from .github/workflows/test_nse_interp.yml rename to .github/workflows/nse_interp_check.yml index d372a0381e..5944361307 100644 --- a/.github/workflows/test_nse_interp.yml +++ b/.github/workflows/nse_interp_check.yml @@ -26,22 +26,22 @@ jobs: sudo apt-get update -y -qq sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0 - - name: Compile, test_nse_interp (NSE, aprox19) + - name: Compile, nse_interp_check (NSE, aprox19) run: | - cd unit_test/test_nse_interp + cd unit_test/nse_interp_check make realclean make -j 4 - - name: Run test_nse_interp (NSE, aprox19) + - name: Run nse_interp_check (NSE, aprox19) run: | - cd unit_test/test_nse_interp + cd unit_test/nse_interp_check ./main3d.gnu.ex amrex.fpe_trap_{invalid,zero,overflow}=1 > test.out - name: Print backtrace - if: ${{ failure() && hashFiles('unit_test/test_nse_interp/Backtrace.0') != '' }} - run: cat unit_test/test_nse_interp/Backtrace.0 + if: ${{ failure() && hashFiles('unit_test/nse_interp_check/Backtrace.0') != '' }} + run: cat unit_test/nse_interp_check/Backtrace.0 - name: Compare to stored output (NSE, aprox19) run: | - cd unit_test/test_nse_interp + cd unit_test/nse_interp_check diff -I "^Initializing AMReX" -I "^AMReX" -I "^reading in reaclib rates" test.out ci-benchmarks/aprox19.out diff --git a/Docs/source/unit_tests.rst b/Docs/source/unit_tests.rst index 675399a4c8..7ef7cd83a0 100644 --- a/Docs/source/unit_tests.rst +++ b/Docs/source/unit_tests.rst @@ -174,7 +174,7 @@ One-zone tests Infrastructure tests ==================== -.. index:: test_linear_algebra, test_nse_interp, test_parameters, test_sdc_vode_rhs +.. index:: test_linear_algebra, nse_interp_check, test_parameters, test_sdc_vode_rhs * ``test_linear_algebra`` : @@ -182,9 +182,9 @@ Infrastructure tests to get $b = Ax$, and then call the linear algebra routines to see if we we recover $x$ from $b$. -* ``test_nse_interp`` : +* ``nse_interp_check`` : - run various tests of the NSE interpolation routines. + run various tests of the tabular NSE interpolation routines. * ``test_parameters`` : diff --git a/unit_test/test_nse_interp/GNUmakefile b/unit_test/nse_interp_check/GNUmakefile similarity index 100% rename from unit_test/test_nse_interp/GNUmakefile rename to unit_test/nse_interp_check/GNUmakefile diff --git a/unit_test/test_nse_interp/Make.package b/unit_test/nse_interp_check/Make.package similarity index 100% rename from unit_test/test_nse_interp/Make.package rename to unit_test/nse_interp_check/Make.package diff --git a/unit_test/test_nse_interp/README.md b/unit_test/nse_interp_check/README.md similarity index 68% rename from unit_test/test_nse_interp/README.md rename to unit_test/nse_interp_check/README.md index 383037be55..f11991c5fa 100644 --- a/unit_test/test_nse_interp/README.md +++ b/unit_test/nse_interp_check/README.md @@ -1,4 +1,4 @@ -# test_nse_interp +# `nse_interp_check` This is a simple test of the NSE cubic interpolation. @@ -9,7 +9,10 @@ Then it does 1-d interpolation in each direction (rho, T, Ye) to make sure that the interpolation there works as expected. This is done just for Abar. -Finally, it calls the full interface that does tricubic interpolation +Next, it calls the full interface that does tricubic interpolation and prints out the interpolated state. +Finally, it chooses a zone where cubic interpolation can violate +monotonicity and interpolates neutrino energy. + This is for the tabular NSE: `USE_NSE_TABLE=TRUE` diff --git a/unit_test/test_nse_interp/_parameters b/unit_test/nse_interp_check/_parameters similarity index 100% rename from unit_test/test_nse_interp/_parameters rename to unit_test/nse_interp_check/_parameters diff --git a/unit_test/nse_interp_check/ci-benchmarks/aprox19.out b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out new file mode 100644 index 0000000000..53a324c615 --- /dev/null +++ b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out @@ -0,0 +1,128 @@ +Initializing AMReX (26.09-119-ge60cdc18711c)... +AMReX (26.09-119-ge60cdc18711c) initialized +starting the single zone burn... +reading the NSE table (C++) ... + +rho, T, Ye = 1230000000 5180000000 0.472 + +density value brackets: 9.05 < 9.089905111 < 9.1 +temperature value brackets: 9.71 < 9.71432976 < 9.72 +Ye value brackets: 0.4725 < 0.472 < 0.47 + +testing 1-d cubic interpolation at input state + + 4 rho values: + 9.000 55.5223665 + 9.050 55.57080116 + 9.100 55.61707675 + 9.150 55.66151541 + + cubic interpolated value: 55.6078922 + + 4 T values: + 9.700 55.5223665 + 9.710 55.39489494 + 9.720 55.22387568 + 9.730 54.99343581 + + cubic interpolated value: 55.32712406 + + 4 Ye values: + 0.475 55.5223665 + 0.472 55.65227675 + 0.470 55.78142906 + 0.468 55.90244371 + + cubic interpolated value: 55.678404 + + +tricubic interpolated values at input state + + abar = 55.60652462 + bea = 8.755119523 + dyedt = -0.004527880871 + dbeadt = 0.000111945508 + e_nu = 1.101641642e+16 + + X(H1) = 0 + X(He3) = 2.723457525e-14 + X(He4) = 0.0003164137294 + X(C12) = 7.975878903e-10 + X(N14) = 3.102372259e-14 + X(O16) = 1.998161464e-09 + X(Ne20) = 2.444604777e-11 + X(Mg24) = 1.029036116e-08 + X(Si28) = 8.173012779e-06 + X(S32) = 1.194076593e-05 + X(Ar36) = 1.243414613e-05 + X(Ca40) = 1.259531492e-05 + X(Ti44) = 0.0001438001413 + X(Cr48) = 0.01146726617 + X(Fe52) = 0.07746625779 + X(Fe54) = 0.9098424639 + X(Ni56) = 0.0007087774945 + X(n) = 2.078852742e-08 + X(p) = 9.843630258e-06 + + +testing temperature derivatives of cubic + + first finite-difference derivatives + dAbar/dT = -1.072562355e-09 + dbea/dT = -6.867059575e-12 + + now using derivative of the interpolant + dAbar/dT = -1.072562604e-09 + dbea/dT = -6.867048589e-12 + +testing density derivatives of cubic + + first finite-difference derivatives + dAbar/drho = 3.987525411e-10 + dbea/drho = 7.619559936e-13 + + now using derivative of the interpolant + dAbar/drho = 3.987522836e-10 + dbea/drho = 7.618831514e-13 + + +EOS T from e consistency check (e should remain unchanged) + + old method: invert EOS without considering NSE changes: + change in e: 1.395278886e+18 1.38844906e+18 + + NSE method: use the nse_T_abar_from_e() function: + updated T: 6394534499 + change in abar: 55.60652462 50.26831386 + change in e: 1.38844906e+18 1.38844906e+18 + +EOS T from p consistency check (p should remain unchanged) + + old method: invert EOS without considering NSE changes: + change in p: 6.622159603e+26 6.577850616e+26 + + NSE method: use the nse_T_abar_from_p() function: + updated T: 6466757500 + change in abar: 55.60652462 49.50320619 + change in p: 6.577850616e+26 6.577850616e+26 + +EOS rho from p consistency check (p should remain unchanged) + + old method: invert EOS without considering NSE changes: + change in p: 6.577758474e+26 6.577850616e+26 + + NSE method: use the nse_rho_abar_from_p() function: + updated T: 5180000000 + change in abar: 55.60652462 55.62494615 + change in p: 6.577850616e+26 6.577850616e+26 + + +explore a troublesome state + + rho = 9.434630294888599e+09, T = 3.667943011533362e+09, Ye = 0.433159096449449 + indices of lowest corner of 4x4x4 cube: ir0 = 58, it0 = 16, ic0 = 26 + + cubic interpolation gives: e_nu = -4.353326582e+12 + bounds of e_nu in the 64 cells used for tricubic: 1.3703405e+11 6.2562515e+14 + +AMReX (26.09-119-ge60cdc18711c) finalized diff --git a/unit_test/test_nse_interp/main.cpp b/unit_test/nse_interp_check/main.cpp similarity index 100% rename from unit_test/test_nse_interp/main.cpp rename to unit_test/nse_interp_check/main.cpp diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H similarity index 57% rename from unit_test/test_nse_interp/nse_cell.H rename to unit_test/nse_interp_check/nse_cell.H index c1eb3cf3f1..3097750bdc 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/nse_interp_check/nse_cell.H @@ -1,6 +1,9 @@ #ifndef NSE_CELL_H #define NSE_CELL_H +#include +#include + #include #include #include @@ -15,10 +18,12 @@ AMREX_INLINE void nse_cell_c() { + std::cout << std::endl; + std::cout << "rho, T, Ye = " << unit_test_rp::density << " " << unit_test_rp::temperature << " " - << unit_test_rp::ye << std::endl; + << unit_test_rp::ye << std::endl << std::endl; // check the indices @@ -43,9 +48,11 @@ void nse_cell_c() std::cout << std::endl; - // now try 1-d interpolation in T -- for cubic interpolation, we + // now try 1-d interpolation in rho, T, and Ye -- for cubic interpolation, we // need to offset the indices one to the left + std::cout << "testing 1-d cubic interpolation at input state\n\n"; + ir0--; ir0 = amrex::max(1, amrex::min(nse_table_size::nden-3, ir0)); @@ -55,25 +62,16 @@ void nse_cell_c() ic0--; ic0 = amrex::max(1, amrex::min(nse_table_size::nye-3, ic0)); + { - amrex::Real rhos[] = {nse_table_logrho(ir0), - nse_table_logrho(ir0+1), - nse_table_logrho(ir0+2), - nse_table_logrho(ir0+3)}; - - amrex::Real Ts[] = {nse_table_logT(it0), - nse_table_logT(it0+1), - nse_table_logT(it0+2), - nse_table_logT(it0+3)}; - - amrex::Real yes[] = {nse_table_ye(ic0), - nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + // density interpolation + amrex::Real rhos[] = {nse_table_logrho(ir0), + nse_table_logrho(ir0+1), + nse_table_logrho(ir0+2), + nse_table_logrho(ir0+3)}; - { - std::cout << "4 rho values: " << std::endl; + std::cout << " 4 rho values: " << std::endl; amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0+1, it0, ic0)), @@ -81,16 +79,24 @@ void nse_cell_c() nse_table::abartab(nse_idx(ir0+3, it0, ic0))}; for (int ii = 0; ii < 4; ++ii) { - std::cout << rhos[ii] << " " << _d[ii] << std::endl; + std::cout << std::format(" {:5.3f} {:15.10g}\n", rhos[ii], _d[ii]); } std::cout << std::endl; - std::cout << "cubic interpolated value: " << + std::cout << " cubic interpolated value: " << cubic(rhos, _d, nse_table_size::dlogrho, logrho) << std::endl << std::endl; } { - std::cout << "4 T values: " << std::endl; + + // temperature interpolation + + amrex::Real Ts[] = {nse_table_logT(it0), + nse_table_logT(it0+1), + nse_table_logT(it0+2), + nse_table_logT(it0+3)}; + + std::cout << " 4 T values: " << std::endl; amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0, it0+1, ic0)), @@ -98,17 +104,25 @@ void nse_cell_c() nse_table::abartab(nse_idx(ir0, it0+3, ic0))}; for (int ii = 0; ii < 4; ++ii) { - std::cout << Ts[ii] << " " << _d[ii] << std::endl; + std::cout << std::format(" {:5.3f} {:15.10g}\n", Ts[ii], _d[ii]); } std::cout << std::endl; - std::cout << "cubic interpolated value: " << + std::cout << " cubic interpolated value: " << cubic(Ts, _d, nse_table_size::dlogT, logT) << std::endl << std::endl; } { - std::cout << "4 Ye values: " << std::endl; + + // Ye interpolation + + amrex::Real yes[] = {nse_table_ye(ic0), + nse_table_ye(ic0+1), + nse_table_ye(ic0+2), + nse_table_ye(ic0+3)}; + + std::cout << " 4 Ye values: " << std::endl; amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0, it0, ic0+1)), @@ -116,16 +130,18 @@ void nse_cell_c() nse_table::abartab(nse_idx(ir0, it0, ic0+3))}; for (int ii = 0; ii < 4; ++ii) { - std::cout << yes[ii] << " " << _d[ii] << std::endl; + std::cout << std::format(" {:5.3f} {:15.10g}\n", yes[ii], _d[ii]); } std::cout << std::endl; - std::cout << "cubic interpolated value: " << + std::cout << " cubic interpolated value: " << cubic(yes, _d, -nse_table_size::dye, unit_test_rp::ye) << std::endl << std::endl; } - std::cout << "tricubic interpolated values: " << std::endl; + std::cout << std::endl; + + std::cout << "tricubic interpolated values at input state\n\n"; nse_table_t nse_state; nse_state.T = unit_test_rp::temperature; @@ -134,22 +150,24 @@ void nse_cell_c() nse_interp(nse_state); - std::cout << "abar = " << nse_state.abar << std::endl; - std::cout << "bea = " << nse_state.bea << std::endl; - std::cout << "dyedt = " << nse_state.dyedt << std::endl; - std::cout << "dbeadt = " << nse_state.dbeadt << std::endl; - std::cout << "e_nu = " << nse_state.e_nu << std::endl; + std::cout << " abar = " << nse_state.abar << std::endl; + std::cout << " bea = " << nse_state.bea << std::endl; + std::cout << " dyedt = " << nse_state.dyedt << std::endl; + std::cout << " dbeadt = " << nse_state.dbeadt << std::endl; + std::cout << " e_nu = " << nse_state.e_nu << std::endl; + std::cout << std::endl; for (int n = 0; n < NumSpec; ++n) { - std::cout << "X(" << short_spec_names_cxx[n] << ") = " << nse_state.X[n] << std::endl; + std::cout << " X(" << short_spec_names_cxx[n] << ") = " << nse_state.X[n] << std::endl; } // temp derivatives std::cout << std::endl; - std::cout << "testing temperature derivatives of cubic" << std::endl; + std::cout << std::endl; + std::cout << "testing temperature derivatives of cubic" << std::endl << std::endl; - std::cout << "first finite-difference derivatives" << std::endl; + std::cout << " first finite-difference derivatives\n"; nse_state.T = unit_test_rp::temperature; nse_state.rho = unit_test_rp::density; @@ -166,10 +184,11 @@ void nse_cell_c() nse_interp(nse_state); - std::cout << "dAbar/dT = " << (nse_state.abar - abar_old) / (nse_state.T - T_old) << std::endl; - std::cout << "dbea/dT = " << (nse_state.bea - bea_old) / (nse_state.T - T_old) << std::endl; + std::cout << " dAbar/dT = " << (nse_state.abar - abar_old) / (nse_state.T - T_old) << std::endl; + std::cout << " dbea/dT = " << (nse_state.bea - bea_old) / (nse_state.T - T_old) << std::endl; - std::cout << "now using derivative of the interpolant" << std::endl; + std::cout << std::endl; + std::cout << " now using derivative of the interpolant" << std::endl; amrex::Real dabardT = nse_interp_dT(unit_test_rp::temperature, unit_test_rp::density, unit_test_rp::ye, nse_table::abartab); @@ -177,17 +196,16 @@ void nse_cell_c() amrex::Real dbeadT = nse_interp_dT(unit_test_rp::temperature, unit_test_rp::density, unit_test_rp::ye, nse_table::beatab); - std::cout << "dAbar/dT = " << dabardT << std::endl; - std::cout << "dbea/dT = " << dbeadT << std::endl; + std::cout << " dAbar/dT = " << dabardT << std::endl; + std::cout << " dbea/dT = " << dbeadT << std::endl; std::cout << std::endl; // dens derivatives - std::cout << std::endl; - std::cout << "testing density derivatives of cubic" << std::endl; + std::cout << "testing density derivatives of cubic" << std::endl << std::endl; - std::cout << "first finite-difference derivatives" << std::endl; + std::cout << " first finite-difference derivatives" << std::endl; nse_state.T = unit_test_rp::temperature; nse_state.rho = unit_test_rp::density; @@ -201,10 +219,11 @@ void nse_cell_c() nse_interp(nse_state); - std::cout << "dAbar/drho = " << (nse_state.abar - abar_old) / (nse_state.rho - rho_old) << std::endl; - std::cout << "dbea/drho = " << (nse_state.bea - bea_old) / (nse_state.rho - rho_old) << std::endl; + std::cout << " dAbar/drho = " << (nse_state.abar - abar_old) / (nse_state.rho - rho_old) << std::endl; + std::cout << " dbea/drho = " << (nse_state.bea - bea_old) / (nse_state.rho - rho_old) << std::endl; - std::cout << "now using derivative of the interpolant" << std::endl; + std::cout << std::endl; + std::cout << " now using derivative of the interpolant" << std::endl; amrex::Real dabardrho = nse_interp_drho(unit_test_rp::temperature, unit_test_rp::density, unit_test_rp::ye, nse_table::abartab); @@ -212,9 +231,10 @@ void nse_cell_c() amrex::Real dbeadrho = nse_interp_drho(unit_test_rp::temperature, unit_test_rp::density, unit_test_rp::ye, nse_table::beatab); - std::cout << "dAbar/drho = " << dabardrho << std::endl; - std::cout << "dbea/drho = " << dbeadrho << std::endl; + std::cout << " dAbar/drho = " << dabardrho << std::endl; + std::cout << " dbea/drho = " << dbeadrho << std::endl; + std::cout << std::endl; std::cout << std::endl; // @@ -243,6 +263,8 @@ void nse_cell_c() { + std::cout << "EOS T from e consistency check (e should remain unchanged)" << std::endl << std::endl; + // first get the abar consistent with our inputs and find e nse_state.T = unit_test_rp::temperature; @@ -278,7 +300,8 @@ void nse_cell_c() eos_state.aux[iabar] = nse_state.abar; eos(eos_input_rt, eos_state); - std::cout << "EOS T from e consistency check (old method): " << eos_state.e << " " << e_new << std::endl; + std::cout << " old method: invert EOS without considering NSE changes:" << std::endl; + std::cout << " change in e: " << eos_state.e << " " << e_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -295,14 +318,16 @@ void nse_cell_c() nse_T_abar_from_e(eos_state.rho, eos_state.e, eos_state.aux[iye], eos_state.T, eos_state.aux[iabar]); - std::cout << "updated T: " << eos_state.T << std::endl; - std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; + std::cout << std::endl; + std::cout << " NSE method: use the nse_T_abar_from_e() function:" << std::endl; + std::cout << " updated T: " << eos_state.T << std::endl; + std::cout << " change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; // now check if we get back the correct e! eos(eos_input_rt, eos_state); - std::cout << "EOS T from e consistency check (new method): " << eos_state.e << " " << e_new << std::endl; + std::cout << " change in e: " << eos_state.e << " " << e_new << std::endl; std::cout << std::endl; } @@ -311,6 +336,8 @@ void nse_cell_c() // now redo it for pressure + std::cout << "EOS T from p consistency check (p should remain unchanged)" << std::endl << std::endl; + nse_state.T = unit_test_rp::temperature; nse_state.rho = unit_test_rp::density; nse_state.Ye = unit_test_rp::ye; @@ -344,7 +371,8 @@ void nse_cell_c() eos_state.aux[iabar] = nse_state.abar; eos(eos_input_rt, eos_state); - std::cout << "EOS T from p consistency check (old method): " << eos_state.p << " " << p_new << std::endl; + std::cout << " old method: invert EOS without considering NSE changes:" << std::endl; + std::cout << " change in p: " << eos_state.p << " " << p_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -361,14 +389,16 @@ void nse_cell_c() nse_T_abar_from_p(eos_state.rho, eos_state.p, eos_state.aux[iye], eos_state.T, eos_state.aux[iabar]); - std::cout << "updated T: " << eos_state.T << std::endl; - std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; + std::cout << std::endl; + std::cout << " NSE method: use the nse_T_abar_from_p() function:" << std::endl; + std::cout << " updated T: " << eos_state.T << std::endl; + std::cout << " change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; // now check if we get back the correct p! eos(eos_input_rt, eos_state); - std::cout << "EOS T from p consistency check (new method): " << eos_state.p << " " << p_new << std::endl; + std::cout << " change in p: " << eos_state.p << " " << p_new << std::endl; std::cout << std::endl; } @@ -379,6 +409,8 @@ void nse_cell_c() { + std::cout << "EOS rho from p consistency check (p should remain unchanged)" << std::endl << std::endl; + nse_state.T = unit_test_rp::temperature; nse_state.rho = unit_test_rp::density; nse_state.Ye = unit_test_rp::ye; @@ -398,7 +430,7 @@ void nse_cell_c() amrex::Real p_orig = eos_state.p; - // now perturb p and find T, then redo NSE to get abar and finally + // now perturb p and find rho, then redo NSE to get abar and finally // see if we get back our new p amrex::Real p_new = eos_state.p * 1.05; @@ -412,7 +444,8 @@ void nse_cell_c() eos_state.aux[iabar] = nse_state.abar; eos(eos_input_rt, eos_state); - std::cout << "EOS rho from p consistency check (old method): " << eos_state.p << " " << p_new << std::endl; + std::cout << " old method: invert EOS without considering NSE changes:" << std::endl; + std::cout << " change in p: " << eos_state.p << " " << p_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -429,17 +462,87 @@ void nse_cell_c() nse_rho_abar_from_p(eos_state.T, eos_state.p, eos_state.aux[iye], eos_state.rho, eos_state.aux[iabar]); - std::cout << "updated T: " << eos_state.T << std::endl; - std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; + std::cout << std::endl; + std::cout << " NSE method: use the nse_rho_abar_from_p() function:" << std::endl; + + std::cout << " updated T: " << eos_state.T << std::endl; + std::cout << " change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; // now check if we get back the correct p! eos(eos_input_rt, eos_state); - std::cout << "EOS rho from p consistency check (new method): " << eos_state.p << " " << p_new << std::endl; + std::cout << " change in p: " << eos_state.p << " " << p_new << std::endl; std::cout << std::endl; } + { + + // now explore a troublesome state where tricubic + // interpolation can lead to changing the sign of the neutrino + // energy. + + const amrex::Real T = 3667943011.533362; + const amrex::Real rho = 9434630294.888599; + const amrex::Real Ye = 0.433159096449449; + + amrex::Real logrho = std::log10(rho); + amrex::Real logT = std::log10(T); + + // shift the indices such that it is the starting value + // for the 4x4x4 cube we use for tricubic + int ir0 = nse_get_logrho_index(logrho) - 1; + ir0 = amrex::Clamp(ir0, 1, nse_table_size::nden-3); + + int it0 = nse_get_logT_index(logT) - 1; + it0 = amrex::Clamp(it0, 1, nse_table_size::ntemp-3); + + int ic0 = nse_get_ye_index(Ye) - 1; + ic0 = amrex::Clamp(ic0, 1, nse_table_size::nye-3); + + std::cout << std::endl; + std::cout << "explore a troublesome state" << std::endl << std::endl; + + std::cout << std::format(" rho = {:20.15e}, T = {:20.15e}, Ye = {}\n", + rho, T, Ye); + std::cout << std::format(" indices of lowest corner of 4x4x4 cube: ir0 = {}, it0 = {}, ic0 = {}\n", + ir0, it0, ic0); + + nse_table_t nse_state; + nse_state.T = T; + nse_state.rho = rho; + nse_state.Ye = Ye; + + nse_interp(nse_state); + + std::cout << std::endl; + std::cout << " cubic interpolation gives: e_nu = " << nse_state.e_nu << std::endl; + + amrex::Real enu_min = std::numeric_limits::max(); + amrex::Real enu_max = std::numeric_limits::lowest(); + + for (int ir = ir0; ir < ir0 + 4; ++ir) { + for (int jt = it0; jt < it0 + 4; ++jt) { + for (int ky = ic0; ky < ic0 + 4; ++ky) { + + auto _enu = nse_table::enutab(nse_idx(ir, jt, ky)); + + enu_min = std::min(enu_min, _enu); + enu_max = std::max(enu_max, _enu); + + } + } + } + + std::cout << " bounds of e_nu in the 64 cells used for tricubic: " << enu_min << " " << enu_max << std::endl; + + + std::cout << std::endl; + + } + + + } #endif diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out deleted file mode 100644 index 68f9dbda7f..0000000000 --- a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out +++ /dev/null @@ -1,92 +0,0 @@ -Initializing AMReX (23.12-21-gef38229189e3)... -AMReX (23.12-21-gef38229189e3) initialized -starting the single zone burn... -reading the NSE table (C++) ... -rho, T, Ye = 1230000000 5180000000 0.472 -density value brackets: 9.05 < 9.089905111 < 9.1 -temperature value brackets: 9.71 < 9.71432976 < 9.72 -Ye value brackets: 0.4725 < 0.472 < 0.47 - -4 rho values: -9 55.5223665 -9.05 55.57080116 -9.1 55.61707675 -9.15 55.66151541 - -cubic interpolated value: 55.6078922 - -4 T values: -9.7 55.5223665 -9.71 55.39489494 -9.72 55.22387568 -9.73 54.99343581 - -cubic interpolated value: 55.32712406 - -4 Ye values: -0.475 55.5223665 -0.4725 55.65227675 -0.47 55.78142906 -0.4675 55.90244371 - -cubic interpolated value: 55.678404 - -tricubic interpolated values: -abar = 55.60652462 -bea = 8.755119523 -dyedt = -0.004527880871 -dbeadt = 0.000111945508 -e_nu = 1.101641642e+16 -X(H1) = 0 -X(He3) = 2.723457525e-14 -X(He4) = 0.0003164137294 -X(C12) = 7.975878903e-10 -X(N14) = 3.102372259e-14 -X(O16) = 1.998161464e-09 -X(Ne20) = 2.444604777e-11 -X(Mg24) = 1.029036116e-08 -X(Si28) = 8.173012779e-06 -X(S32) = 1.194076593e-05 -X(Ar36) = 1.243414613e-05 -X(Ca40) = 1.259531492e-05 -X(Ti44) = 0.0001438001413 -X(Cr48) = 0.01146726617 -X(Fe52) = 0.07746625779 -X(Fe54) = 0.9098424639 -X(Ni56) = 0.0007087774945 -X(n) = 2.078852742e-08 -X(p) = 9.843630258e-06 - -testing temperature derivatives of cubic -first finite-difference derivatives -dAbar/dT = -1.072562355e-09 -dbea/dT = -6.867059575e-12 -now using derivative of the interpolant -dAbar/dT = -1.072562604e-09 -dbea/dT = -6.867048589e-12 - - -testing density derivatives of cubic -first finite-difference derivatives -dAbar/drho = 3.987525411e-10 -dbea/drho = 7.619559936e-13 -now using derivative of the interpolant -dAbar/drho = 3.987522836e-10 -dbea/drho = 7.618831514e-13 - -EOS T from e consistency check (old method): 1.395278886e+18 1.38844906e+18 -updated T: 6394534499 -change in abar: 55.60652462 50.26831386 -EOS T from e consistency check (new method): 1.38844906e+18 1.38844906e+18 - -EOS T from p consistency check (old method): 6.622159603e+26 6.577850616e+26 -updated T: 6466757500 -change in abar: 55.60652462 49.50320619 -EOS T from p consistency check (new method): 6.577850616e+26 6.577850616e+26 - -EOS rho from p consistency check (old method): 6.577758474e+26 6.577850616e+26 -updated T: 5180000000 -change in abar: 55.60652462 55.62494615 -EOS rho from p consistency check (new method): 6.577850616e+26 6.577850616e+26 - -AMReX (23.12-21-gef38229189e3) finalized