add NSE EOS interfaces for finding a consistent density (#1430)

these mirror the existing interfaces for temperature

nse_tabular/nse_eos.H

@@ -10,15 +10,16 @@
 #include <nse_table.H>
 
 
+///
+/// This function inverts this form of the EOS to find the T
+/// that satisfies the EOS and NSE given an input e and rho.
 ///
 /// if we are in NSE, then the entire thermodynamic state is just
 /// a function of rho, T, Ye.  We can write the energy as:
 ///
 ///    e = e(rho, T, Y_e, Abar(rho, T, Ye))
 ///
 /// where we note that Abar is a function of those same inputs.
-/// This function inverts this form of the EOS to find the T
-/// that satisfies the EOS and NSE.
 ///
 /// The basic idea is that Abar and Zbar are both functions of
 /// rho, T, Ye through the NSE table, so we express the energy
@@ -100,15 +101,107 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
 }
 
 
+///
+/// This function inverts this form of the EOS to find the rho
+/// that satisfies the EOS and NSE given an input e and T.
+///
+/// if we are in NSE, then the entire thermodynamic state is just
+/// a function of rho, T, Ye.  We can write the energy as:
+///
+///    e = e(rho, T, Y_e, Abar(rho, T, Ye))
+///
+/// where we note that Abar is a function of those same inputs.
+///
+/// The basic idea is that Abar and Zbar are both functions of
+/// rho, T, Ye through the NSE table, so we express the energy
+/// as:
+///
+///      e = e(rho, T, Abar(rho, T, Ye), Zbar(rho, T, Ye)
+///
+/// and NR on that.  Note that Zbar = Ye Abar, so we can group
+/// those derivative terms together.
+///
+/// rho and abar come in as initial guesses and are updated
+/// on output
+///
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void
+nse_rho_abar_from_e(const Real T, const Real e_in, const Real Ye,
+                    Real& rho, Real& abar) {
+
+    using namespace amrex;
+    using namespace AuxZero;
+
+    const Real dtol{1.e-6_rt};
+    const int max_iter{100};
+
+    // we need the full EOS type, since we need de/dA
+    //eos_extra_t eos_state;
+
+    bool converged{false};
+
+    int iter{};
+
+    nse_table_t nse_state;
+
+    while (not converged && iter < max_iter) {
+
+        // call NSE table to get Abar
+        nse_state.T = T;
+        nse_state.rho = rho;
+        nse_state.Ye = Ye;
+
+        constexpr bool skip_X_fill{true};
+        nse_interp(nse_state, skip_X_fill);
+        Real abar_old = nse_state.abar;
+
+        // call the EOS with the initial guess for rho
+        eos_extra_t eos_state;
+        eos_state.rho = rho;
+        eos_state.T = T;
+        eos_state.aux[iye] = Ye;
+        eos_state.aux[iabar] = abar_old;
+        eos(eos_input_rt, eos_state);
+
+        // f is the quantity we want to zero
+
+        Real f = eos_state.e - e_in;
+
+        Real dabar_drho = nse_interp_drho(T, rho, Ye, nse_table::abartab);
+
+        // compute the correction to our guess
+
+        Real drho = -f / (eos_state.dedr + eos_state.dedA * dabar_drho
+                                    + Ye * eos_state.dedZ * dabar_drho);
+
+        // update the density
+        rho = std::clamp(rho + drho, 0.25 * rho, 4.0 * rho);
+
+        // check convergence
+
+        if (std::abs(drho) < dtol * rho) {
+            converged = true;
+        }
+        iter++;
+    }
+
+    // rho is set to the last rho
+    // we just need to save abar for output
+    abar = nse_state.abar;
+
+}
+
+
+///
+/// This function inverts this form of the EOS to find the T
+/// that satisfies the EOS and NSE given an input p and rho.
 ///
 /// if we are in NSE, then the entire thermodynamic state is just
 /// a function of rho, T, Ye.  We can write the pressure as:
 ///
 ///    p = [(rho, T, Y_e, Abar(rho, T, Ye))
 ///
 /// where we note that Abar is a function of those same inputs.
-/// This function inverts this form of the EOS to find the T
-/// that satisfies the EOS and NSE.
 ///
 /// The basic idea is that Abar and Zbar are both functions of
 /// rho, T, Ye through the NSE table, so we express the pressure
@@ -189,4 +282,95 @@ nse_T_abar_from_p(const Real rho, const Real p_in, const Real Ye,
 
 }
 
