diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H index 9bee82b25b..e484f923df 100644 --- a/nse_tabular/nse_eos.H +++ b/nse_tabular/nse_eos.H @@ -20,7 +20,16 @@ /// This function inverts this form of the EOS to find the T /// that satisfies the EOS and NSE. /// -/// Note: T and abar come in as initial guesses and are updated +/// 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 @@ -31,7 +40,7 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye, using namespace amrex; using namespace AuxZero; - const Real ttol{1.e-8_rt}; + const Real ttol{1.e-6_rt}; const int max_iter{100}; const Real eps{1.e-8_rt}; @@ -67,18 +76,15 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye, Real f = eos_state.e - e_in; - // perturb T and call NSE again to get perturbed Abar - nse_state.T *= (1.0_rt + eps); - nse_interp(nse_state, skip_X_fill); - - // estimate dAbar/dT - Real dabar_dT = (nse_state.abar - abar_old) / (eps * T); + 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); + + 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.5 * T, 2.0 * T); + T = std::clamp(T + dT, 0.25 * T, 4.0 * T); // check convergence diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out index 42a1b41381..be47c3708c 100644 --- a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out +++ b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out @@ -65,4 +65,9 @@ dbea/dT = -6.886297725e-12 now using derivative of the interpolant dAbar/dT = -1.070778122e-09 dbea/dT = -6.886272826e-12 + +EOS e consistency check: 1.395273834e+18 1.388449367e+18 +updated T: 6394612134 +change in abar: 55.60621897 50.27196857 +EOS e consistency check: 1.388449367e+18 1.388449367e+18 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 246bc0f6df..39b1c60b9d 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -239,10 +239,13 @@ void nse_cell_c() 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!