an EOS interface with NSE (#1405)

solving for T(e) or T(p) in NSE is trickier than our standard interface because Abar and Zbar are
also a function of rho, T, Ye through the NSE table. This does a self-consistent inversion of the
EOS together with the NSE table.

nse_tabular/nse_eos.H

@@ -0,0 +1,192 @@
+#ifndef NSE_EOS_H
+#define NSE_EOS_H
+
+#include <AMReX_REAL.H>
+
+#include <eos.H>
+
+#include <extern_parameters.H>
+#include <nse_table_type.H>
+#include <nse_table.H>
+
+
+///
+/// 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
+/// 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.
+///
+/// T and abar come in as initial guesses and are updated
+/// on output
+///
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void
+nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
+                  Real& T, Real& abar) {
+
+    using namespace amrex;
+    using namespace AuxZero;
+
+    const Real ttol{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 T
+        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_dT = nse_interp_dT(T, rho, Ye, nse_table::abartab);
+
+        // compute the correction to our guess
+
+        Real dT = -f / (eos_state.dedT + eos_state.dedA * dabar_dT
+                                  + Ye * eos_state.dedZ * dabar_dT);
+
+        // update the temperature
+        T = std::clamp(T + dT, 0.25 * T, 4.0 * T);
+
+        // check convergence
+
+        if (std::abs(dT) < ttol * T) {
+            converged = true;
+        }
+        iter++;
+    }
+
+    // T is set to the last T
+    // we just need to save abar for output
+    abar = nse_state.abar;
+
+}
+
+
+///
+/// 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
+/// 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.
+///
+/// T and abar come in as initial guesses and are updated
+/// on output
+///
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void
+nse_T_abar_from_p(const Real rho, const Real p_in, const Real Ye,
+                  Real& T, Real& abar) {
+
+    using namespace amrex;
+    using namespace AuxZero;
+
+    const Real ttol{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 T
+        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_dT = nse_interp_dT(T, rho, Ye, nse_table::abartab);
+
+        // compute the correction to our guess
+
+        Real dT = -f / (eos_state.dpdT + eos_state.dpdA * dabar_dT
+                                  + Ye * eos_state.dpdZ * dabar_dT);
+
+        // update the temperature
+        T = std::clamp(T + dT, 0.25 * T, 4.0 * T);
+
+        // check convergence
+
+        if (std::abs(dT) < ttol * T) {
+            converged = true;
+        }
+        iter++;
+    }
+
+    // T is set to the last T
+    // we just need to save abar for output
+    abar = nse_state.abar;
+
+}
+
+#endif

unit_test/test_nse_interp/ci-benchmarks/aprox19.out

@@ -1,5 +1,5 @@
-Initializing AMReX (23.12-8-g43d71da32fa4)...
-AMReX (23.12-8-g43d71da32fa4) initialized
+Initializing AMReX (23.11-5-gd36463103dae)...
+AMReX (23.11-5-gd36463103dae) initialized
 starting the single zone burn...
 reading the NSE table (C++) ...
 rho, T, Ye = 1230000000 5180000000 0.472
@@ -65,4 +65,14 @@ 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
+
+EOS e consistency check (old method): 1.395273834e+18 1.388449367e+18
+updated T: 6394612134
+change in abar: 55.60621897 50.27196857
+EOS e consistency check (new method): 1.388449367e+18 1.388449367e+18
+
+EOS p consistency check (old method): 6.622132093e+26 6.57785215e+26
+updated T: 6466848772
+change in abar: 55.60621897 49.50670173
+EOS p consistency check (new method): 6.57785215e+26 6.57785215e+26
+AMReX (23.11-5-gd36463103dae) finalized

unit_test/test_nse_interp/main.cpp

@@ -23,7 +23,7 @@ int main(int argc, char *argv[]) {
   init_unit_test();
 
   // C++ EOS initialization (must be done after Fortran eos_init and init_extern_parameters)
-  eos_init(small_temp, small_dens);
+  eos_init(unit_test_rp::small_temp, unit_test_rp::small_dens);
 
   // C++ Network, RHS, screening, rates initialization
   network_init();