From d9d42bdd6022d6effb4711b9cb92c21def9d3893 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 12 Sep 2026 08:00:43 -0400 Subject: [PATCH 01/11] rename test_nse_interp -> nse_interp_check this is consistent with our naming convention that test_* tests use a grid of cells. --- .github/workflows/test_nse_interp.yml | 14 +++++++------- .../GNUmakefile | 0 .../Make.package | 0 .../README.md | 0 .../_parameters | 0 .../ci-benchmarks/aprox19.out | 0 .../{test_nse_interp => nse_interp_check}/main.cpp | 0 .../nse_cell.H | 0 8 files changed, 7 insertions(+), 7 deletions(-) rename unit_test/{test_nse_interp => nse_interp_check}/GNUmakefile (100%) rename unit_test/{test_nse_interp => nse_interp_check}/Make.package (100%) rename unit_test/{test_nse_interp => nse_interp_check}/README.md (100%) rename unit_test/{test_nse_interp => nse_interp_check}/_parameters (100%) rename unit_test/{test_nse_interp => nse_interp_check}/ci-benchmarks/aprox19.out (100%) rename unit_test/{test_nse_interp => nse_interp_check}/main.cpp (100%) rename unit_test/{test_nse_interp => nse_interp_check}/nse_cell.H (100%) diff --git a/.github/workflows/test_nse_interp.yml b/.github/workflows/test_nse_interp.yml index d372a0381e..5944361307 100644 --- a/.github/workflows/test_nse_interp.yml +++ b/.github/workflows/test_nse_interp.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/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 100% rename from unit_test/test_nse_interp/README.md rename to unit_test/nse_interp_check/README.md 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/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out similarity index 100% rename from unit_test/test_nse_interp/ci-benchmarks/aprox19.out rename to unit_test/nse_interp_check/ci-benchmarks/aprox19.out 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 100% rename from unit_test/test_nse_interp/nse_cell.H rename to unit_test/nse_interp_check/nse_cell.H From d538a0efa1669231bf578b932fe13ece73d4d5dc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 12 Sep 2026 09:10:46 -0400 Subject: [PATCH 02/11] updates to style --- Docs/source/unit_tests.rst | 6 +++--- unit_test/nse_interp_check/README.md | 7 +++++-- unit_test/nse_interp_check/nse_cell.H | 25 ++++++++++++++++--------- 3 files changed, 24 insertions(+), 14 deletions(-) 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/nse_interp_check/README.md b/unit_test/nse_interp_check/README.md index 383037be55..f11991c5fa 100644 --- a/unit_test/nse_interp_check/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/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H index c1eb3cf3f1..60ff2165c6 100644 --- a/unit_test/nse_interp_check/nse_cell.H +++ b/unit_test/nse_interp_check/nse_cell.H @@ -15,10 +15,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 @@ -46,6 +48,8 @@ void nse_cell_c() // now try 1-d interpolation in T -- for cubic interpolation, we // need to offset the indices one to the left + std::cout << "testing 1-d cubic interpolation at input state" << std::endl; + ir0--; ir0 = amrex::max(1, amrex::min(nse_table_size::nden-3, ir0)); @@ -125,7 +129,9 @@ void nse_cell_c() } - std::cout << "tricubic interpolated values: " << std::endl; + std::cout << std::endl; + + std::cout << "tricubic interpolated values at input state: " << std::endl; nse_table_t nse_state; nse_state.T = unit_test_rp::temperature; @@ -146,6 +152,7 @@ void nse_cell_c() // temp derivatives + std::cout << std::endl; std::cout << std::endl; std::cout << "testing temperature derivatives of cubic" << std::endl; @@ -184,7 +191,6 @@ void nse_cell_c() // dens derivatives - std::cout << std::endl; std::cout << "testing density derivatives of cubic" << std::endl; std::cout << "first finite-difference derivatives" << std::endl; @@ -215,6 +221,7 @@ void nse_cell_c() std::cout << "dAbar/drho = " << dabardrho << std::endl; std::cout << "dbea/drho = " << dbeadrho << std::endl; + std::cout << std::endl; std::cout << std::endl; // @@ -278,7 +285,7 @@ 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 << "EOS T from e consistency check (old method); change in e: " << eos_state.e << " " << e_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -302,7 +309,7 @@ void nse_cell_c() 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 << "EOS T from e consistency check (new method); change in e: " << eos_state.e << " " << e_new << std::endl; std::cout << std::endl; } @@ -344,7 +351,7 @@ 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 << "EOS T from p consistency check (old method); change in p: " << eos_state.p << " " << p_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -368,7 +375,7 @@ void nse_cell_c() 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 << "EOS T from p consistency check (new method); change in p: " << eos_state.p << " " << p_new << std::endl; std::cout << std::endl; } @@ -412,7 +419,7 @@ 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 << "EOS rho from p consistency check (old method); change in p: " << eos_state.p << " " << p_new << std::endl; // attempt 2: // now we try the new interface. This effectively does: @@ -436,7 +443,7 @@ void nse_cell_c() 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 << "EOS rho from p consistency check (new method); change in p: " << eos_state.p << " " << p_new << std::endl; std::cout << std::endl; } From 57c3e0c271a07906fb170efed8b6b7e4d27f3483 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 12 Sep 2026 10:32:25 -0400 Subject: [PATCH 03/11] more formatting --- .../ci-benchmarks/aprox19.out | 196 +++++++++------- unit_test/nse_interp_check/nse_cell.H | 216 +++++++++++++----- 2 files changed, 272 insertions(+), 140 deletions(-) diff --git a/unit_test/nse_interp_check/ci-benchmarks/aprox19.out b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out index 68f9dbda7f..53a324c615 100644 --- a/unit_test/nse_interp_check/ci-benchmarks/aprox19.out +++ b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out @@ -1,92 +1,128 @@ -Initializing AMReX (23.12-21-gef38229189e3)... -AMReX (23.12-21-gef38229189e3) initialized +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 -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 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 + 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 + + 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/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H index 60ff2165c6..3097750bdc 100644 --- a/unit_test/nse_interp_check/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 @@ -45,10 +48,10 @@ 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" << std::endl; + 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)); @@ -59,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)}; + // density interpolation - amrex::Real yes[] = {nse_table_ye(ic0), - nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + 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)), @@ -85,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)), @@ -102,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)), @@ -120,18 +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 << std::endl; - std::cout << "tricubic interpolated values at input state: " << std::endl; + std::cout << "tricubic interpolated values at input state\n\n"; nse_table_t nse_state; nse_state.T = unit_test_rp::temperature; @@ -140,13 +150,14 @@ 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; } @@ -154,9 +165,9 @@ void nse_cell_c() std::cout << std::endl; std::cout << std::endl; - std::cout << "testing temperature derivatives of cubic" << 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; @@ -173,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); @@ -184,16 +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 << "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; @@ -207,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); @@ -218,8 +231,8 @@ 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; @@ -250,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; @@ -285,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); change in e: " << 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: @@ -302,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); change in e: " << eos_state.e << " " << e_new << std::endl; + std::cout << " change in e: " << eos_state.e << " " << e_new << std::endl; std::cout << std::endl; } @@ -318,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; @@ -351,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); change in p: " << 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: @@ -368,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); change in p: " << eos_state.p << " " << p_new << std::endl; + std::cout << " change in p: " << eos_state.p << " " << p_new << std::endl; std::cout << std::endl; } @@ -386,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; @@ -405,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; @@ -419,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); change in p: " << 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: @@ -436,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); change in p: " << 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 From 11940e85acc980ef3c123cb217193daf5805aee5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 12 Sep 2026 12:43:59 -0400 Subject: [PATCH 04/11] implement monotonic Hermite polynomial cubic interp for NSE --- nse_tabular/nse_table.H | 115 ++++++++++++++++++++++++++++++++-------- 1 file changed, 93 insertions(+), 22 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 26aeecf3e7..16027442ec 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -175,26 +175,74 @@ int nse_get_ye_index(const amrex::Real ye) { return ic0 + 1; } + +// Slope Limiter to preserve monotonicity +// See C. Moler, Numerical Computing with Matlab, 2004. +// :doi:`10.1137/1.9780898717952`, Chapter 3.4 +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real +limit_slope(const amrex::Real dm, const amrex::Real dp) { + + // note: we assume that dx is constant in our interval + + if (dm * dp <= 0.0_rt) { + return 0.0_rt; + } + + return 2.0_rt / (1.0_rt / dm + 1.0_rt / dp); +} + + /// /// given 4 points (xs, fs), with spacing dx, return the interplated -/// value of f at point x by fitting a cubic to the points +/// value of f at point x by fitting a monotone cubic Hermite +/// polynomial to the points: +/// +/// H(x) = f0 H_{1,0}(x) + f1 H_{1,1}(x) + f'0 Ĥ_{1,0}(x) + f'1 Ĥ_{1,0}(x) +/// +/// The slope f' is limited to preserve monotonicity. /// AMREX_GPU_HOST_DEVICE AMREX_INLINE -amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, const amrex::Real dx, const amrex::Real x) { +amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, + const amrex::Real dx, const amrex::Real x) { + + // given the 4 points xs, the point x is in [xs[1], xs[2]] + // we'll use i0 = 1 as the reference for the interpolation + + // we have already guaranteed that all of the xs are in the + // domain, so we don't need to worry about boundaries. - // fit a cubic of the form - // f(x) = a (x - x_i)**3 + b (x - x_i)**2 + c (x - x_i) + d - // to the data (xs, fs) - // we take x_i to be x[1] + // Compute the slope at i0 and i0 + 1 + const auto i0 = 1; - amrex::Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * amrex::Math::powi<3>(dx)); - amrex::Real b = (-2 * fs[1] + fs[2] + fs[0]) / (2 * dx * dx); - amrex::Real c = (-3 * fs[1] + 6 * fs[2] - fs[3] - 2 * fs[0]) / (6 * dx); - amrex::Real d = fs[1]; + const auto dx_inv = 1.0_rt / dx; - return a * amrex::Math::powi<3>(x - xs[1]) + - b * amrex::Math::powi<2>(x - xs[1]) + c * (x - xs[1]) + d; + amrex::Real m0, m1; + m1 = m0 = (fs[i0+1] - fs[i0]) * dx_inv; + // Compute the slope from the other side and limit the slope to + // preserve monotonicity + amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + m0 = limit_slope(d_im1, m0); + + amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + m1 = limit_slope(m1, d_ip1); + + // Compute cubic Hermite polynomial basis + const amrex::Real t = (x - xs[i0]) * dx_inv; + const amrex::Real t2 = t * t; + const amrex::Real t3 = t2 * t; + + const amrex::Real H1 = -2.0_rt * t3 + 3.0_rt * t2; + const amrex::Real H0 = -H1 + 1.0_rt; + const amrex::Real Hhat1 = dx * (t3 - t2); + const amrex::Real Hhat0 = dx * (t3 - 2.0_rt * t2 + t); + + // Construct the interpolant + // H_3(x) = f0 H_{1,0}(x) + f1 H_{1,1}(x) + f'0 Ĥ_{1,0}(x) + f'1 Ĥ_{1,0}(x) + + const amrex::Real fv = fs[i0] * H0 + fs[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; + return fv; } /// @@ -205,18 +253,41 @@ amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, const amrex::Rea AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, const amrex::Real dx, const amrex::Real x) { - // fit a cubic of the form - // f(x) = a (x - x_i)**3 + b (x - x_i)**2 + c (x - x_i) + d - // to the data (xs, fs) - // we take x_i to be x[1] - // then return dfdx = 3 a (x - x_i)**2 + 2 b (x - x_i) + c + // given the 4 points xs, the point x is in [xs[1], xs[2]] + // we'll use i0 = 1 as the reference for the interpolation + + // we have already guaranteed that all of the xs are in the + // domain, so we don't need to worry about boundaries. + + // Compute the slope at i0 and i0 + 1 + const auto i0 = 1; + + const auto dx_inv = 1.0_rt / dx; + + amrex::Real m0, m1; + m1 = m0 = (fs[i0+1] - fs[i0]) * dx_inv; + + // Compute the slope from the other side and limit the slope to + // preserve monotonicity + amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + m0 = limit_slope(d_im1, m0); + + amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + m1 = limit_slope(m1, d_ip1); + + // Compute cubic Hermite polynomial basis + const amrex::Real t = (x - xs[i0]) * dx_inv; + const amrex::Real t2 = t * t; + + const amrex::Real dH1dx = -6.0_rt * (t2 - t) * dx_inv; + const amrex::Real dH0dx = -dH1dx; + const amrex::Real dHhat1dx = 3.0_rt * t2 - 2.0_rt * t; + const amrex::Real dHhat0dx = 3.0_rt * t2 - 4.0_rt * t + 1.0_rt; - amrex::Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * amrex::Math::powi<3>(dx)); - amrex::Real b = (-2 * fs[1] + fs[2] + fs[0]) / (2 * dx * dx); - amrex::Real c = (-3 * fs[1] + 6 * fs[2] - fs[3] - 2 * fs[0]) / (6 * dx); - //amrex::Real d = fs[1]; + const amrex::Real dfdx = fs[i0] * dH0dx + fs[i0+1] * dH1dx + + m0 * dHhat0dx + m1 * dHhat1dx; - return 3.0_rt * a * amrex::Math::powi<2>(x - xs[1]) + 2.0_rt * b * (x - xs[1]) + c; + return dfdx; } From db71e0bc2af44fc6ebae5eb54c7024eab3fef057 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 12 Sep 2026 16:15:26 -0400 Subject: [PATCH 05/11] update --- nse_tabular/nse_table.H | 220 ++++++++++++++++---------- unit_test/nse_interp_check/nse_cell.H | 52 +++--- 2 files changed, 156 insertions(+), 116 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 16027442ec..b82099590b 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -204,7 +204,8 @@ limit_slope(const amrex::Real dm, const amrex::Real dp) { /// AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, - const amrex::Real dx, const amrex::Real x) { + const amrex::Real dx, const amrex::Real x, + const bool lo_bnd, const bool hi_bnd) { // given the 4 points xs, the point x is in [xs[1], xs[2]] // we'll use i0 = 1 as the reference for the interpolation @@ -222,11 +223,15 @@ amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, // Compute the slope from the other side and limit the slope to // preserve monotonicity - amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; - m0 = limit_slope(d_im1, m0); + if (! lo_bnd) { + amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + m0 = limit_slope(d_im1, m0); + } - amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; - m1 = limit_slope(m1, d_ip1); + if (! hi_bnd) { + amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + m1 = limit_slope(m1, d_ip1); + } // Compute cubic Hermite polynomial basis const amrex::Real t = (x - xs[i0]) * dx_inv; @@ -251,7 +256,9 @@ amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, /// points and differentiating the interpolant /// AMREX_GPU_HOST_DEVICE AMREX_INLINE -amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, const amrex::Real dx, const amrex::Real x) { +amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, + const amrex::Real dx, const amrex::Real x, + const bool lo_bnd, const bool hi_bnd) { // given the 4 points xs, the point x is in [xs[1], xs[2]] // we'll use i0 = 1 as the reference for the interpolation @@ -269,11 +276,15 @@ amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, const amre // Compute the slope from the other side and limit the slope to // preserve monotonicity - amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; - m0 = limit_slope(d_im1, m0); + if (! lo_bnd) { + amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + m0 = limit_slope(d_im1, m0); + } - amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; - m1 = limit_slope(m1, d_ip1); + if (! hi_bnd) { + amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + m1 = limit_slope(m1, d_ip1); + } // Compute cubic Hermite polynomial basis const amrex::Real t = (x - xs[i0]) * dx_inv; @@ -338,37 +349,48 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real tricubic(const int ir0, const int it0, const int ic0, const amrex::Real rho, const amrex::Real temp, const amrex::Real ye, const T& data) { - const amrex::Real yes[] = {nse_table_ye(ic0), + // as we come in, we expect temp to be in it0:it0+1, and we also + // are assured that it0 and it0+1 don't go out of bounds. Likewise + // for rho and ye. Our 1D cubic interpolation routine will handle + // the boundary conditions for us. + + const amrex::Real yes[] = {nse_table_ye(std::max(1, ic0-1)), + nse_table_ye(ic0), nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + nse_table_ye(std::min(nse_table_size::nye, ic0+2))}; - const amrex::Real Ts[] = {nse_table_logT(it0), + const amrex::Real Ts[] = {nse_table_logT(std::max(1, it0-1)), + nse_table_logT(it0), nse_table_logT(it0+1), - nse_table_logT(it0+2), - nse_table_logT(it0+3)}; + nse_table_logT(std::min(nse_table_size::ntemp, it0+2))}; - const amrex::Real rhos[] = {nse_table_logrho(ir0), + const amrex::Real rhos[] = {nse_table_logrho(std::max(1, ir0-1)), + nse_table_logrho(ir0), nse_table_logrho(ir0+1), - nse_table_logrho(ir0+2), - nse_table_logrho(ir0+3)}; + nse_table_logrho(std::min(nse_table_size::nden, ir0+2))}; // first do the 16 ye interpolations // the first index will be rho and the second will be T amrex::Real d1[4][4]; + bool lo_bnd = (ic0 == 1); + bool hi_bnd = (ic0 + 1 == nse_table_size::nye); + for (int ii = 0; ii < 4; ++ii) { + int idx_r = std::max(1, std::min(ir0 - 1 + ii, nse_table_size::nden)); + for (int jj = 0; jj < 4; ++jj) { + int idx_t = std::max(1, std::min(it0 - 1 + jj, nse_table_size::ntemp)); - const amrex::Real _d[] = {data(nse_idx(ir0+ii, it0+jj, ic0)), - data(nse_idx(ir0+ii, it0+jj, ic0+1)), - data(nse_idx(ir0+ii, it0+jj, ic0+2)), - data(nse_idx(ir0+ii, it0+jj, ic0+3))}; + const amrex::Real _d[] = {data(nse_idx(idx_r, idx_t, std::max(1, ic0-1))), + data(nse_idx(idx_r, idx_t, ic0)), + data(nse_idx(idx_r, idx_t, ic0+1)), + data(nse_idx(idx_r, idx_t, std::min(ic0+2, nse_table_size::nye)))}; // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye); + d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); } } @@ -376,15 +398,20 @@ amrex::Real tricubic(const int ir0, const int it0, const int ic0, amrex::Real d2[4]; - for (int ii = 0; ii < 4; ++ii) { + lo_bnd = (it0 == 1); + hi_bnd = (it0 + 1 == nse_table_size::ntemp); + for (int ii = 0; ii < 4; ++ii) { const amrex::Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]}; - d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp); + d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); } // finally do the remaining interpolation over rho - amrex::Real val = cubic(rhos, d2, nse_table_size::dlogrho, rho); + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); + + amrex::Real val = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); return val; @@ -399,37 +426,48 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, const amrex::Real rho, const amrex::Real temp, const amrex::Real ye, const T& data) { - const amrex::Real yes[] = {nse_table_ye(ic0), + // as we come in, we expect temp to be in it0:it0+1, and we also + // are assured that it0 and it0+1 don't go out of bounds. Likewise + // for rho and ye. Our 1D cubic interpolation routine will handle + // the boundary conditions for us. + + const amrex::Real yes[] = {nse_table_ye(std::max(1, ic0-1)), + nse_table_ye(ic0), nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + nse_table_ye(std::min(nse_table_size::nye, ic0+2))}; - const amrex::Real Ts[] = {nse_table_logT(it0), + const amrex::Real Ts[] = {nse_table_logT(std::max(1, it0-1)), + nse_table_logT(it0), nse_table_logT(it0+1), - nse_table_logT(it0+2), - nse_table_logT(it0+3)}; + nse_table_logT(std::min(nse_table_size::ntemp, it0+2))}; - const amrex::Real rhos[] = {nse_table_logrho(ir0), + const amrex::Real rhos[] = {nse_table_logrho(std::max(1, ir0-1)), + nse_table_logrho(ir0), nse_table_logrho(ir0+1), - nse_table_logrho(ir0+2), - nse_table_logrho(ir0+3)}; + nse_table_logrho(std::min(nse_table_size::nden, ir0+2))}; // first do the 16 ye interpolations // the first index will be rho and the second will be T amrex::Real d1[4][4]; + bool lo_bnd = (ic0 == 1); + bool hi_bnd = (ic0 + 1 == nse_table_size::nye); + for (int ii = 0; ii < 4; ++ii) { + int idx_r = std::max(1, std::min(ir0 - 1 + ii, nse_table_size::nden)); + for (int jj = 0; jj < 4; ++jj) { + int idx_t = std::max(1, std::min(it0 - 1 + jj, nse_table_size::ntemp)); - const amrex::Real _d[] = {data(nse_idx(ir0+ii, it0+jj, ic0)), - data(nse_idx(ir0+ii, it0+jj, ic0+1)), - data(nse_idx(ir0+ii, it0+jj, ic0+2)), - data(nse_idx(ir0+ii, it0+jj, ic0+3))}; + const amrex::Real _d[] = {data(nse_idx(idx_r, idx_t, std::max(1, ic0-1))), + data(nse_idx(idx_r, idx_t, ic0)), + data(nse_idx(idx_r, idx_t, ic0+1)), + data(nse_idx(idx_r, idx_t, std::min(ic0+2, nse_table_size::nye)))}; // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye); + d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); } } @@ -437,16 +475,21 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, amrex::Real d2[4]; - for (int jj = 0; jj < 4; ++jj) { + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); + for (int jj = 0; jj < 4; ++jj) { const amrex::Real _d[] = {d1[0][jj], d1[1][jj], d1[2][jj], d1[3][jj]}; - d2[jj] = cubic(rhos, _d, nse_table_size::dlogrho, rho); + d2[jj] = cubic(rhos, _d, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); } // finally do the remaining interpolation over T, but return // the derivative of the interpolant - amrex::Real val = cubic_deriv(Ts, d2, nse_table_size::dlogT, temp); + lo_bnd = (it0 == 1); + hi_bnd = (it0 + 1 == nse_table_size::ntemp); + + amrex::Real val = cubic_deriv(Ts, d2, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); return val; @@ -462,37 +505,48 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real tricubic_drho(const int ir0, const int it0, const int ic0, const amrex::Real rho, const amrex::Real temp, const amrex::Real ye, const T& data) { - const amrex::Real yes[] = {nse_table_ye(ic0), + // as we come in, we expect temp to be in it0:it0+1, and we also + // are assured that it0 and it0+1 don't go out of bounds. Likewise + // for rho and ye. Our 1D cubic interpolation routine will handle + // the boundary conditions for us. + + const amrex::Real yes[] = {nse_table_ye(std::max(1, ic0-1)), + nse_table_ye(ic0), nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + nse_table_ye(std::min(nse_table_size::nye, ic0+2))}; - const amrex::Real Ts[] = {nse_table_logT(it0), + const amrex::Real Ts[] = {nse_table_logT(std::max(1, it0-1)), + nse_table_logT(it0), nse_table_logT(it0+1), - nse_table_logT(it0+2), - nse_table_logT(it0+3)}; + nse_table_logT(std::min(nse_table_size::ntemp, it0+2))}; - const amrex::Real rhos[] = {nse_table_logrho(ir0), + const amrex::Real rhos[] = {nse_table_logrho(std::max(1, ir0-1)), + nse_table_logrho(ir0), nse_table_logrho(ir0+1), - nse_table_logrho(ir0+2), - nse_table_logrho(ir0+3)}; + nse_table_logrho(std::min(nse_table_size::nden, ir0+2))}; // first do the 16 ye interpolations // the first index will be rho and the second will be T amrex::Real d1[4][4]; + bool lo_bnd = (ic0 == 1); + bool hi_bnd = (ic0 + 1 == nse_table_size::nye); + for (int ii = 0; ii < 4; ++ii) { + int idx_r = std::max(1, std::min(ir0 - 1 + ii, nse_table_size::nden)); + for (int jj = 0; jj < 4; ++jj) { + int idx_t = std::max(1, std::min(it0 - 1 + jj, nse_table_size::ntemp)); - const amrex::Real _d[] = {data(nse_idx(ir0+ii, it0+jj, ic0)), - data(nse_idx(ir0+ii, it0+jj, ic0+1)), - data(nse_idx(ir0+ii, it0+jj, ic0+2)), - data(nse_idx(ir0+ii, it0+jj, ic0+3))}; + const amrex::Real _d[] = {data(nse_idx(idx_r, idx_t, std::max(1, ic0-1))), + data(nse_idx(idx_r, idx_t, ic0)), + data(nse_idx(idx_r, idx_t, ic0+1)), + data(nse_idx(idx_r, idx_t, std::min(ic0+2, nse_table_size::nye)))}; // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye); + d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); } } @@ -500,16 +554,23 @@ amrex::Real tricubic_drho(const int ir0, const int it0, const int ic0, amrex::Real d2[4]; + lo_bnd = (it0 == 1); + hi_bnd = (it0 + 1 == nse_table_size::ntemp); + for (int ii = 0; ii < 4; ++ii) { const amrex::Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]}; - d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp); + d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); } // finally do the remaining interpolation over rho, but return // the derivative of the interpolant - amrex::Real val = cubic_deriv(rhos, d2, nse_table_size::dlogrho, rho); + + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); + + amrex::Real val = cubic_deriv(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); return val; @@ -573,18 +634,15 @@ void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) { } else { // for a cubic interpolant, we need 4 points that span the data value - // for temperature, these will be it0, it0+1, it0+2, it0+3 - // with the idea that the temperature we want is between it0+1 and it0+2 - // so we offset one to the left and also ensure that we don't go off the table - - int ir0 = nse_get_logrho_index(rholog) - 1; - ir0 = amrex::Clamp(ir0, 1, nse_table_size::nden-3); + // for temperature, these will be it0-1, it0, it0+1, it0+2 + // with the idea that the temperature we want is between it0 and it0+1 - int it0 = nse_get_logT_index(tlog) - 1; - it0 = amrex::Clamp(it0, 1, nse_table_size::ntemp-3); + // we will handle the case where it0 or it0+1 is at the bounds of + // the table in the caller. - int ic0 = nse_get_ye_index(yet) - 1; - ic0 = amrex::Clamp(ic0, 1, nse_table_size::nye-3); + int ir0 = nse_get_logrho_index(rholog); + int it0 = nse_get_logT_index(tlog); + int ic0 = nse_get_ye_index(yet); nse_state.abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab); nse_state.bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab); @@ -640,14 +698,9 @@ nse_interp_dT(const amrex::Real temp, const amrex::Real rho, const amrex::Real y yet = amrex::Clamp(yet, yemin, yemax); } - int ir0 = nse_get_logrho_index(rholog) - 1; - ir0 = amrex::Clamp(ir0, 1, nse_table_size::nden-3); - - int it0 = nse_get_logT_index(tlog) - 1; - it0 = amrex::Clamp(it0, 1, nse_table_size::ntemp-3); - - int ic0 = nse_get_ye_index(yet) - 1; - ic0 = amrex::Clamp(ic0, 1, nse_table_size::nye-3); + int ir0 = nse_get_logrho_index(rholog); + int it0 = nse_get_logT_index(tlog); + int ic0 = nse_get_ye_index(yet); // note: this is returning the derivative wrt log10(T), so we need to // convert to d/dT @@ -693,14 +746,9 @@ nse_interp_drho(const amrex::Real temp, const amrex::Real rho, const amrex::Real yet = amrex::Clamp(yet, yemin, yemax); } - int ir0 = nse_get_logrho_index(rholog) - 1; - ir0 = amrex::Clamp(ir0, 1, nse_table_size::nden-3); - - int it0 = nse_get_logT_index(tlog) - 1; - it0 = amrex::Clamp(it0, 1, nse_table_size::ntemp-3); - - int ic0 = nse_get_ye_index(yet) - 1; - ic0 = amrex::Clamp(ic0, 1, nse_table_size::nye-3); + int ir0 = nse_get_logrho_index(rholog); + int it0 = nse_get_logT_index(tlog); + int ic0 = nse_get_ye_index(yet); // note: this is returning the derivative wrt log10(rho), so we need to // convert to d/drho diff --git a/unit_test/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H index 3097750bdc..5e084aa864 100644 --- a/unit_test/nse_interp_check/nse_cell.H +++ b/unit_test/nse_interp_check/nse_cell.H @@ -53,30 +53,22 @@ void nse_cell_c() 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)); - - it0--; - it0 = amrex::max(1, amrex::min(nse_table_size::ntemp-3, it0)); - - ic0--; - ic0 = amrex::max(1, amrex::min(nse_table_size::nye-3, ic0)); - { // density interpolation + // note: assuming we are not at a boundary - amrex::Real rhos[] = {nse_table_logrho(ir0), + amrex::Real rhos[] = {nse_table_logrho(ir0-1), + nse_table_logrho(ir0), nse_table_logrho(ir0+1), - nse_table_logrho(ir0+2), - nse_table_logrho(ir0+3)}; + nse_table_logrho(ir0+2)}; std::cout << " 4 rho values: " << std::endl; - amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), + amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0-1, it0, ic0)), + nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0+1, it0, ic0)), - nse_table::abartab(nse_idx(ir0+2, it0, ic0)), - nse_table::abartab(nse_idx(ir0+3, it0, ic0))}; + nse_table::abartab(nse_idx(ir0+2, it0, ic0))}; for (int ii = 0; ii < 4; ++ii) { std::cout << std::format(" {:5.3f} {:15.10g}\n", rhos[ii], _d[ii]); @@ -84,24 +76,24 @@ void nse_cell_c() std::cout << std::endl; std::cout << " cubic interpolated value: " << - cubic(rhos, _d, nse_table_size::dlogrho, logrho) << std::endl << std::endl; + cubic(rhos, _d, nse_table_size::dlogrho, logrho, false, false) << std::endl << std::endl; } { // temperature interpolation - amrex::Real Ts[] = {nse_table_logT(it0), + amrex::Real Ts[] = {nse_table_logT(it0-1), + nse_table_logT(it0), nse_table_logT(it0+1), - nse_table_logT(it0+2), - nse_table_logT(it0+3)}; + nse_table_logT(it0+2)}; std::cout << " 4 T values: " << std::endl; - amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), + amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0-1, ic0)), + nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0, it0+1, ic0)), - nse_table::abartab(nse_idx(ir0, it0+2, ic0)), - nse_table::abartab(nse_idx(ir0, it0+3, ic0))}; + nse_table::abartab(nse_idx(ir0, it0+2, ic0))}; for (int ii = 0; ii < 4; ++ii) { std::cout << std::format(" {:5.3f} {:15.10g}\n", Ts[ii], _d[ii]); @@ -109,7 +101,7 @@ void nse_cell_c() std::cout << std::endl; std::cout << " cubic interpolated value: " << - cubic(Ts, _d, nse_table_size::dlogT, logT) << std::endl << std::endl; + cubic(Ts, _d, nse_table_size::dlogT, logT, false, false) << std::endl << std::endl; } @@ -117,17 +109,17 @@ void nse_cell_c() // Ye interpolation - amrex::Real yes[] = {nse_table_ye(ic0), + amrex::Real yes[] = {nse_table_ye(ic0-1), + nse_table_ye(ic0), nse_table_ye(ic0+1), - nse_table_ye(ic0+2), - nse_table_ye(ic0+3)}; + nse_table_ye(ic0+2)}; std::cout << " 4 Ye values: " << std::endl; - amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0)), + amrex::Real _d[] = {nse_table::abartab(nse_idx(ir0, it0, ic0-1)), + nse_table::abartab(nse_idx(ir0, it0, ic0)), nse_table::abartab(nse_idx(ir0, it0, ic0+1)), - nse_table::abartab(nse_idx(ir0, it0, ic0+2)), - nse_table::abartab(nse_idx(ir0, it0, ic0+3))}; + nse_table::abartab(nse_idx(ir0, it0, ic0+2))}; for (int ii = 0; ii < 4; ++ii) { std::cout << std::format(" {:5.3f} {:15.10g}\n", yes[ii], _d[ii]); @@ -135,7 +127,7 @@ void nse_cell_c() std::cout << std::endl; std::cout << " cubic interpolated value: " << - cubic(yes, _d, -nse_table_size::dye, unit_test_rp::ye) << std::endl << std::endl; + cubic(yes, _d, -nse_table_size::dye, unit_test_rp::ye, false, false) << std::endl << std::endl; } From f11124a1b5975d8cb66af2578f23429d8d84f540 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 13 Sep 2026 09:23:38 -0400 Subject: [PATCH 06/11] more interface work --- nse_tabular/nse_table.H | 84 +++++++++------------------ unit_test/nse_interp_check/nse_cell.H | 14 +++-- 2 files changed, 34 insertions(+), 64 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index b82099590b..fc5a3f28c1 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -203,9 +204,10 @@ limit_slope(const amrex::Real dm, const amrex::Real dp) { /// The slope f' is limited to preserve monotonicity. /// AMREX_GPU_HOST_DEVICE AMREX_INLINE -amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, - const amrex::Real dx, const amrex::Real x, - const bool lo_bnd, const bool hi_bnd) { +std::pair +cubic(const amrex::Real* xs, const amrex::Real* fs, + const amrex::Real dx, const amrex::Real x, + const bool lo_bnd, const bool hi_bnd) { // given the 4 points xs, the point x is in [xs[1], xs[2]] // we'll use i0 = 1 as the reference for the interpolation @@ -247,49 +249,8 @@ amrex::Real cubic(const amrex::Real* xs, const amrex::Real* fs, // H_3(x) = f0 H_{1,0}(x) + f1 H_{1,1}(x) + f'0 Ĥ_{1,0}(x) + f'1 Ĥ_{1,0}(x) const amrex::Real fv = fs[i0] * H0 + fs[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; - return fv; -} - -/// -/// given 4 points (xs, fs), with spacing dx between the xs, return -/// the derivative of f at point x by fitting a cubic to the -/// points and differentiating the interpolant -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, - const amrex::Real dx, const amrex::Real x, - const bool lo_bnd, const bool hi_bnd) { - - // given the 4 points xs, the point x is in [xs[1], xs[2]] - // we'll use i0 = 1 as the reference for the interpolation - - // we have already guaranteed that all of the xs are in the - // domain, so we don't need to worry about boundaries. - - // Compute the slope at i0 and i0 + 1 - const auto i0 = 1; - - const auto dx_inv = 1.0_rt / dx; - - amrex::Real m0, m1; - m1 = m0 = (fs[i0+1] - fs[i0]) * dx_inv; - - // Compute the slope from the other side and limit the slope to - // preserve monotonicity - if (! lo_bnd) { - amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; - m0 = limit_slope(d_im1, m0); - } - - if (! hi_bnd) { - amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; - m1 = limit_slope(m1, d_ip1); - } - - // Compute cubic Hermite polynomial basis - const amrex::Real t = (x - xs[i0]) * dx_inv; - const amrex::Real t2 = t * t; + // now the derivative const amrex::Real dH1dx = -6.0_rt * (t2 - t) * dx_inv; const amrex::Real dH0dx = -dH1dx; const amrex::Real dHhat1dx = 3.0_rt * t2 - 2.0_rt * t; @@ -298,7 +259,7 @@ amrex::Real cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, const amrex::Real dfdx = fs[i0] * dH0dx + fs[i0+1] * dH1dx + m0 * dHhat0dx + m1 * dHhat1dx; - return dfdx; + return {fv, dfdx}; } @@ -390,7 +351,8 @@ amrex::Real tricubic(const int ir0, const int it0, const int ic0, // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; } } @@ -403,7 +365,8 @@ amrex::Real tricubic(const int ir0, const int it0, const int ic0, for (int ii = 0; ii < 4; ++ii) { const amrex::Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]}; - d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + d2[ii] = fv; } // finally do the remaining interpolation over rho @@ -411,9 +374,9 @@ amrex::Real tricubic(const int ir0, const int it0, const int ic0, lo_bnd = (ir0 == 1); hi_bnd = (ir0 + 1 == nse_table_size::nden); - amrex::Real val = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); - return val; + return fv; } @@ -467,7 +430,8 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; } } @@ -480,7 +444,8 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, for (int jj = 0; jj < 4; ++jj) { const amrex::Real _d[] = {d1[0][jj], d1[1][jj], d1[2][jj], d1[3][jj]}; - d2[jj] = cubic(rhos, _d, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(rhos, _d, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + d2[jj] = fv; } // finally do the remaining interpolation over T, but return @@ -489,9 +454,9 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, lo_bnd = (it0 == 1); hi_bnd = (it0 + 1 == nse_table_size::ntemp); - amrex::Real val = cubic_deriv(Ts, d2, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(Ts, d2, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); - return val; + return dfdx; } @@ -546,7 +511,9 @@ amrex::Real tricubic_drho(const int ir0, const int it0, const int ic0, // note that the ye values are monotonically decreasing, // so the "dx" needs to be negative - d1[ii][jj] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; + } } @@ -560,7 +527,8 @@ amrex::Real tricubic_drho(const int ir0, const int it0, const int ic0, for (int ii = 0; ii < 4; ++ii) { const amrex::Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]}; - d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + d2[ii] = fv; } // finally do the remaining interpolation over rho, but return @@ -570,9 +538,9 @@ amrex::Real tricubic_drho(const int ir0, const int it0, const int ic0, lo_bnd = (ir0 == 1); hi_bnd = (ir0 + 1 == nse_table_size::nden); - amrex::Real val = cubic_deriv(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); - return val; + return dfdx; } diff --git a/unit_test/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H index 5e084aa864..2b4ea4ce03 100644 --- a/unit_test/nse_interp_check/nse_cell.H +++ b/unit_test/nse_interp_check/nse_cell.H @@ -75,8 +75,9 @@ void nse_cell_c() } std::cout << std::endl; - std::cout << " cubic interpolated value: " << - cubic(rhos, _d, nse_table_size::dlogrho, logrho, false, false) << std::endl << std::endl; + const auto [fv, dfdx] = cubic(rhos, _d, nse_table_size::dlogrho, logrho, false, false); + + std::cout << " cubic interpolated value: " << fv << std::endl << std::endl; } { @@ -100,8 +101,8 @@ void nse_cell_c() } std::cout << std::endl; - std::cout << " cubic interpolated value: " << - cubic(Ts, _d, nse_table_size::dlogT, logT, false, false) << std::endl << std::endl; + const auto [fv, dfdx] = cubic(Ts, _d, nse_table_size::dlogT, logT, false, false); + std::cout << " cubic interpolated value: " << fv << std::endl << std::endl; } @@ -126,8 +127,9 @@ void nse_cell_c() } std::cout << std::endl; - std::cout << " cubic interpolated value: " << - cubic(yes, _d, -nse_table_size::dye, unit_test_rp::ye, false, false) << std::endl << std::endl; + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, unit_test_rp::ye, false, false); + + std::cout << " cubic interpolated value: " << fv << std::endl << std::endl; } From 04a7e5e366976ecaadf7577f4a5cceedf3c7a95f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 13 Sep 2026 09:44:31 -0400 Subject: [PATCH 07/11] more work --- nse_tabular/nse_table.H | 108 ++++++++++++++++++++++++++++++++++------ 1 file changed, 94 insertions(+), 14 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index fc5a3f28c1..4834f425e9 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -264,6 +264,73 @@ cubic(const amrex::Real* xs, const amrex::Real* fs, } +/// +/// an alternate version of the monotone cubic interpolation. +/// Here we use fs to compute the limiter but dfs_dx is what +/// we interpolate. +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +std::pair +cubic_split_limiter(const amrex::Real* xs, + const amrex::Real* fs, const amrex::Real* dfs_dx, + const amrex::Real dx, const amrex::Real x, + const bool lo_bnd, const bool hi_bnd) { + + // given the 4 points xs, the point x is in [xs[1], xs[2]] + // we'll use i0 = 1 as the reference for the interpolation + + // we have already guaranteed that all of the xs are in the + // domain, so we don't need to worry about boundaries. + + // Compute the slope at i0 and i0 + 1 + const auto i0 = 1; + + const auto dx_inv = 1.0_rt / dx; + + amrex::Real m0, m1; + m1 = m0 = (fs[i0+1] - fs[i0]) * dx_inv; + + // Compute the slope from the other side and limit the slope to + // preserve monotonicity + if (! lo_bnd) { + amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + m0 = limit_slope(d_im1, m0); + } + + if (! hi_bnd) { + amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + m1 = limit_slope(m1, d_ip1); + } + + // Compute cubic Hermite polynomial basis + const amrex::Real t = (x - xs[i0]) * dx_inv; + const amrex::Real t2 = t * t; + const amrex::Real t3 = t2 * t; + + const amrex::Real H1 = -2.0_rt * t3 + 3.0_rt * t2; + const amrex::Real H0 = -H1 + 1.0_rt; + const amrex::Real Hhat1 = dx * (t3 - t2); + const amrex::Real Hhat0 = dx * (t3 - 2.0_rt * t2 + t); + + // Construct the interpolant + // H_3(x) = f0 H_{1,0}(x) + f1 H_{1,1}(x) + f'0 Ĥ_{1,0}(x) + f'1 Ĥ_{1,0}(x) + + const amrex::Real fv = dfs_dx[i0] * H0 + dfs_dx[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; + + // now the derivative + const amrex::Real dH1dx = -6.0_rt * (t2 - t) * dx_inv; + const amrex::Real dH0dx = -dH1dx; + const amrex::Real dHhat1dx = 3.0_rt * t2 - 2.0_rt * t; + const amrex::Real dHhat0dx = 3.0_rt * t2 - 4.0_rt * t + 1.0_rt; + + const amrex::Real dfdx = dfs_dx[i0] * dH0dx + dfs_dx[i0+1] * dH1dx + + m0 * dHhat0dx + m1 * dHhat1dx; + + return {fv, dfdx}; + +} + + template AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::Real trilinear(const int ir1, const int it1, const int ic1, @@ -394,6 +461,13 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, // for rho and ye. Our 1D cubic interpolation routine will handle // the boundary conditions for us. + // note, our cubic is nonlinear in the function values (because of + // the limiter), so we need to do the interpolations in the same + // order (Ye -> T -> rho) since they don't commute. For a variable f, + // we want to store the interpolated results of f and df/dT, and then + // for the final rho interpolation, we want to use f to compute the + // limiter for the cubic but interpolate on the values df/dT. + const amrex::Real yes[] = {nse_table_ye(std::max(1, ic0-1)), nse_table_ye(ic0), nse_table_ye(ic0+1), @@ -435,28 +509,34 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, } } - // now do the 4 rho interpolations (one in each T plane) + // now do the 4 T interpolations (one in each rho plane) - amrex::Real d2[4]; + // we need to store both q and dq/dT interpolated in the T + // direction in all 4 rho planes. - lo_bnd = (ir0 == 1); - hi_bnd = (ir0 + 1 == nse_table_size::nden); + amrex::Real d2[4]; + amrex::Real d2_deriv[4]; + lo_bnd = (it0 == 1); + hi_bnd = (it0 + 1 == nse_table_size::ntemp); - for (int jj = 0; jj < 4; ++jj) { - const amrex::Real _d[] = {d1[0][jj], d1[1][jj], d1[2][jj], d1[3][jj]}; - const auto [fv, dfdx] = cubic(rhos, _d, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); - d2[jj] = fv; + for (int ii = 0; ii < 4; ++ii) { + const amrex::Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]}; + const auto [fv, dfdx] = cubic(Ts, _d, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + d2[ii] = fv; + d2_deriv[ii] = dfdx; } - // finally do the remaining interpolation over T, but return - // the derivative of the interpolant + // finally do the remaining interpolation over rho, but we want to limit on + // q but interpolate dq/dT. - lo_bnd = (it0 == 1); - hi_bnd = (it0 + 1 == nse_table_size::ntemp); + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); - const auto [fv, dfdx] = cubic(Ts, d2, nse_table_size::dlogT, temp, lo_bnd, hi_bnd); + const auto [fv, dfdx] = cubic_split_limiter(rhos, d2, d2_deriv, + nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); - return dfdx; + // we return fv now -- it is the interpolation of dq/dT over rho. + return fv; } From 591fb3cbc3799304450a048dbd8f316689d1a6ae Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 13 Sep 2026 09:52:37 -0400 Subject: [PATCH 08/11] this fixes the T derivative --- nse_tabular/nse_table.H | 60 ++++++++++++++------- unit_test/nse_interp_check/nse_cell.H | 78 +++++++++++++++++++++++++++ 2 files changed, 120 insertions(+), 18 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 4834f425e9..dcfcc3ec3a 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -265,14 +265,34 @@ cubic(const amrex::Real* xs, const amrex::Real* fs, /// -/// an alternate version of the monotone cubic interpolation. -/// Here we use fs to compute the limiter but dfs_dx is what -/// we interpolate. +/// Differentiate the slope limiter with respect to a parameter on which +/// its input slopes depend. Select the branch using the original slopes. +/// At a branch transition, use the zero-branch derivative. +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real +limit_slope_deriv(const amrex::Real dm, const amrex::Real dp, + const amrex::Real ddm, const amrex::Real ddp) { + + if (dm * dp <= 0.0_rt) { + return 0.0_rt; + } + + const amrex::Real wm = dm / (dm + dp); + const amrex::Real wp = dp / (dm + dp); + return 2.0_rt * (wp * wp * ddm + wm * wm * ddp); +} + + +/// +/// Differentiate the cubic with respect to a parameter q at fixed x. +/// dfs_dq contains the derivatives of fs (in our application q = log10(T)). +/// Return dfdq and its spatial derivative, propagating through the limiter. /// AMREX_GPU_HOST_DEVICE AMREX_INLINE std::pair cubic_split_limiter(const amrex::Real* xs, - const amrex::Real* fs, const amrex::Real* dfs_dx, + const amrex::Real* fs, const amrex::Real* dfs_dq, const amrex::Real dx, const amrex::Real x, const bool lo_bnd, const bool hi_bnd) { @@ -287,19 +307,23 @@ cubic_split_limiter(const amrex::Real* xs, const auto dx_inv = 1.0_rt / dx; - amrex::Real m0, m1; - m1 = m0 = (fs[i0+1] - fs[i0]) * dx_inv; + const amrex::Real dc = (fs[i0+1] - fs[i0]) * dx_inv; + const amrex::Real ddc = (dfs_dq[i0+1] - dfs_dq[i0]) * dx_inv; + amrex::Real m0 = ddc; + amrex::Real m1 = ddc; - // Compute the slope from the other side and limit the slope to - // preserve monotonicity + // Differentiate the limited slopes. At a domain endpoint the slope + // is the interval secant, whose derivative is already stored above. if (! lo_bnd) { - amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; - m0 = limit_slope(d_im1, m0); + const amrex::Real d_im1 = (fs[i0] - fs[i0-1]) * dx_inv; + const amrex::Real dd_im1 = (dfs_dq[i0] - dfs_dq[i0-1]) * dx_inv; + m0 = limit_slope_deriv(d_im1, dc, dd_im1, ddc); } if (! hi_bnd) { - amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; - m1 = limit_slope(m1, d_ip1); + const amrex::Real d_ip1 = (fs[i0+2] - fs[i0+1]) * dx_inv; + const amrex::Real dd_ip1 = (dfs_dq[i0+2] - dfs_dq[i0+1]) * dx_inv; + m1 = limit_slope_deriv(dc, d_ip1, ddc, dd_ip1); } // Compute cubic Hermite polynomial basis @@ -315,7 +339,7 @@ cubic_split_limiter(const amrex::Real* xs, // Construct the interpolant // H_3(x) = f0 H_{1,0}(x) + f1 H_{1,1}(x) + f'0 Ĥ_{1,0}(x) + f'1 Ĥ_{1,0}(x) - const amrex::Real fv = dfs_dx[i0] * H0 + dfs_dx[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; + const amrex::Real fv = dfs_dq[i0] * H0 + dfs_dq[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; // now the derivative const amrex::Real dH1dx = -6.0_rt * (t2 - t) * dx_inv; @@ -323,7 +347,7 @@ cubic_split_limiter(const amrex::Real* xs, const amrex::Real dHhat1dx = 3.0_rt * t2 - 2.0_rt * t; const amrex::Real dHhat0dx = 3.0_rt * t2 - 4.0_rt * t + 1.0_rt; - const amrex::Real dfdx = dfs_dx[i0] * dH0dx + dfs_dx[i0+1] * dH1dx + + const amrex::Real dfdx = dfs_dq[i0] * dH0dx + dfs_dq[i0+1] * dH1dx + m0 * dHhat0dx + m1 * dHhat1dx; return {fv, dfdx}; @@ -465,8 +489,8 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, // the limiter), so we need to do the interpolations in the same // order (Ye -> T -> rho) since they don't commute. For a variable f, // we want to store the interpolated results of f and df/dT, and then - // for the final rho interpolation, we want to use f to compute the - // limiter for the cubic but interpolate on the values df/dT. + // for the final rho interpolation, we differentiate the limited + // slopes using both f and df/dT. const amrex::Real yes[] = {nse_table_ye(std::max(1, ic0-1)), nse_table_ye(ic0), @@ -526,8 +550,8 @@ amrex::Real tricubic_dT(const int ir0, const int it0, const int ic0, d2_deriv[ii] = dfdx; } - // finally do the remaining interpolation over rho, but we want to limit on - // q but interpolate dq/dT. + // Finally propagate dq/dT through the rho interpolation, including + // the temperature dependence of its limited slopes. lo_bnd = (ir0 == 1); hi_bnd = (ir0 + 1 == nse_table_size::nden); diff --git a/unit_test/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H index 2b4ea4ce03..979ae21d5e 100644 --- a/unit_test/nse_interp_check/nse_cell.H +++ b/unit_test/nse_interp_check/nse_cell.H @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -17,6 +18,83 @@ AMREX_INLINE void nse_cell_c() { + // Test the derivative of the 1D interpolant when its sample values + // depend on another parameter q, while the interpolation coordinate x + // stays fixed. In the NSE temperature derivative, x is log10(rho) and + // q is log10(T). These synthetic data let us test the chain rule through + // the slope limiter independently of the NSE table. + { + // Each row supplies four neighboring samples. The first two rows + // are monotone with unequal secants, so they exercise the harmonic + // mean and its derivative. The last two contain extrema, where + // opposite-sign secants cause the limiter to select zero slopes. + const amrex::Real samples[][4] = { + {0.0_rt, 1.0_rt, 3.0_rt, 6.0_rt}, + {6.0_rt, 3.0_rt, 1.0_rt, 0.0_rt}, + {0.0_rt, 2.0_rt, 1.0_rt, 3.0_rt}, + {3.0_rt, 1.0_rt, 2.0_rt, 0.0_rt}}; + // Define f_i(q) = samples_i + q * dfs_dq_i, with q = 0 as the + // reference state. Different derivatives at each sample make the + // secants, and therefore the limited slopes, depend on q. + const amrex::Real dfs_dq[] = {0.7_rt, -0.4_rt, 0.3_rt, 1.1_rt}; + + // A centered difference has O(h^2) truncation error and O(eps/h) + // roundoff error. h = eps^(1/3) balances the two for these O(1) data. + const amrex::Real h = std::cbrt(std::numeric_limits::epsilon()); + + for (const auto& fs : samples) { + // Exercise increasing coordinates and decreasing coordinates + // (the latter are used for Ye in the NSE table). + for (const auto dx : {1.0_rt, -1.0_rt}) { + const amrex::Real xs[] = {0.0_rt, dx, 2.0_rt * dx, 3.0_rt * dx}; + // Each flag replaces the corresponding endpoint's limited + // slope with the interval secant. Test neither, either, or + // both endpoints using that fallback. An outer sample is + // ignored when its endpoint flag is set. + for (const bool lo_bnd : {false, true}) { + for (const bool hi_bnd : {false, true}) { + // Evaluate at both ends of the middle interval and + // at an interior point, where slope terms contribute + // to the interpolated value as well as its derivative. + for (const auto t : {0.0_rt, 0.4_rt, 1.0_rt}) { + const amrex::Real x = xs[1] + t * dx; + // Perturb q, not x. Recomputing the full cubic at + // q = +/-h includes the change in its limiter. + // These samples stay away from limiter branch + // transitions, where a derivative may not exist. + amrex::Real fp[4], fm[4]; + for (int i = 0; i < 4; ++i) { + fp[i] = fs[i] + h * dfs_dq[i]; + fm[i] = fs[i] - h * dfs_dq[i]; + } + const auto plus = cubic(xs, fp, dx, x, lo_bnd, hi_bnd); + const auto minus = cubic(xs, fm, dx, x, lo_bnd, hi_bnd); + // Propagate the known sample derivatives through + // the unperturbed interpolant analytically. + const auto deriv = cubic_split_limiter(xs, fs, dfs_dq, + dx, x, lo_bnd, hi_bnd); + // cubic returns {f, df/dx}, so differencing its + // two outputs in q gives df/dq and d(df/dx)/dq. + // Both must agree with cubic_split_limiter. + const amrex::Real fd = (plus.first - minus.first) / (2.0_rt * h); + const amrex::Real fd_dx = (plus.second - minus.second) / (2.0_rt * h); + // Allow a margin over the expected O(h^2) error, + // with an absolute tolerance near zero and a + // relative tolerance for larger derivatives. + const auto agrees = [=] (amrex::Real a, amrex::Real b) { + return std::abs(a - b) <= 100.0_rt * h * h * + std::max(1.0_rt, std::abs(b)); + }; + if (!agrees(deriv.first, fd) || !agrees(deriv.second, fd_dx)) { + amrex::Abort("NSE limited-slope derivative regression failed"); + } + } + } + } + } + } + std::cout << "limited-slope derivative checks passed\n"; + } std::cout << std::endl; From 5f86d4082989ac2c9b8a166ac144548c5562c9c0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 14 Sep 2026 09:33:33 -0400 Subject: [PATCH 09/11] add option to force linear for species --- nse_tabular/_parameters | 3 +++ nse_tabular/nse_table.H | 9 ++++++++- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/nse_tabular/_parameters b/nse_tabular/_parameters index 3bd9a03caf..5358ae196d 100644 --- a/nse_tabular/_parameters +++ b/nse_tabular/_parameters @@ -29,3 +29,6 @@ nse_relax_factor real 1.0 # do we do tri-linear or tri-cubic interpolation on the table? nse_table_interp_linear bool 0 + +# do we always do tri-linear for species, even when doing cubic for other quantities? +nse_force_linear_interp_species bool 0 diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index dcfcc3ec3a..9b29af6034 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -727,8 +727,15 @@ void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) { if (! skip_X_fill) { for (int n = 1; n <= NumSpec; n++) { - amrex::Real _X = tricubic(ir0, it0, ic0, rholog, tlog, yet, + amrex::Real _X{}; + if (network_rp::nse_force_linear_interp_species) { + _X = trilinear(ir0, it0, ic0, rholog, tlog, yet, [=] (const int i) {return massfractab(n, i);}); + + } else { + _X = tricubic(ir0, it0, ic0, rholog, tlog, yet, + [=] (const int i) {return massfractab(n, i);}); + } nse_state.X[n-1] = amrex::Clamp(_X, 0.0_rt, 1.0_rt); } } From fee6e5164ed4f50c2e4912d7e2e0022177703be5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 14 Sep 2026 09:45:08 -0400 Subject: [PATCH 10/11] update benchmark --- .../aprox19_NSE_state_over_time.txt | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_NSE_state_over_time.txt b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_NSE_state_over_time.txt index d8bc2940d4..1b367f7b6d 100644 --- a/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_NSE_state_over_time.txt +++ b/unit_test/burn_cell_sdc/ci-benchmarks/aprox19_NSE_state_over_time.txt @@ -1,13 +1,13 @@ # Time Density Temperature H1 He3 He4 C12 N14 O16 Ne20 Mg24 Si28 S32 Ar36 Ca40 Ti44 Cr48 Fe52 Fe54 Ni56 n p 0 1e+07 5e+09 0.1 0.025 0.8 0.025 0.025 0.025 0 0 0 0 0 0 0 0 0 0 0 0 0 0.000909091 1e+07 5.16542e+09 0.10647 9.90099e-31 0.788419 1.34585e-05 0.0164549 0.0274406 0.0173537 0.0125208 0.0259179 0.00434822 0.00094121 0.000118497 1.87907e-06 1.24233e-07 5.57196e-09 9.31398e-11 7.88442e-15 3.34056e-14 3.48303e-12 -0.00115544 1e+07 5.38284e+09 1e-30 3.03736e-10 0.135293 1.10373e-06 8.70507e-11 2.82863e-06 4.35126e-08 9.70415e-06 0.00582948 0.00645566 0.00419711 0.00477662 0.00131915 0.0375917 0.0582282 0.669506 0.0501732 1.57018e-06 0.0266143 -0.00146854 1e+07 5.38284e+09 1e-30 3.03735e-10 0.135293 1.10373e-06 8.70506e-11 2.82862e-06 4.35126e-08 9.70415e-06 0.00582948 0.00645566 0.00419711 0.00477663 0.00131915 0.0375917 0.0582282 0.669506 0.0501733 1.57018e-06 0.0266143 -0.00186649 1e+07 5.38284e+09 1e-30 3.03735e-10 0.135293 1.10373e-06 8.70504e-11 2.82862e-06 4.35125e-08 9.70414e-06 0.00582948 0.00645566 0.00419711 0.00477663 0.00131915 0.0375917 0.0582281 0.669506 0.0501733 1.57018e-06 0.0266143 -0.00237227 1e+07 5.38284e+09 1e-30 3.03735e-10 0.135293 1.10373e-06 8.70502e-11 2.82862e-06 4.35124e-08 9.70413e-06 0.00582948 0.00645566 0.00419711 0.00477663 0.00131915 0.0375917 0.0582281 0.669507 0.0501734 1.57017e-06 0.0266143 -0.00301511 1e+07 5.38284e+09 1e-30 3.03734e-10 0.135292 1.10373e-06 8.705e-11 2.82861e-06 4.35124e-08 9.70412e-06 0.00582948 0.00645566 0.00419712 0.00477664 0.00131915 0.0375917 0.0582281 0.669507 0.0501736 1.57017e-06 0.0266143 -0.00383215 1e+07 5.38284e+09 1e-30 3.03733e-10 0.135292 1.10373e-06 8.70497e-11 2.82861e-06 4.35123e-08 9.70411e-06 0.00582948 0.00645566 0.00419712 0.00477664 0.00131915 0.0375917 0.0582281 0.669507 0.0501737 1.57016e-06 0.0266143 - 0.0048706 1e+07 5.38284e+09 1e-30 3.03732e-10 0.135292 1.10372e-06 8.70493e-11 2.8286e-06 4.35121e-08 9.70409e-06 0.00582948 0.00645567 0.00419712 0.00477665 0.00131915 0.0375916 0.0582281 0.669507 0.0501739 1.57015e-06 0.0266143 -0.00619044 1e+07 5.38284e+09 1e-30 3.03731e-10 0.135292 1.10372e-06 8.70487e-11 2.82859e-06 4.3512e-08 9.70407e-06 0.00582947 0.00645567 0.00419712 0.00477666 0.00131915 0.0375916 0.0582281 0.669507 0.0501742 1.57014e-06 0.0266143 -0.00786793 1e+07 5.38284e+09 1e-30 3.03729e-10 0.135291 1.10371e-06 8.70481e-11 2.82858e-06 4.35118e-08 9.70404e-06 0.00582947 0.00645567 0.00419713 0.00477667 0.00131915 0.0375916 0.0582281 0.669507 0.0501745 1.57013e-06 0.0266142 - 0.01 1e+07 5.38284e+09 1e-30 3.03727e-10 0.135291 1.1037e-06 8.70473e-11 2.82856e-06 4.35115e-08 9.704e-06 0.00582947 0.00645568 0.00419714 0.00477669 0.00131914 0.0375915 0.058228 0.669507 0.0501749 1.57011e-06 0.0266142 +0.00115544 1e+07 5.38284e+09 9.99801e-31 3.03612e-10 0.135259 1.10329e-06 8.70011e-11 2.82759e-06 4.34945e-08 9.70124e-06 0.00582711 0.00645595 0.00419692 0.0047762 0.0013188 0.0375803 0.0582025 0.669578 0.0501823 1.57046e-06 0.026608 +0.00146854 1e+07 5.38284e+09 9.99801e-31 3.03611e-10 0.135259 1.10329e-06 8.7001e-11 2.82759e-06 4.34944e-08 9.70124e-06 0.00582711 0.00645595 0.00419693 0.0047762 0.0013188 0.0375803 0.0582025 0.669578 0.0501824 1.57046e-06 0.026608 +0.00186649 1e+07 5.38284e+09 9.99801e-31 3.03611e-10 0.135258 1.10328e-06 8.70008e-11 2.82759e-06 4.34944e-08 9.70123e-06 0.00582711 0.00645595 0.00419693 0.0047762 0.0013188 0.0375803 0.0582025 0.669578 0.0501825 1.57045e-06 0.026608 +0.00237227 1e+07 5.38284e+09 9.99801e-31 3.03611e-10 0.135258 1.10328e-06 8.70006e-11 2.82758e-06 4.34943e-08 9.70122e-06 0.00582711 0.00645595 0.00419693 0.0047762 0.0013188 0.0375803 0.0582025 0.669578 0.0501826 1.57045e-06 0.026608 +0.00301511 1e+07 5.38284e+09 9.99801e-31 3.0361e-10 0.135258 1.10328e-06 8.70004e-11 2.82758e-06 4.34943e-08 9.70121e-06 0.00582711 0.00645596 0.00419693 0.00477621 0.0013188 0.0375803 0.0582025 0.669578 0.0501827 1.57044e-06 0.026608 +0.00383215 1e+07 5.38284e+09 9.99801e-31 3.03609e-10 0.135258 1.10328e-06 8.70001e-11 2.82757e-06 4.34942e-08 9.7012e-06 0.00582711 0.00645596 0.00419693 0.00477621 0.00131879 0.0375803 0.0582024 0.669578 0.0501829 1.57044e-06 0.026608 + 0.0048706 1e+07 5.38284e+09 9.99801e-31 3.03608e-10 0.135258 1.10327e-06 8.69997e-11 2.82757e-06 4.3494e-08 9.70118e-06 0.00582712 0.00645596 0.00419694 0.00477622 0.00131879 0.0375802 0.0582024 0.669578 0.0501831 1.57043e-06 0.026608 +0.00619044 1e+07 5.38284e+09 9.99801e-31 3.03607e-10 0.135257 1.10327e-06 8.69992e-11 2.82756e-06 4.34939e-08 9.70116e-06 0.00582712 0.00645596 0.00419694 0.00477623 0.00131879 0.0375802 0.0582024 0.669578 0.0501833 1.57042e-06 0.026608 +0.00786793 1e+07 5.38284e+09 9.99801e-31 3.03605e-10 0.135257 1.10327e-06 8.69986e-11 2.82755e-06 4.34937e-08 9.70114e-06 0.00582712 0.00645597 0.00419695 0.00477624 0.00131879 0.0375802 0.0582024 0.669579 0.0501836 1.5704e-06 0.026608 + 0.01 1e+07 5.38284e+09 9.99801e-31 3.03603e-10 0.135256 1.10326e-06 8.69978e-11 2.82753e-06 4.34934e-08 9.7011e-06 0.00582712 0.00645598 0.00419695 0.00477626 0.00131879 0.0375802 0.0582024 0.669579 0.050184 1.57039e-06 0.0266079 From abbb87c389ff422e1d8ef316601e8f84bb070fb3 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 14 Sep 2026 09:52:36 -0400 Subject: [PATCH 11/11] update benchmark --- .../ci-benchmarks/aprox19.out | 117 +++++++++--------- 1 file changed, 59 insertions(+), 58 deletions(-) diff --git a/unit_test/nse_interp_check/ci-benchmarks/aprox19.out b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out index 53a324c615..eb68805dd4 100644 --- a/unit_test/nse_interp_check/ci-benchmarks/aprox19.out +++ b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out @@ -2,6 +2,7 @@ Initializing AMReX (26.09-119-ge60cdc18711c)... AMReX (26.09-119-ge60cdc18711c) initialized starting the single zone burn... reading the NSE table (C++) ... +limited-slope derivative checks passed rho, T, Ye = 1230000000 5180000000 0.472 @@ -12,109 +13,109 @@ 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 + 9.000 55.53344914 + 9.050 55.58956045 + 9.100 55.64268998 + 9.150 55.69324501 - cubic interpolated value: 55.6078922 + cubic interpolated value: 55.63218008 4 T values: - 9.700 55.5223665 - 9.710 55.39489494 - 9.720 55.22387568 - 9.730 54.99343581 + 9.700 55.69817527 + 9.710 55.58956045 + 9.720 55.44641534 + 9.730 55.25607232 - cubic interpolated value: 55.32712406 + cubic interpolated value: 55.53246868 4 Ye values: - 0.475 55.5223665 - 0.472 55.65227675 - 0.470 55.78142906 - 0.468 55.90244371 + 0.475 55.45473797 + 0.472 55.58956045 + 0.470 55.72339147 + 0.468 55.85036642 - cubic interpolated value: 55.678404 + cubic interpolated value: 55.61650246 tricubic interpolated values at input state - abar = 55.60652462 - bea = 8.755119523 - dyedt = -0.004527880871 - dbeadt = 0.000111945508 - e_nu = 1.101641642e+16 + abar = 55.60635196 + bea = 8.755115136 + dyedt = -0.004529430118 + dbeadt = 0.0001120934856 + e_nu = 1.102105121e+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 + X(He3) = 2.773928062e-14 + X(He4) = 0.0003164637043 + X(C12) = 8.009439378e-10 + X(N14) = 3.158308123e-14 + X(O16) = 2.004728182e-09 + X(Ne20) = 2.46116577e-11 + X(Mg24) = 1.031696084e-08 + X(Si28) = 8.185902305e-06 + X(S32) = 1.196632498e-05 + X(Ar36) = 1.246184748e-05 + X(Ca40) = 1.267413754e-05 + X(Ti44) = 0.0001438375571 + X(Cr48) = 0.01147226672 + X(Fe52) = 0.07746064144 + X(Fe54) = 0.9098518462 + X(Ni56) = 0.000725722646 + X(n) = 2.087009888e-08 + X(p) = 9.852136026e-06 testing temperature derivatives of cubic first finite-difference derivatives - dAbar/dT = -1.072562355e-09 - dbea/dT = -6.867059575e-12 + dAbar/dT = -1.075182173e-09 + dbea/dT = -6.863870363e-12 now using derivative of the interpolant - dAbar/dT = -1.072562604e-09 - dbea/dT = -6.867048589e-12 + dAbar/dT = -1.075182436e-09 + dbea/dT = -6.863845759e-12 testing density derivatives of cubic first finite-difference derivatives - dAbar/drho = 3.987525411e-10 - dbea/drho = 7.619559936e-13 + dAbar/drho = 3.988131972e-10 + dbea/drho = 7.618115743e-13 now using derivative of the interpolant - dAbar/drho = 3.987522836e-10 - dbea/drho = 7.618831514e-13 + dAbar/drho = 3.988138616e-10 + dbea/drho = 7.618659448e-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 + change in e: 1.395277041e+18 1.388449234e+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 + updated T: 6394534595 + change in abar: 55.60635196 50.26816014 + change in e: 1.388449234e+18 1.388449234e+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 + change in p: 6.622157915e+26 6.577851483e+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 + updated T: 6466780918 + change in abar: 55.60635196 49.50402292 + change in p: 6.577851483e+26 6.577851483e+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 + change in p: 6.577759303e+26 6.577851483e+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 + change in abar: 55.60635196 55.62478096 + change in p: 6.577851483e+26 6.577851483e+26 explore a troublesome state @@ -122,7 +123,7 @@ 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 + cubic interpolation gives: e_nu = 5.481706469e+12 bounds of e_nu in the 64 cells used for tricubic: 1.3703405e+11 6.2562515e+14 AMReX (26.09-119-ge60cdc18711c) finalized