diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index e2a5eded22..bee26352e2 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -17,6 +17,7 @@ #include #ifdef NSE_TABLE +#include #include #endif #ifdef NSE_NET @@ -85,18 +86,14 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // get the current NSE state from the table -- this will be used // to compute the NSE evolution sources - Real abar_out; - Real dq_out; - Real dyedt; - Real dabardt; - Real dbeadt; - Real enu; - Real X[NumSpec]; + nse_table_t nse_state; + nse_state.T = T_in; + nse_state.rho = state.rho; + nse_state.Ye = ye_inl - nse_interp(T_in, state.rho, ye_in, - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); + nse_interp(nse_state); - Real dyedt_old = dyedt; + Real dyedt_old = nse_state.dyedt; // predict the U^{n+1,*} state with only estimates of the aux @@ -107,7 +104,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // initial aux_sources Real aux_source[NumAux] = {0.0_rt}; - aux_source[iye] = rho_old * dyedt; + aux_source[iye] = rho_old * nse_state.dyedt; Real rho_aux_new[NumAux]; @@ -147,8 +144,11 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // call the NSE table using the * state to get the t^{n+1} // source estimates - nse_interp(T_new, eos_state.rho, eos_state.aux[iye], - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); + nse_state.T = T_new; + nse_state.rho = eos_state.rho; + nse_state.Ye = eos_state.aux[iye]; + + nse_interp(nse_state); // note that the abar, dq (B/A), and X here have all already // seen advection implicitly @@ -156,27 +156,27 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // compute the energy release -- we need to remove the advective part - Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * state.ydot_a[SFX+ibea]; + Real rho_bea_tilde = state.y[SRHO] * nse_state.bea - dt * state.ydot_a[SFX+ibea]; Real rho_dBEA = rho_bea_tilde - rho_bea_old; // this is MeV / nucleon * g / cm**3 // convert the energy to erg / cm**3 rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; // now get the updated neutrino term - zbar = abar_out * eos_state.aux[iye]; + zbar = nse_state.abar * eos_state.aux[iye]; #ifdef NEUTRINOS - sneut5(T_new, eos_state.rho, abar_out, zbar, + sneut5(T_new, eos_state.rho, nse_state.abar, zbar, snu, dsnudt, dsnudd, dsnuda, dsnudz); #endif rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu); // update the new state for the next pass - aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + dyedt); + aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + nse_state.dyedt); rho_aux_new[iye] = state.y[SFX+iye] + dt * state.ydot_a[SFX+iye] + dt * aux_source[iye]; - rho_aux_new[iabar] = state.y[SRHO] * abar_out; - rho_aux_new[ibea] = state.y[SRHO] * dq_out; + rho_aux_new[iabar] = state.y[SRHO] * nse_state.abar; + rho_aux_new[ibea] = state.y[SRHO] * nse_state.bea; } @@ -185,24 +185,24 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // the new mass fractions are just those that come from the table // make sure they are normalized Real sum_X{0.0_rt}; - for (auto & xn : X) { - xn = amrex::max(small_x, amrex::min(1.0_rt, xn)); + for (auto & xn : nse_state.X) { + xn = std::clamp(xn, small_x, 1.0_rt); sum_X += xn; } - for (auto & xn : X) { + for (auto & xn : nse_state.X) { xn /= sum_X; } for (int n = 0; n < NumSpec; ++n) { - state.y[SFS+n] = state.y[SRHO] * X[n]; + state.y[SFS+n] = state.y[SRHO] * nse_state.X[n]; } // aux data comes from the table (except Ye, which we predicted) state.y[SFX+iye] = eos_state.rho * eos_state.aux[iye]; - state.y[SFX+iabar] = eos_state.rho * abar_out; - state.y[SFX+ibea] = eos_state.rho * dq_out; + state.y[SFX+iabar] = eos_state.rho * nse_state.abar; + state.y[SFX+ibea] = eos_state.rho * nse_state.bea; // density and momenta have already been updated diff --git a/nse_tabular/Make.package b/nse_tabular/Make.package index 11845e7125..8f786aca2b 100644 --- a/nse_tabular/Make.package +++ b/nse_tabular/Make.package @@ -2,3 +2,4 @@ CEXE_headers += nse_table.H CEXE_headers += nse_table_check.H CEXE_headers += nse_table_data.H CEXE_sources += nse_table_data.cpp +CEXE_headers += nse_table_type.H \ No newline at end of file diff --git a/nse_tabular/nse_table.H b/nse_tabular/nse_table.H index 933b05f8dd..0a2b530da9 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -10,9 +10,13 @@ #include #include +#include + #include #include #include +#include + using namespace amrex; using namespace network_rp; @@ -235,37 +239,35 @@ Real tricubic(const int ir0, const int it0, const int ic0, } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void nse_interp(const Real T, const Real rho, const Real ye, - Real& abar, Real& bea, Real& dyedt, Real& dabardt, Real& dbeadt, Real& e_nu, - Real* X, bool skip_X_fill=false) { +void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false) { // if skip_X_fill = true then we don't fill X[] with the mass fractions. using namespace nse_table; using namespace AuxZero; - Real rholog = std::log10(rho); + Real rholog = std::log10(nse_state.rho); { Real rmin = nse_table_size::logrho_min; Real rmax = nse_table_size::logrho_max; - rholog = std::max(rmin, std::min(rmax, rholog)); + rholog = std::clamp(rholog, rmin, rmax); } - Real tlog = std::log10(T); + Real tlog = std::log10(nse_state.T); { Real tmin = nse_table_size::logT_min; Real tmax = nse_table_size::logT_max; - tlog = std::max(tmin, std::min(tmax, tlog)); + tlog = std::clamp(tlog, tmin, tmax); } - Real yet = ye; + Real yet = nse_state.Ye; { Real yemin = nse_table_size::ye_min; Real yemax = nse_table_size::ye_max; - yet = std::max(yemin, std::min(yemax, yet)); + yet = std::clamp(yet, yemin, yemax); } if (nse_table_interp_linear) { @@ -274,12 +276,12 @@ void nse_interp(const Real T, const Real rho, const Real ye, int it1 = nse_get_logT_index(tlog); int ic1 = nse_get_ye_index(yet); - abar = trilinear(ir1, it1, ic1, rholog, tlog, yet, abartab); - bea = trilinear(ir1, it1, ic1, rholog, tlog, yet, beatab); - dyedt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dyedttab); - dabardt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dabardttab); - dbeadt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dbeadttab); - e_nu = trilinear(ir1, it1, ic1, rholog, tlog, yet, enutab); + nse_state.abar = trilinear(ir1, it1, ic1, rholog, tlog, yet, abartab); + nse_state.bea = trilinear(ir1, it1, ic1, rholog, tlog, yet, beatab); + nse_state.dyedt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dyedttab); + nse_state.dabardt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dabardttab); + nse_state.dbeadt = trilinear(ir1, it1, ic1, rholog, tlog, yet, dbeadttab); + nse_state.e_nu = trilinear(ir1, it1, ic1, rholog, tlog, yet, enutab); // massfractab is 2-d, so we wrap the access in a lambda already // indexing the component @@ -288,7 +290,7 @@ void nse_interp(const Real T, const Real rho, const Real ye, for (int n = 1; n <= NumSpec; n++) { Real _X = trilinear(ir1, it1, ic1, rholog, tlog, yet, [=] (const int i) {return massfractab(n, i);}); - X[n-1] = amrex::max(0.0_rt, amrex::min(1.0_rt, _X)); + nse_state.X[n-1] = std::clamp(_X, 0.0_rt, 1.0_rt); } } @@ -308,12 +310,12 @@ void nse_interp(const Real T, const Real rho, const Real ye, int ic0 = nse_get_ye_index(yet) - 1; ic0 = amrex::max(1, amrex::min(nse_table_size::nye-3, ic0)); - abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab); - bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab); - dyedt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dyedttab); - dabardt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dabardttab); - dbeadt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dbeadttab); - e_nu = tricubic(ir0, it0, ic0, rholog, tlog, yet, enutab); + nse_state.abar = tricubic(ir0, it0, ic0, rholog, tlog, yet, abartab); + nse_state.bea = tricubic(ir0, it0, ic0, rholog, tlog, yet, beatab); + nse_state.dyedt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dyedttab); + nse_state.dabardt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dabardttab); + nse_state.dbeadt = tricubic(ir0, it0, ic0, rholog, tlog, yet, dbeadttab); + nse_state.e_nu = tricubic(ir0, it0, ic0, rholog, tlog, yet, enutab); // massfractab is 2-d, so we wrap the access in a lambda already // indexing the component @@ -322,11 +324,89 @@ void nse_interp(const Real T, const Real rho, const Real ye, for (int n = 1; n <= NumSpec; n++) { Real _X = tricubic(ir0, it0, ic0, rholog, tlog, yet, [=] (const int i) {return massfractab(n, i);}); - X[n-1] = amrex::max(0.0_rt, amrex::min(1.0_rt, _X)); + nse_state.X[n-1] = std::clamp(_X, 0.0_rt, 1.0_rt); } } } } + +/// +/// 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. +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +Real +nse_T_from_e(const Real rho, const Real e_in, const Real Ye, const Real T_guess) { + + using namespace AuxZero; + + const Real ttol{1.e-8_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; + + bool converged{false}; + + Real T{T_guess}; + int iter{}; + + while (not converged && iter < max_iter) { + + // call NSE table to get Abar + nse_table_t nse_state; + 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_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; + + // 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); + + // compute the correction to our guess + Real dT = f / (eos_state.dedT + eos_state.dedA * dabar_dT); + + // update the temperature + T = std::clamp(T + dT, 0.5 * T, 2.0 * T); + + // check convergence + + if (std::abs(dT) < ttol * T) { + converged = true; + } + iter++; + } + + return T; + +} + #endif