From da880f02194405e1f5093f215a3244428b5a70e1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 18 Dec 2023 11:57:20 -0500 Subject: [PATCH] add a pressure interface --- nse_tabular/nse_eos.H | 91 ++++++++++- .../test_nse_interp/ci-benchmarks/aprox19.out | 15 +- unit_test/test_nse_interp/nse_cell.H | 146 +++++++++++++----- 3 files changed, 208 insertions(+), 44 deletions(-) diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H index e484f923df..54c29af304 100644 --- a/nse_tabular/nse_eos.H +++ b/nse_tabular/nse_eos.H @@ -42,7 +42,6 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye, const Real ttol{1.e-6_rt}; const int max_iter{100}; - const Real eps{1.e-8_rt}; // we need the full EOS type, since we need de/dA //eos_extra_t eos_state; @@ -100,4 +99,94 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye, } + +/// +/// 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 be47c3708c..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 @@ -66,8 +66,13 @@ 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 +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: 1.388449367e+18 1.388449367e+18 -AMReX (23.12-8-g43d71da32fa4) finalized +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/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index 39b1c60b9d..07ea37c4e1 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -193,64 +193,134 @@ void nse_cell_c() // // 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; + // first get the abar consistent with our inputs and find e - nse_interp(nse_state); + nse_state.T = temperature; + nse_state.rho = density; + nse_state.Ye = ye; - Real abar_orig = nse_state.abar; + nse_interp(nse_state); - eos_t eos_state; + Real abar_orig = nse_state.abar; - eos_state.T = temperature; - eos_state.rho = density; - eos_state.aux[iye] = nse_state.Ye; - eos_state.aux[iabar] = nse_state.abar; + eos_t eos_state; - eos(eos_input_rt, eos_state); + eos_state.T = temperature; + eos_state.rho = density; + eos_state.aux[iye] = nse_state.Ye; + eos_state.aux[iabar] = nse_state.abar; - Real e_orig = eos_state.e; + eos(eos_input_rt, eos_state); - // now perturn e and find T, then redo NSE to get abar and finally - // see if we get back our new e + Real e_orig = eos_state.e; - Real e_new = eos_state.e * 1.05; + // now perturn e and find T, then redo NSE to get abar and finally + // see if we get back our new e - eos_state.e = e_new; - eos(eos_input_re, eos_state); + Real e_new = eos_state.e * 1.05; - nse_state.T = eos_state.T; - nse_interp(nse_state); + 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; - eos_state.aux[iabar] = nse_state.abar; - eos(eos_input_rt, eos_state); + // attempt 2: + // now we try the new interface. This effectively does: + // e', rho, Ye -> Abar', T' - std::cout << "EOS e consistency check: " << eos_state.e << " " << e_new << std::endl; + eos_state.T = temperature; + eos_state.e = e_new; + eos_state.rho = density; + eos_state.aux[iye] = ye; + eos_state.aux[iabar] = abar_orig; - // now we try the new interface. This effectively does: - // e', rho, Ye -> Abar', T' + Real abar_start = eos_state.aux[iabar]; - eos_state.T = temperature; - eos_state.e = e_new; - eos_state.rho = density; - eos_state.aux[iye] = ye; - eos_state.aux[iabar] = abar_orig; + 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; + } + + + { - Real abar_start = eos_state.aux[iabar]; + // now redo it for pressure - nse_T_abar_from_e(eos_state.rho, eos_state.e, eos_state.aux[iye], - eos_state.T, eos_state.aux[iabar]); + nse_state.T = temperature; + nse_state.rho = density; + nse_state.Ye = ye; - std::cout << "updated T: " << eos_state.T << std::endl; - std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl; + nse_interp(nse_state); - // now check if we get back the correct e! + Real abar_orig = nse_state.abar; - eos(eos_input_rt, eos_state); + eos_t eos_state; - std::cout << "EOS e consistency check: " << eos_state.e << " " << e_new << std::endl; + 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 perturn 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; + } }