From 7de212422501a88c1c1f536b93217721e3dc6701 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 20 Dec 2023 14:29:03 -0500 Subject: [PATCH] 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 | 192 ++++++++++++++++++ .../test_nse_interp/ci-benchmarks/aprox19.out | 16 +- unit_test/test_nse_interp/main.cpp | 2 +- unit_test/test_nse_interp/nse_cell.H | 159 ++++++++++++++- 4 files changed, 363 insertions(+), 6 deletions(-) create mode 100644 nse_tabular/nse_eos.H diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H new file mode 100644 index 0000000000..54c29af304 --- /dev/null +++ b/nse_tabular/nse_eos.H @@ -0,0 +1,192 @@ +#ifndef NSE_EOS_H +#define NSE_EOS_H + +#include + +#include + +#include +#include +#include + + +/// +/// 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 diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out index 42a1b41381..655488226f 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.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 diff --git a/unit_test/test_nse_interp/main.cpp b/unit_test/test_nse_interp/main.cpp index 880fc83255..df7e627507 100644 --- a/unit_test/test_nse_interp/main.cpp +++ b/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(); diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index 1ab91c2c87..db923ea117 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -8,12 +8,12 @@ #include #include - -using namespace unit_test_rp; +#include void nse_cell_c() { + using namespace unit_test_rp; std::cout << "rho, T, Ye = " << density << " " << temperature << " " << ye << std::endl; @@ -167,5 +167,160 @@ void nse_cell_c() std::cout << "dAbar/dT = " << dabardT << std::endl; std::cout << "dbea/dT = " << dbeadT << 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 = temperature; + nse_state.rho = density; + nse_state.Ye = ye; + + nse_interp(nse_state); + + Real abar_orig = nse_state.abar; + + eos_t eos_state; + + eos_state.T = temperature; + eos_state.rho = density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; + + eos(eos_input_rt, eos_state); + + 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 + + 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 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 = temperature; + eos_state.e = e_new; + eos_state.rho = density; + eos_state.aux[iye] = ye; + eos_state.aux[iabar] = abar_orig; + + 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 e consistency check (new method): " << eos_state.e << " " << e_new << std::endl; + std::cout << std::endl; + } + + + { + + // now redo it for pressure + + nse_state.T = temperature; + nse_state.rho = density; + nse_state.Ye = ye; + + nse_interp(nse_state); + + Real abar_orig = nse_state.abar; + + eos_t eos_state; + + eos_state.T = temperature; + eos_state.rho = density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; + + eos(eos_input_rt, eos_state); + + 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 + + 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 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 = temperature; + eos_state.p = p_new; + eos_state.rho = density; + eos_state.aux[iye] = ye; + eos_state.aux[iabar] = abar_orig; + + 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 p consistency check (new method): " << eos_state.p << " " << p_new << std::endl; + } }