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/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/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 26aeecf3e7..9b29af6034 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -175,48 +176,181 @@ 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) { +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 + + // 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); + } - // 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] + 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 = fs[i0] * H0 + fs[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; - 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]; + // 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; - return a * amrex::Math::powi<3>(x - xs[1]) + - b * amrex::Math::powi<2>(x - xs[1]) + c * (x - xs[1]) + d; + const amrex::Real dfdx = fs[i0] * dH0dx + fs[i0+1] * dH1dx + + m0 * dHhat0dx + m1 * dHhat1dx; + + return {fv, dfdx}; } + /// -/// 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 +/// 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 cubic_deriv(const amrex::Real* xs, const amrex::Real* fs, const amrex::Real dx, const amrex::Real x) { +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_dq, + 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; + + 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; + + // Differentiate the limited slopes. At a domain endpoint the slope + // is the interval secant, whose derivative is already stored above. + if (! lo_bnd) { + 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) { + 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 + 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_dq[i0] * H0 + dfs_dq[i0+1] * H1 + m0 * Hhat0 + m1 * Hhat1; - // 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 + // 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; - 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 = dfs_dq[i0] * dH0dx + dfs_dq[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 {fv, dfdx}; } @@ -267,37 +401,49 @@ 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); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; } } @@ -305,17 +451,23 @@ 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); + 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 - amrex::Real val = cubic(rhos, d2, nse_table_size::dlogrho, rho); + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); - return val; + const auto [fv, dfdx] = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + + return fv; } @@ -328,56 +480,87 @@ 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. + + // 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 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), 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); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; } } - // now do the 4 rho interpolations (one in each T plane) + // now do the 4 T interpolations (one in each rho plane) + + // we need to store both q and dq/dT interpolated in the T + // direction in all 4 rho planes. 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]}; - d2[jj] = cubic(rhos, _d, nse_table_size::dlogrho, rho); + 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 propagate dq/dT through the rho interpolation, including + // the temperature dependence of its limited slopes. - amrex::Real val = cubic_deriv(Ts, d2, nse_table_size::dlogT, temp); + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); - return val; + const auto [fv, dfdx] = cubic_split_limiter(rhos, d2, d2_deriv, + nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + + // we return fv now -- it is the interpolation of dq/dT over rho. + return fv; } @@ -391,37 +574,50 @@ 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); + const auto [fv, dfdx] = cubic(yes, _d, -nse_table_size::dye, ye, lo_bnd, hi_bnd); + d1[ii][jj] = fv; + } } @@ -429,18 +625,26 @@ 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); + 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 // the derivative of the interpolant - amrex::Real val = cubic_deriv(rhos, d2, nse_table_size::dlogrho, rho); - return val; + lo_bnd = (ir0 == 1); + hi_bnd = (ir0 + 1 == nse_table_size::nden); + + const auto [fv, dfdx] = cubic(rhos, d2, nse_table_size::dlogrho, rho, lo_bnd, hi_bnd); + + return dfdx; } @@ -502,18 +706,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 + // 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 ir0 = nse_get_logrho_index(rholog) - 1; - ir0 = amrex::Clamp(ir0, 1, nse_table_size::nden-3); + // we will handle the case where it0 or it0+1 is at the bounds of + // the table in the caller. - 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); nse_state.abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab); nse_state.bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab); @@ -526,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); } } @@ -569,14 +777,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 @@ -622,14 +825,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/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 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..eb68805dd4 --- /dev/null +++ b/unit_test/nse_interp_check/ci-benchmarks/aprox19.out @@ -0,0 +1,129 @@ +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 + +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.53344914 + 9.050 55.58956045 + 9.100 55.64268998 + 9.150 55.69324501 + + cubic interpolated value: 55.63218008 + + 4 T values: + 9.700 55.69817527 + 9.710 55.58956045 + 9.720 55.44641534 + 9.730 55.25607232 + + cubic interpolated value: 55.53246868 + + 4 Ye values: + 0.475 55.45473797 + 0.472 55.58956045 + 0.470 55.72339147 + 0.468 55.85036642 + + cubic interpolated value: 55.61650246 + + +tricubic interpolated values at input state + + abar = 55.60635196 + bea = 8.755115136 + dyedt = -0.004529430118 + dbeadt = 0.0001120934856 + e_nu = 1.102105121e+16 + + X(H1) = 0 + 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.075182173e-09 + dbea/dT = -6.863870363e-12 + + now using derivative of the interpolant + dAbar/dT = -1.075182436e-09 + dbea/dT = -6.863845759e-12 + +testing density derivatives of cubic + + first finite-difference derivatives + dAbar/drho = 3.988131972e-10 + dbea/drho = 7.618115743e-13 + + now using derivative of the interpolant + 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.395277041e+18 1.388449234e+18 + + NSE method: use the nse_T_abar_from_e() function: + 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.622157915e+26 6.577851483e+26 + + NSE method: use the nse_T_abar_from_p() function: + 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.577759303e+26 6.577851483e+26 + + NSE method: use the nse_rho_abar_from_p() function: + updated T: 5180000000 + change in abar: 55.60635196 55.62478096 + change in p: 6.577851483e+26 6.577851483e+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 = 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 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/nse_interp_check/nse_cell.H b/unit_test/nse_interp_check/nse_cell.H new file mode 100644 index 0000000000..979ae21d5e --- /dev/null +++ b/unit_test/nse_interp_check/nse_cell.H @@ -0,0 +1,620 @@ +#ifndef NSE_CELL_H +#define NSE_CELL_H + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +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; + + std::cout << "rho, T, Ye = " + << unit_test_rp::density << " " + << unit_test_rp::temperature << " " + << unit_test_rp::ye << std::endl << std::endl; + + // check the indices + + amrex::Real logrho = std::log10(unit_test_rp::density); + amrex::Real logT = std::log10(unit_test_rp::temperature); + + int ir0 = nse_get_logrho_index(logrho); + int it0 = nse_get_logT_index(logT); + int ic0 = nse_get_ye_index(unit_test_rp::ye); + + std::cout << "density value brackets: " + << nse_table_logrho(ir0) << " < " << logrho << " < " + << nse_table_logrho(ir0+1) << std::endl; + + std::cout << "temperature value brackets: " + << nse_table_logT(it0) << " < " << logT << " < " + << nse_table_logT(it0+1) << std::endl; + + std::cout << "Ye value brackets: " + << nse_table_ye(ic0) << " < " << unit_test_rp::ye << " < " + << nse_table_ye(ic0+1) << std::endl; + + std::cout << std::endl; + + // 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"; + + { + + // density interpolation + // note: assuming we are not at a boundary + + amrex::Real rhos[] = {nse_table_logrho(ir0-1), + nse_table_logrho(ir0), + nse_table_logrho(ir0+1), + nse_table_logrho(ir0+2)}; + + std::cout << " 4 rho values: " << std::endl; + + 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))}; + + for (int ii = 0; ii < 4; ++ii) { + std::cout << std::format(" {:5.3f} {:15.10g}\n", rhos[ii], _d[ii]); + } + std::cout << 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; + } + + { + + // temperature interpolation + + amrex::Real Ts[] = {nse_table_logT(it0-1), + nse_table_logT(it0), + nse_table_logT(it0+1), + nse_table_logT(it0+2)}; + + std::cout << " 4 T values: " << std::endl; + + 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))}; + + for (int ii = 0; ii < 4; ++ii) { + std::cout << std::format(" {:5.3f} {:15.10g}\n", Ts[ii], _d[ii]); + } + std::cout << 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; + + } + + { + + // Ye interpolation + + amrex::Real yes[] = {nse_table_ye(ic0-1), + nse_table_ye(ic0), + nse_table_ye(ic0+1), + nse_table_ye(ic0+2)}; + + std::cout << " 4 Ye values: " << std::endl; + + 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))}; + + for (int ii = 0; ii < 4; ++ii) { + std::cout << std::format(" {:5.3f} {:15.10g}\n", yes[ii], _d[ii]); + } + std::cout << 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; + + } + + 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; + nse_state.rho = unit_test_rp::density; + nse_state.Ye = unit_test_rp::ye; + + 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 << std::endl; + for (int n = 0; n < NumSpec; ++n) { + std::cout << " X(" << short_spec_names_cxx[n] << ") = " << nse_state.X[n] << std::endl; + } + + + // temp derivatives + + std::cout << std::endl; + std::cout << std::endl; + std::cout << "testing temperature derivatives of cubic" << std::endl << std::endl; + + std::cout << " first finite-difference derivatives\n"; + + nse_state.T = unit_test_rp::temperature; + nse_state.rho = unit_test_rp::density; + nse_state.Ye = unit_test_rp::ye; + + nse_interp(nse_state); + + amrex::Real abar_old = nse_state.abar; + amrex::Real bea_old = nse_state.bea; + amrex::Real T_old = nse_state.T; + + const amrex::Real eps = 1.e-8_rt; + nse_state.T *= (1.0_rt + eps); + + 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 << 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); + + 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 << std::endl; + + // dens derivatives + + std::cout << "testing density derivatives of cubic" << std::endl << std::endl; + + std::cout << " first finite-difference derivatives" << std::endl; + + nse_state.T = unit_test_rp::temperature; + nse_state.rho = unit_test_rp::density; + nse_state.Ye = unit_test_rp::ye; + + nse_interp(nse_state); + + amrex::Real rho_old = nse_state.rho; + + nse_state.rho *= (1.0_rt + eps); + + 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 << 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); + + 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 << std::endl; + std::cout << std::endl; + + // + // EOS testing + // + + // now we test the EOS, in particular, we want to ensure that + // given e, rho, Y_e, we can find a T that is consistent with both + // the EOS and the NSE table + + // attempt 1: using the normal EOS interfaces + // + // we'll do: + // + // T, rho, Ye -> Abar + // T, rho, Ye, Abar -> e + // + // then perturb e to e' and do: + // + // e', rho, Ye -> T' + // T', rho, Ye -> Abar' + // + // and finally check + // + // T', rho, Ye, Abar' -> e' ??? + + { + + 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; + nse_state.rho = unit_test_rp::density; + nse_state.Ye = unit_test_rp::ye; + + nse_interp(nse_state); + + amrex::Real abar_orig = nse_state.abar; + + eos_t eos_state; + + eos_state.T = unit_test_rp::temperature; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; + + eos(eos_input_rt, eos_state); + + amrex::Real e_orig = eos_state.e; + + // now perturb e and find T, then redo NSE to get abar and finally + // see if we get back our new e + + amrex::Real e_new = eos_state.e * 1.05; + + eos_state.e = e_new; + eos(eos_input_re, eos_state); + + nse_state.T = eos_state.T; + nse_interp(nse_state); + + eos_state.aux[iabar] = nse_state.abar; + eos(eos_input_rt, eos_state); + + 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: + // e', rho, Ye -> Abar', T' + + eos_state.T = unit_test_rp::temperature; + eos_state.e = e_new; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = unit_test_rp::ye; + eos_state.aux[iabar] = abar_orig; + + amrex::Real abar_start = eos_state.aux[iabar]; + + nse_T_abar_from_e(eos_state.rho, eos_state.e, eos_state.aux[iye], + eos_state.T, eos_state.aux[iabar]); + + 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 << " change in e: " << eos_state.e << " " << e_new << std::endl; + std::cout << std::endl; + } + + + { + + // 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; + + nse_interp(nse_state); + + amrex::Real abar_orig = nse_state.abar; + + eos_t eos_state; + + eos_state.T = unit_test_rp::temperature; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; + + eos(eos_input_rt, eos_state); + + amrex::Real p_orig = eos_state.p; + + // now perturb p and find T, 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; + + eos_state.p = p_new; + eos(eos_input_rp, eos_state); + + nse_state.T = eos_state.T; + nse_interp(nse_state); + + eos_state.aux[iabar] = nse_state.abar; + eos(eos_input_rt, eos_state); + + 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: + // p', rho, Ye -> Abar', T' + + eos_state.T = unit_test_rp::temperature; + eos_state.p = p_new; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = unit_test_rp::ye; + eos_state.aux[iabar] = abar_orig; + + amrex::Real abar_start = eos_state.aux[iabar]; + + nse_T_abar_from_p(eos_state.rho, eos_state.p, eos_state.aux[iye], + eos_state.T, eos_state.aux[iabar]); + + 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 << " change in p: " << eos_state.p << " " << p_new << std::endl; + std::cout << std::endl; + } + + + // now we test the EOS inversion for finding rho given e or p. The idea + // is the same as above. We don't have an eos_input_te, so we will only + // check p + + { + + 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; + + nse_interp(nse_state); + + amrex::Real abar_orig = nse_state.abar; + + eos_t eos_state; + + eos_state.T = unit_test_rp::temperature; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; + + eos(eos_input_rt, eos_state); + + amrex::Real p_orig = eos_state.p; + + // 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; + + eos_state.p = p_new; + eos(eos_input_tp, eos_state); + + nse_state.rho = eos_state.rho; + nse_interp(nse_state); + + eos_state.aux[iabar] = nse_state.abar; + eos(eos_input_rt, eos_state); + + 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: + // p', T, Ye -> Abar', rho' + + eos_state.T = unit_test_rp::temperature; + eos_state.p = p_new; + eos_state.rho = unit_test_rp::density; + eos_state.aux[iye] = unit_test_rp::ye; + eos_state.aux[iabar] = abar_orig; + + amrex::Real abar_start = eos_state.aux[iabar]; + + nse_rho_abar_from_p(eos_state.T, eos_state.p, eos_state.aux[iye], + eos_state.rho, eos_state.aux[iabar]); + + 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 << " 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 diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H deleted file mode 100644 index c1eb3cf3f1..0000000000 --- a/unit_test/test_nse_interp/nse_cell.H +++ /dev/null @@ -1,445 +0,0 @@ -#ifndef NSE_CELL_H -#define NSE_CELL_H - -#include -#include -#include -#include -#include - -#include -#include -#include - -AMREX_INLINE -void nse_cell_c() -{ - - std::cout << "rho, T, Ye = " - << unit_test_rp::density << " " - << unit_test_rp::temperature << " " - << unit_test_rp::ye << std::endl; - - // check the indices - - amrex::Real logrho = std::log10(unit_test_rp::density); - amrex::Real logT = std::log10(unit_test_rp::temperature); - - int ir0 = nse_get_logrho_index(logrho); - int it0 = nse_get_logT_index(logT); - int ic0 = nse_get_ye_index(unit_test_rp::ye); - - std::cout << "density value brackets: " - << nse_table_logrho(ir0) << " < " << logrho << " < " - << nse_table_logrho(ir0+1) << std::endl; - - std::cout << "temperature value brackets: " - << nse_table_logT(it0) << " < " << logT << " < " - << nse_table_logT(it0+1) << std::endl; - - std::cout << "Ye value brackets: " - << nse_table_ye(ic0) << " < " << unit_test_rp::ye << " < " - << nse_table_ye(ic0+1) << std::endl; - - std::cout << std::endl; - - // now try 1-d interpolation in T -- for cubic interpolation, we - // need to offset the indices one to the left - - 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)); - - - 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)}; - - - { - 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)), - nse_table::abartab(nse_idx(ir0+2, it0, ic0)), - 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::endl; - - std::cout << "cubic interpolated value: " << - cubic(rhos, _d, nse_table_size::dlogrho, logrho) << std::endl << std::endl; - } - - { - 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)), - nse_table::abartab(nse_idx(ir0, it0+2, ic0)), - 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::endl; - - std::cout << "cubic interpolated value: " << - cubic(Ts, _d, nse_table_size::dlogT, logT) << std::endl << std::endl; - - } - - { - 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)), - nse_table::abartab(nse_idx(ir0, it0, ic0+2)), - 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::endl; - - 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; - - nse_table_t nse_state; - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - 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; - for (int n = 0; n < NumSpec; ++n) { - 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 << "first finite-difference derivatives" << std::endl; - - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - nse_interp(nse_state); - - amrex::Real abar_old = nse_state.abar; - amrex::Real bea_old = nse_state.bea; - amrex::Real T_old = nse_state.T; - - const amrex::Real eps = 1.e-8_rt; - nse_state.T *= (1.0_rt + eps); - - 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 << "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); - - 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 << std::endl; - - // dens derivatives - - std::cout << std::endl; - std::cout << "testing density derivatives of cubic" << std::endl; - - std::cout << "first finite-difference derivatives" << std::endl; - - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - nse_interp(nse_state); - - amrex::Real rho_old = nse_state.rho; - - nse_state.rho *= (1.0_rt + eps); - - 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 << "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); - - 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 << std::endl; - - // - // EOS testing - // - - // now we test the EOS, in particular, we want to ensure that - // given e, rho, Y_e, we can find a T that is consistent with both - // the EOS and the NSE table - - // attempt 1: using the normal EOS interfaces - // - // we'll do: - // - // T, rho, Ye -> Abar - // T, rho, Ye, Abar -> e - // - // then perturb e to e' and do: - // - // e', rho, Ye -> T' - // T', rho, Ye -> Abar' - // - // and finally check - // - // T', rho, Ye, Abar' -> e' ??? - - { - - // first get the abar consistent with our inputs and find e - - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - nse_interp(nse_state); - - amrex::Real abar_orig = nse_state.abar; - - eos_t eos_state; - - eos_state.T = unit_test_rp::temperature; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = nse_state.Ye; - eos_state.aux[iabar] = nse_state.abar; - - eos(eos_input_rt, eos_state); - - amrex::Real e_orig = eos_state.e; - - // now perturb e and find T, then redo NSE to get abar and finally - // see if we get back our new e - - amrex::Real e_new = eos_state.e * 1.05; - - eos_state.e = e_new; - eos(eos_input_re, eos_state); - - nse_state.T = eos_state.T; - nse_interp(nse_state); - - 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; - - // attempt 2: - // now we try the new interface. This effectively does: - // e', rho, Ye -> Abar', T' - - eos_state.T = unit_test_rp::temperature; - eos_state.e = e_new; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = unit_test_rp::ye; - eos_state.aux[iabar] = abar_orig; - - amrex::Real abar_start = eos_state.aux[iabar]; - - 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; - - // 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 << std::endl; - } - - - { - - // now redo it for pressure - - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - nse_interp(nse_state); - - amrex::Real abar_orig = nse_state.abar; - - eos_t eos_state; - - eos_state.T = unit_test_rp::temperature; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = nse_state.Ye; - eos_state.aux[iabar] = nse_state.abar; - - eos(eos_input_rt, eos_state); - - amrex::Real p_orig = eos_state.p; - - // now perturb p and find T, 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; - - eos_state.p = p_new; - eos(eos_input_rp, eos_state); - - nse_state.T = eos_state.T; - nse_interp(nse_state); - - 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; - - // attempt 2: - // now we try the new interface. This effectively does: - // p', rho, Ye -> Abar', T' - - eos_state.T = unit_test_rp::temperature; - eos_state.p = p_new; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = unit_test_rp::ye; - eos_state.aux[iabar] = abar_orig; - - amrex::Real abar_start = eos_state.aux[iabar]; - - 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; - - // 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 << std::endl; - } - - - // now we test the EOS inversion for finding rho given e or p. The idea - // is the same as above. We don't have an eos_input_te, so we will only - // check p - - { - - nse_state.T = unit_test_rp::temperature; - nse_state.rho = unit_test_rp::density; - nse_state.Ye = unit_test_rp::ye; - - nse_interp(nse_state); - - amrex::Real abar_orig = nse_state.abar; - - eos_t eos_state; - - eos_state.T = unit_test_rp::temperature; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = nse_state.Ye; - eos_state.aux[iabar] = nse_state.abar; - - eos(eos_input_rt, eos_state); - - amrex::Real p_orig = eos_state.p; - - // now perturb p and find T, 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; - - eos_state.p = p_new; - eos(eos_input_tp, eos_state); - - nse_state.rho = eos_state.rho; - nse_interp(nse_state); - - 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; - - // attempt 2: - // now we try the new interface. This effectively does: - // p', T, Ye -> Abar', rho' - - eos_state.T = unit_test_rp::temperature; - eos_state.p = p_new; - eos_state.rho = unit_test_rp::density; - eos_state.aux[iye] = unit_test_rp::ye; - eos_state.aux[iabar] = abar_orig; - - amrex::Real abar_start = eos_state.aux[iabar]; - - 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; - - // 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 << std::endl; - } - - -} -#endif