+
+///
+/// This function inverts this form of the EOS to find the rho
+/// that satisfies the EOS and NSE given an input p and T.
+///
+/// if we are in NSE, then the entire thermodynamic state is just
+/// a function of rho, T, Ye.  We can write the pressure as:
+///
+///    p = [(rho, T, Y_e, Abar(rho, T, Ye))
+///
+/// where we note that Abar is a function of those same inputs.
+///
+/// The basic idea is that Abar and Zbar are both functions of
+/// rho, T, Ye through the NSE table, so we express the pressure
+/// as:
+///
+///      p = p(rho, T, Abar(rho, T, Ye), Zbar(rho, T, Ye)
+///
+/// and NR on that.  Note that Zbar = Ye Abar, so we can group
+/// those derivative terms together.
+///
+/// rho and abar come in as initial guesses and are updated
+/// on output
+///
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void
+nse_rho_abar_from_p(const Real T, const Real p_in, const Real Ye,
+                    Real& rho, Real& abar) {
+
+    using namespace amrex;
+    using namespace AuxZero;
+
+    const Real dtol{1.e-6_rt};
+    const int max_iter{100};
+
+    // we need the full EOS type, since we need de/dA
+    //eos_extra_t eos_state;
+
+    bool converged{false};
+
+    int iter{};
+
+    nse_table_t nse_state;
+
+    while (not converged && iter < max_iter) {
+
+        // call NSE table to get Abar
+        nse_state.T = T;
+        nse_state.rho = rho;
+        nse_state.Ye = Ye;
+
+        constexpr bool skip_X_fill{true};
+        nse_interp(nse_state, skip_X_fill);
+        Real abar_old = nse_state.abar;
+
+        // call the EOS with the initial guess for rho
+        eos_extra_t eos_state;
+        eos_state.rho = rho;
+        eos_state.T = T;
+        eos_state.aux[iye] = Ye;
+        eos_state.aux[iabar] = abar_old;
+        eos(eos_input_rt, eos_state);
+
+        // f is the quantity we want to zero
+
+        Real f = eos_state.p - p_in;
+
+        Real dabar_drho = nse_interp_drho(T, rho, Ye, nse_table::abartab);
+
+        // compute the correction to our guess
+
+        Real drho = -f / (eos_state.dpdr + eos_state.dpdA * dabar_drho
+                                    + Ye * eos_state.dpdZ * dabar_drho);
+
+        // update the density
+        rho = std::clamp(rho + drho, 0.25 * rho, 4.0 * rho);
+
+        // check convergence
+
+        if (std::abs(drho) < dtol * rho) {
+            converged = true;
+        }
+        iter++;
+    }
+
+    // rho is set to the last rho
+    // we just need to save abar for output
+    abar = nse_state.abar;
+
+}
+
 #endif

nse_tabular/nse_table.H

@@ -327,6 +327,69 @@ Real tricubic_dT(const int ir0, const int it0, const int ic0,
 
 }
 
+
+///
+/// take the density derivative of a table quantity by differentiating
+/// the cubic interpolant
+///
+template <typename T>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+Real tricubic_drho(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 T interpolations (one in each rho plane)
+
+    Real d2[4];
+
+    for (int ii = 0; ii < 4; ++ii) {
+
+        Real _d[] = {d1[ii][0], d1[ii][1], d1[ii][2], d1[ii][3]};
+        d2[ii] = cubic(Ts, _d, nse_table_size::dlogT, temp);
+    }
+
+    // finally do the remaining interpolation over rho, but return
+    // the derivative of the interpolant
+
+    Real val = cubic_deriv(rhos, d2, nse_table_size::dlogrho, rho);
+
+    return val;
+
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) {
 
@@ -470,6 +533,59 @@ nse_interp_dT(const Real temp, const Real rho, const Real ye, const T& data) {
 
     return ddatadT;
 
+
+}
+
+
+///
+/// compute the density 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_drho(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(rho), so we need to
+    // convert to d/drho
+
+    Real ddatadrho = tricubic_drho(ir0, it0, ic0, rholog, tlog, yet, data) / (std::log(10.0_rt) * rho);
+
+    return ddatadrho;
+
 }
 
 #endif

unit_test/test_nse_interp/ci-benchmarks/aprox19.out

@@ -66,13 +66,28 @@ now using derivative of the interpolant
 dAbar/dT = -1.072562604e-09
 dbea/dT = -6.867048589e-12
 
-EOS e consistency check (old method): 1.395278886e+18 1.38844906e+18
+
+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 e consistency check (new method): 1.38844906e+18 1.38844906e+18
+EOS T from e consistency check (new method): 1.38844906e+18 1.38844906e+18
 
-EOS p consistency check (old method): 6.622159603e+26 6.577850616e+26
+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 p consistency check (new method): 6.577850616e+26 6.577850616e+26
+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