add a method to compute the T derivative of an NSE table quantity

nse_tabular/nse_table.H

@@ -121,6 +121,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) {
 
@@ -129,15 +133,40 @@ 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;
 
 }
 
+///
+/// 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 * 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 * amrex::Math::powi<2>(x - xs[1]) + 2.0_rt * b * (x - xs[1]) + c;
+
+}
+
+
 template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 Real trilinear(const int ir1, const int it1, const int ic1,
@@ -236,6 +265,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 <typename T>
+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(nse_table_t& nse_state, bool skip_X_fill=false) {
 
@@ -300,13 +391,13 @@ void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) {
         // 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);
 
         nse_state.abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab);
         nse_state.bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab);
@@ -329,4 +420,56 @@ void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) {
 
 }
 
+
+///
+/// compute the temperature derivative of the table quantity data
+/// at the point T, rho, ye by using cubic interpolation
+///
+template <typename T>
+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 = tricubic_dT(ir0, it0, ic0, rholog, tlog, yet, data) / (std::log(10.0_rt) * temp);
+
+    return ddatadT;
+
+}
+
 #endif

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 = -6.886272826e-12
+AMReX (23.12-8-g43d71da32fa4) finalized

unit_test/test_nse_interp/nse_cell.H

@@ -139,4 +139,33 @@ void nse_cell_c()
     }
 
 
+    std::cout << std::endl;
+    std::cout << "testing temperature derivatives of cubic" << 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 = " << dbeadT << std::endl;
+
+
 }