From 23eb8907781d6a5f1c025e812ecd6ca65df9260c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 11 Dec 2023 12:48:25 -0500 Subject: [PATCH 1/4] start working on temperature derivatives --- nse_tabular/nse_table.H | 148 ++++++++++++++++++++++++++- unit_test/test_nse_interp/nse_cell.H | 6 ++ 2 files changed, 151 insertions(+), 3 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 933b05f8dd..4368e65984 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -119,6 +119,10 @@ int nse_get_ye_index(const Real ye) { return ic0 + 1; } +/// +/// given 4 points (xs, fs), with spacing dx, return the interplated +/// value of f at point x by fitting a cubic to the points +/// AMREX_GPU_HOST_DEVICE AMREX_INLINE Real cubic(const Real* xs, const Real* fs, const Real dx, const Real x) { @@ -136,6 +140,30 @@ Real cubic(const Real* xs, const Real* fs, const Real dx, const Real x) { } +/// +/// 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 +Real cubic_deriv(const Real* xs, const Real* fs, const Real dx, const 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 + + Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * std::pow(dx, 3)); + Real b = (-2 * fs[1] + fs[2] + fs[0]) / (2 * dx * dx); + Real c = (-3 * fs[1] + 6 * fs[2] - fs[3] - 2 * fs[0]) / (6 * dx); + //Real d = fs[1]; + + return 3.0_rt * a * std::pow(x - xs[1], 2) + 2.0_rt * b * (x - xs[1]) + c; + +} + + template AMREX_GPU_HOST_DEVICE AMREX_INLINE Real trilinear(const int ir1, const int it1, const int ic1, @@ -234,6 +262,68 @@ Real tricubic(const int ir0, const int it0, const int ic0, } +/// +/// take the temperature derivative of a table quantity by differentiating +/// the cubic interpolant +/// +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +Real tricubic_dT(const int ir0, const int it0, const int ic0, + const Real rho, const Real temp, const Real ye, const T& data) { + + Real yes[] = {nse_table_ye(ic0), + nse_table_ye(ic0+1), + nse_table_ye(ic0+2), + nse_table_ye(ic0+3)}; + + Real Ts[] = {nse_table_logT(it0), + nse_table_logT(it0+1), + nse_table_logT(it0+2), + nse_table_logT(it0+3)}; + + Real rhos[] = {nse_table_logrho(ir0), + nse_table_logrho(ir0+1), + nse_table_logrho(ir0+2), + nse_table_logrho(ir0+3)}; + + // first do the 16 ye interpolations + + // the first index will be rho and the second will be T + Real d1[4][4]; + + for (int ii = 0; ii < 4; ++ii) { + for (int jj = 0; jj < 4; ++jj) { + + 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))}; + + // 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); + } + } + + // now do the 4 rho interpolations (one in each T plane) + + Real d2[4]; + + for (int jj = 0; jj < 4; ++jj) { + + Real _d[] = {d1[0][jj], d1[1][jj], d1[2][jj], d1[3][jj]}; + d2[jj] = cubic(rhos, _d, nse_table_size::dlogrho, rho); + } + + // finally do the remaining interpolation over T, but return + // the derivative of the interpolant + + Real val = cubic_deriv(Ts, d2, nse_table_size::dlogT, temp); + + return val; + +} + AMREX_GPU_HOST_DEVICE AMREX_INLINE void nse_interp(const Real T, const Real rho, const Real ye, Real& abar, Real& bea, Real& dyedt, Real& dabardt, Real& dbeadt, Real& e_nu, @@ -300,13 +390,13 @@ void nse_interp(const Real T, const Real rho, const Real ye, // 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::max(1, amrex::min(nse_table_size::nden-3, ir0)); + ir0 = std::clamp(ir0, 1, nse_table_size::nden-3); int it0 = nse_get_logT_index(tlog) - 1; - it0 = amrex::max(1, amrex::min(nse_table_size::ntemp-3, it0)); + it0 = std::clamp(it0, 1, nse_table_size::ntemp-3); int ic0 = nse_get_ye_index(yet) - 1; - ic0 = amrex::max(1, amrex::min(nse_table_size::nye-3, ic0)); + ic0 = std::clamp(ic0, 1, nse_table_size::nye-3); abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab); bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab); @@ -329,4 +419,56 @@ void nse_interp(const Real T, const Real rho, const Real ye, } + +/// +/// compute the temperature derivative of the table quantity data +/// at the point T, rho, ye by using cubic interpolation +/// +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +Real +nse_interp_dT(const Real temp, const Real rho, const Real ye, const T& data) { + + Real rholog = std::log10(rho); + { + Real rmin = nse_table_size::logrho_min; + Real rmax = nse_table_size::logrho_max; + + rholog = std::clamp(rholog, rmin, rmax); + } + + Real tlog = std::log10(temp); + { + Real tmin = nse_table_size::logT_min; + Real tmax = nse_table_size::logT_max; + + tlog = std::clamp(tlog, tmin, tmax); + } + + Real yet = ye; + { + Real yemin = nse_table_size::ye_min; + Real yemax = nse_table_size::ye_max; + + yet = std::clamp(yet, yemin, yemax); + } + + int ir0 = nse_get_logrho_index(rholog) - 1; + ir0 = std::clamp(ir0, 1, nse_table_size::nden-3); + + int it0 = nse_get_logT_index(tlog) - 1; + it0 = std::clamp(it0, 1, nse_table_size::ntemp-3); + + int ic0 = nse_get_ye_index(yet) - 1; + ic0 = std::clamp(ic0, 1, nse_table_size::nye-3); + + // note: this is returning the derivative wrt log10(T), so we need to + // convert to d/dT + + Real ddatadT = std::log(10.0_rt) * temp * tricubic_deriv(ir0, it0, ic0, rholog, tlog, yet, data); + + return ddatadT; + +} + #endif diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index d7d779ec26..108f638352 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -142,4 +142,10 @@ void nse_cell_c() } + std::cout << "testing temperature derivatives of cubic" << std::endl; + + std::cout << "first manual derivative of dAbar/dt" << std::endl; + + + } From ca7d025229c766b233d45934fc9b996abde82a4b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 11 Dec 2023 19:48:41 -0500 Subject: [PATCH 2/4] add a method to compute the T derivative of an NSE table quantity this works by differentiating the cubic interpolant --- nse_tabular/nse_table.H | 2 +- .../test_nse_interp/ci-benchmarks/aprox19.out | 14 ++++++++--- unit_test/test_nse_interp/nse_cell.H | 25 ++++++++++++++++++- 3 files changed, 36 insertions(+), 5 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 3f098a39de..fdf6d582fb 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -465,7 +465,7 @@ nse_interp_dT(const Real temp, const Real rho, const Real ye, const T& data) { // note: this is returning the derivative wrt log10(T), so we need to // convert to d/dT - Real ddatadT = std::log(10.0_rt) * temp * tricubic_deriv(ir0, it0, ic0, rholog, tlog, yet, data); + Real ddatadT = tricubic_dT(ir0, it0, ic0, rholog, tlog, yet, data) / (std::log(10.0_rt) * temp); return ddatadT; diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out index a881ca2f75..95a5ad1edc 100644 --- a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out +++ b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out @@ -1,5 +1,5 @@ -Initializing AMReX (23.11-5-gd36463103dae)... -AMReX (23.11-5-gd36463103dae) initialized +Initializing AMReX (23.12-8-g43d71da32fa4)... +AMReX (23.12-8-g43d71da32fa4) initialized starting the single zone burn... reading the NSE table (C++) ... rho, T, Ye = 1230000000 5180000000 0.472 @@ -57,4 +57,12 @@ X(Fe54) = 0.9104458693 X(Ni56) = 0.0003104223608 X(n) = 2.004738472e-08 X(p) = 9.620335872e-06 -AMReX (23.11-5-gd36463103dae) finalized + +testing temperature derivatives of cubic +first finite-difference derivatives +dAbar/dT = -1.070778042e-09 +dbea/dT = -6.886297725e-12 +now using derivative of the interpolant +dAbar/dT = -1.070778122e-09 +dbea/dT = -1.070778122e-09 +AMReX (23.12-8-g43d71da32fa4) finalized diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index 8708d1e161..9286ac27de 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -139,10 +139,33 @@ void nse_cell_c() } + std::cout << std::endl; std::cout << "testing temperature derivatives of cubic" << std::endl; - std::cout << "first manual derivative of dAbar/dt" << std::endl; + std::cout << "first finite-difference derivatives" << std::endl; + + Real abar_old = nse_state.abar; + Real bea_old = nse_state.bea; + Real T_old = nse_state.T; + + const 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; + + Real dabardT = nse_interp_dT(temperature, density, ye, + nse_table::abartab); + + Real dbeadT = nse_interp_dT(temperature, density, ye, + nse_table::beatab); + std::cout << "dAbar/dT = " << dabardT << std::endl; + std::cout << "dbea/dT = " << dabardT << std::endl; } From a7b63f2157fce658c7c60060ace193a0ed6c8162 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 11 Dec 2023 19:52:44 -0500 Subject: [PATCH 3/4] fix output --- unit_test/test_nse_interp/ci-benchmarks/aprox19.out | 2 +- unit_test/test_nse_interp/nse_cell.H | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out index 95a5ad1edc..42a1b41381 100644 --- a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out +++ b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out @@ -64,5 +64,5 @@ dAbar/dT = -1.070778042e-09 dbea/dT = -6.886297725e-12 now using derivative of the interpolant dAbar/dT = -1.070778122e-09 -dbea/dT = -1.070778122e-09 +dbea/dT = -6.886272826e-12 AMReX (23.12-8-g43d71da32fa4) finalized diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index 9286ac27de..1ab91c2c87 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -165,7 +165,7 @@ void nse_cell_c() nse_table::beatab); std::cout << "dAbar/dT = " << dabardT << std::endl; - std::cout << "dbea/dT = " << dabardT << std::endl; + std::cout << "dbea/dT = " << dbeadT << std::endl; } From 2ee64555110968900bab92693ac2df99a3b8069f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 17 Dec 2023 12:02:36 -0500 Subject: [PATCH 4/4] switch from std::pow() to amrex::Math::powi --- nse_tabular/nse_table.H | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index fdf6d582fb..360465aa15 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -133,12 +133,13 @@ Real cubic(const Real* xs, const Real* fs, const Real dx, const Real x) { // to the data (xs, fs) // we take x_i to be x[1] - Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * std::pow(dx, 3)); + Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * amrex::Math::powi<3>(dx)); Real b = (-2 * fs[1] + fs[2] + fs[0]) / (2 * dx * dx); Real c = (-3 * fs[1] + 6 * fs[2] - fs[3] - 2 * fs[0]) / (6 * dx); Real d = fs[1]; - return a * std::pow(x - xs[1], 3) + b * std::pow(x - xs[1], 2) + c * (x - xs[1]) + d; + return a * amrex::Math::powi<3>(x - xs[1]) + + b * amrex::Math::powi<2>(x - xs[1]) + c * (x - xs[1]) + d; } @@ -156,12 +157,12 @@ Real cubic_deriv(const Real* xs, const Real* fs, const Real dx, const Real x) { // we take x_i to be x[1] // then return dfdx = 3 a (x - x_i)**2 + 2 b (x - x_i) + c - Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * std::pow(dx, 3)); + Real a = (3 * fs[1] - 3 * fs[2] + fs[3] - fs[0]) / (6 * amrex::Math::powi<3>(dx)); Real b = (-2 * fs[1] + fs[2] + fs[0]) / (2 * dx * dx); Real c = (-3 * fs[1] + 6 * fs[2] - fs[3] - 2 * fs[0]) / (6 * dx); //Real d = fs[1]; - return 3.0_rt * a * std::pow(x - xs[1], 2) + 2.0_rt * b * (x - xs[1]) + c; + return 3.0_rt * a * amrex::Math::powi<2>(x - xs[1]) + 2.0_rt * b * (x - xs[1]) + c; }