diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index e2a5eded22..9c8345e8d2 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_in; - 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/integration/nse_update_strang.H b/integration/nse_update_strang.H index 4bbbc5d8e4..8e594caefa 100644 --- a/integration/nse_update_strang.H +++ b/integration/nse_update_strang.H @@ -17,6 +17,7 @@ #include #ifdef NSE_TABLE +#include #include #endif #ifdef NSE_NET @@ -37,13 +38,6 @@ void nse_burn(BurnT& state, const Real dt) { using namespace AuxZero; // use the NSE table - Real abar_out; - Real dq_out; - Real dyedt; - Real dabardt; - Real dbeadt; - Real enu; - Real X[NumSpec]; // get the energy -- we are assuming that rho, T are valid on input @@ -58,25 +52,29 @@ void nse_burn(BurnT& state, const Real dt) { // source estimates. The thermodynamnics here is specified // in terms of the auxiliary composition, Ye, abar, and B/A - nse_interp(T_in, state.rho, state.aux[iye], - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); + nse_table_t nse_state; + nse_state.T = T_in; + nse_state.rho = state.rho; + nse_state.Ye = state.aux[iye]; + + nse_interp(nse_state); // update Ye - state.aux[iye] += dt * dyedt; + state.aux[iye] += dt * nse_state.dyedt; // now get the composition from the table using the updated Ye - nse_interp(T_in, state.rho, state.aux[iye], - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); + nse_state.Ye = state.aux[iye]; + nse_interp(nse_state); // this is MeV / nucleon -- here aux has not yet been updated, so we // access the old binding energy - Real deltaq = dq_out - state.aux[ibea]; + Real deltaq = nse_state.bea - state.aux[ibea]; state.aux[ibea] += deltaq; - state.aux[iabar] = abar_out; + state.aux[iabar] = nse_state.abar; #elif defined(NSE_NET) @@ -116,7 +114,7 @@ void nse_burn(BurnT& state, const Real dt) { #if defined(NSE_TABLE) for (int n = 0; n < NumSpec; ++n) { - state.xn[n] = X[n]; + state.xn[n] = nse_state.X[n]; } #elif defined(NSE_NET) for (int n = 0; n < NumSpec; ++n) { 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..7b9cfb6e99 100644 --- a/nse_tabular/nse_table.H +++ b/nse_tabular/nse_table.H @@ -13,6 +13,8 @@ #include #include #include +#include + using namespace amrex; using namespace network_rp; @@ -235,37 +237,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 +274,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 +288,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 +308,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,7 +322,7 @@ 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); } } } diff --git a/nse_tabular/nse_table_type.H b/nse_tabular/nse_table_type.H new file mode 100644 index 0000000000..82cebe863d --- /dev/null +++ b/nse_tabular/nse_table_type.H @@ -0,0 +1,23 @@ +#ifndef NSE_TABLE_TYPE_H +#define NSE_TABLE_TYPE_H + +#include +#include + +struct nse_table_t { + + amrex::Real T; + amrex::Real rho; + amrex::Real Ye; + amrex::Real abar; + amrex::Real bea; + amrex::Real dyedt; + amrex::Real dabardt; + amrex::Real dbeadt; + amrex::Real e_nu; + amrex::Real X[NumSpec]; + +}; + + +#endif diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H index d7d779ec26..e6ea2e3421 100644 --- a/unit_test/test_nse_interp/nse_cell.H +++ b/unit_test/test_nse_interp/nse_cell.H @@ -7,6 +7,7 @@ #include #include +#include using namespace unit_test_rp; @@ -120,25 +121,21 @@ void nse_cell_c() std::cout << "tricubic interpolated values: " << std::endl; - Real abar; - Real bea; - Real dyedt; - Real dabardt; - Real dbeadt; - Real e_nu; - Real X[NumSpec]; - - nse_interp(temperature, density, ye, - abar, bea, dyedt, dabardt, dbeadt, e_nu, X); - - std::cout << "abar = " << abar << std::endl; - std::cout << "bea = " << bea << std::endl; - std::cout << "dyedt = " << dyedt << std::endl; - std::cout << "dabardt = " << dabardt << std::endl; - std::cout << "dbeadt = " << dbeadt << std::endl; - std::cout << "e_nu = " << e_nu << std::endl; + nse_table_t nse_state; + nse_state.T = temperature; + nse_state.rho = density; + nse_state.Ye = ye; + + nse_interp(nse_state); + + std::cout << "abar = " << nse_state.abar << std::endl; + std::cout << "bea = " << nse_state.bea << std::endl; + std::cout << "dyedt = " << nse_state.dyedt << std::endl; + std::cout << "dabardt = " << nse_state.dabardt << std::endl; + std::cout << "dbeadt = " << nse_state.dbeadt << std::endl; + std::cout << "e_nu = " << nse_state.e_nu << std::endl; for (int n = 0; n < NumSpec; ++n) { - std::cout << "X(" << short_spec_names_cxx[n] << ") = " << X[n] << std::endl; + std::cout << "X(" << short_spec_names_cxx[n] << ") = " << nse_state.X[n] << std::endl; }