From ca6d25613a12dc8c7741ec58e770b49242aabf69 Mon Sep 17 00:00:00 2001 From: Max Katz Date: Sun, 27 Nov 2022 17:18:01 -0500 Subject: [PATCH] Construct some containers for data in rhs.H (#1054) Containers for state data and rate data are constructed. These significantly clean up interfaces (and therefore make them more extensible in the future) and also help quiet unused variable warnings. --- interfaces/Make.package | 1 + interfaces/rhs_type.H | 71 ++++++ {rates => interfaces}/tfactors.H | 2 + networks/aprox13/actual_network.H | 120 +++++----- networks/aprox19/actual_network.H | 322 +++++++++++++------------ networks/aprox21/actual_network.H | 358 ++++++++++++++-------------- networks/iso7/actual_network.H | 58 +++-- networks/rhs.H | 375 ++++++++++++------------------ rates/Make.package | 1 - 9 files changed, 634 insertions(+), 674 deletions(-) rename {rates => interfaces}/tfactors.H (99%) diff --git a/interfaces/Make.package b/interfaces/Make.package index 53fab28ddc..a4f813be8c 100644 --- a/interfaces/Make.package +++ b/interfaces/Make.package @@ -9,6 +9,7 @@ CEXE_sources += eos_data.cpp CEXE_headers += network.H CEXE_headers += rhs_type.H +CEXE_headers += tfactors.H CEXE_sources += network_initialization.cpp ifeq ($(USE_CONDUCTIVITY), TRUE) diff --git a/interfaces/rhs_type.H b/interfaces/rhs_type.H index d1f595c9f1..2e9e4d95b7 100644 --- a/interfaces/rhs_type.H +++ b/interfaces/rhs_type.H @@ -2,6 +2,11 @@ #define rhs_type_H #include +#include +#ifdef SCREENING +#include +#endif +#include namespace RHS { @@ -43,6 +48,72 @@ struct rhs_t { int additional_reaction_3; }; +constexpr amrex::Real tab_tlo = 6.0e0_rt; +constexpr amrex::Real tab_thi = 10.0e0_rt; +constexpr int tab_per_decade = 2000; +constexpr int nrattab = static_cast(tab_thi - tab_tlo) * tab_per_decade + 1; +constexpr int tab_imax = static_cast(tab_thi - tab_tlo) * tab_per_decade + 1; +constexpr amrex::Real tab_tstp = (tab_thi - tab_tlo) / static_cast(tab_imax - 1); + +extern AMREX_GPU_MANAGED Array1D ttab; + +// A struct that contains the terms needed to evaluate a tabulated rate. +struct rate_tab_t +{ + amrex::Real alfa, beta, gama, delt; + int iat; + + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void initialize (amrex::Real temp) { + // hash locate + constexpr int mp = 4; + iat = static_cast((std::log10(temp) - tab_tlo) / tab_tstp) + 1; + iat = amrex::max(1, amrex::min(iat - 1, tab_imax - mp + 1)); + + // setup the lagrange interpolation coefficients for a cubic + amrex::Real x = temp; + amrex::Real x1 = ttab(iat); + amrex::Real x2 = ttab(iat+1); + amrex::Real x3 = ttab(iat+2); + amrex::Real x4 = ttab(iat+3); + amrex::Real a = x - x1; + amrex::Real b = x - x2; + amrex::Real c = x - x3; + amrex::Real d = x - x4; + amrex::Real e = x1 - x2; + amrex::Real f = x1 - x3; + amrex::Real g = x1 - x4; + amrex::Real h = x2 - x3; + amrex::Real p = x2 - x4; + amrex::Real q = x3 - x4; + alfa = b * c * d / (e * f * g); + beta = -a * c * d / (e * h * p); + gama = a * b * d / (f * h * q); + delt = -a * b * c / (g * p * q); + } +}; + +struct rhs_state_t +{ + amrex::Real rho; + tf_t tf; + rate_tab_t tab; +#ifdef SCREENING + plasma_state_t pstate; +#endif + amrex::Real y_e; + amrex::Real eta; + Array1D y; +}; + +struct rate_t +{ + amrex::Real fr; + amrex::Real rr; + amrex::Real frdt; + amrex::Real rrdt; +}; + } // namespace RHS #endif diff --git a/rates/tfactors.H b/interfaces/tfactors.H similarity index 99% rename from rates/tfactors.H rename to interfaces/tfactors.H index f2b2b1be45..1221490490 100644 --- a/rates/tfactors.H +++ b/interfaces/tfactors.H @@ -4,6 +4,8 @@ #include #include +using namespace amrex; + struct tf_t { amrex::Real temp; amrex::Real t9; diff --git a/networks/aprox13/actual_network.H b/networks/aprox13/actual_network.H index 6e89a97723..b359dd97ac 100644 --- a/networks/aprox13/actual_network.H +++ b/networks/aprox13/actual_network.H @@ -712,125 +712,119 @@ namespace RHS { template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void evaluate_analytical_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt) + void evaluate_analytical_rate (const rhs_state_t& state, rate_t& rates) { using namespace Species; using namespace Rates; if constexpr (rate == C12_He4_to_O16) { - rate_c12ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He4_He4_He4_to_C12) { - rate_triplealf(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_triplealf(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_C12_to_Ne20_He4) { - rate_c12c12(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12c12(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Mg24_He4) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Si28) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_Si28_He4) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_P31_P) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_S32) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_He4_to_Ne20) { - rate_o16ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ne20_He4_to_Mg24) { - rate_ne20ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ne20ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Si28) { - rate_mg24ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_S32) { - rate_si28ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Ar36) { - rate_s32ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_Ca40) { - rate_ar36ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Ti44) { - rate_ca40ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_Cr48) { - rate_ti44ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Fe52) { - rate_cr48ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Ni56) { - rate_fe52ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Al27_P) { - rate_mg24ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Al27_P_to_Si28) { - rate_al27pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_al27pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_P31_P) { - rate_si28ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == P31_P_to_S32) { - rate_p31pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_p31pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Cl35_P) { - rate_s32ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cl35_P_to_Ar36) { - rate_cl35pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cl35pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_K39_P) { - rate_ar36ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == K39_P_to_Ca40) { - rate_k39pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_k39pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Sc43_P) { - rate_ca40ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Sc43_P_to_Ti44) { - rate_sc43pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_sc43pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_V47_P) { - rate_ti44ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == V47_P_to_Cr48) { - rate_v47pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_v47pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Mn51_P) { - rate_cr48ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mn51_P_to_Fe52) { - rate_mn51pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mn51pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Co55_P) { - rate_fe52ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Co55_P_to_Ni56) { - rate_co55pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_co55pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } } template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void postprocess_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt, - Real& fr1, Real& fr1dt, Real& rr1, Real& rr1dt, - Real& fr2, Real& fr2dt, Real& rr2, Real& rr2dt, - Real& fr3, Real& fr3dt, Real& rr3, Real& rr3dt) + void postprocess_rate ([[maybe_unused]] const rhs_state_t& state, rate_t& rates, + rate_t& rates1, rate_t& rates2, rate_t& rates3) { using namespace Species; using namespace Rates; @@ -856,20 +850,20 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr += fr1 * (1.0_rt - rate_ir); - frdt += fr1dt * (1.0_rt - rate_ir) - fr1 * dratedt_ir; + rates.fr += rates1.fr * (1.0_rt - rate_ir); + rates.frdt += rates1.frdt * (1.0_rt - rate_ir) - rates1.fr * dratedt_ir; - rr += rr2 * rate_ir; - rrdt += rr2dt * rate_ir + rr2 * dratedt_ir; + rates.rr += rates2.rr * rate_ir; + rates.rrdt += rates2.rrdt * rate_ir + rates2.rr * dratedt_ir; } // The O16+O16->Si28+He4 reaction has an additional contribution from the O16+O16->P31+P @@ -878,17 +872,17 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr += fr3 * rate_ir; - frdt += fr3dt * rate_ir + fr3 * dratedt_ir; + rates.fr += rates3.fr * rate_ir; + rates.frdt += rates3.frdt * rate_ir + rates3.fr * dratedt_ir; } // The O16+O16->S32 reaction has an additional contribution from the O16+O16->P31+P @@ -897,17 +891,17 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr += fr3 * (1.0_rt - rate_ir); - frdt += fr3dt * (1.0_rt - rate_ir) - fr3 * dratedt_ir; + rates.fr += rates3.fr * (1.0_rt - rate_ir); + rates.frdt += rates3.frdt * (1.0_rt - rate_ir) - rates3.fr * dratedt_ir; } } diff --git a/networks/aprox19/actual_network.H b/networks/aprox19/actual_network.H index 6c955ade1c..58e36686e7 100644 --- a/networks/aprox19/actual_network.H +++ b/networks/aprox19/actual_network.H @@ -1263,164 +1263,162 @@ namespace RHS { template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void evaluate_analytical_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt) + void evaluate_analytical_rate (const rhs_state_t& state, rate_t& rates) { using namespace Species; using namespace Rates; if constexpr (rate == H1_H1_to_He3) { Real tmp1, tmp2; - rate_pp(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_pp(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == P_to_N) { - if (tf.temp >= 1.0e6_rt && rho >= 1.0e-9_rt) { + if (state.tf.temp >= 1.0e6_rt && state.rho >= 1.0e-9_rt) { Real tmp1, tmp2; - ecapnuc(eta, tf.temp, fr, rr, tmp1, tmp2); + ecapnuc(state.eta, state.tf.temp, rates.fr, rates.rr, tmp1, tmp2); } } else if constexpr (rate == P_N_to_H2) { - rate_png(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_png(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == H2_P_to_He3) { - rate_dpg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_dpg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He3_N_to_He4) { - rate_he3ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_he3ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He3_He3_to_He4_2H1) { Real tmp1, tmp2; - rate_he3he3(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_he3he3(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == He3_He4_H1_to_2He4) { Real tmp1, tmp2; - rate_he3he4(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_he3he4(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == C12_He4_to_O16) { - rate_c12ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He4_He4_He4_to_C12) { - rate_triplealf(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_triplealf(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_H1_to_N13) { - rate_c12pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_C12_to_Ne20_He4) { - rate_c12c12(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12c12(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Mg24_He4) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N14_H1_to_O15) { - rate_n14pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n14pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N15_H1_to_O16) { - rate_n15pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n15pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N15_H1_to_C12_He4) { - rate_n15pa(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n15pa(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N14_He4_to_F18) { - rate_n14ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n14ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_H1_to_F17) { - rate_o16pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_P31_P) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_He4_to_Ne20) { - rate_o16ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ne20_He4_to_Mg24) { - rate_ne20ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ne20ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Si28) { - rate_mg24ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_S32) { - rate_si28ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Ar36) { - rate_s32ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_Ca40) { - rate_ar36ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Ti44) { - rate_ca40ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_Cr48) { - rate_ti44ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Fe52) { - rate_cr48ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Ni56) { - rate_fe52ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Al27_P) { - rate_mg24ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Al27_P_to_Si28) { - rate_al27pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_al27pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_P31_P) { - rate_si28ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == P31_P_to_S32) { - rate_p31pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_p31pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Cl35_P) { - rate_s32ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cl35_P_to_Ar36) { - rate_cl35pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cl35pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_K39_P) { - rate_ar36ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == K39_P_to_Ca40) { - rate_k39pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_k39pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Sc43_P) { - rate_ca40ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Sc43_P_to_Ti44) { - rate_sc43pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_sc43pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_V47_P) { - rate_ti44ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == V47_P_to_Cr48) { - rate_v47pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_v47pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Mn51_P) { - rate_cr48ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mn51_P_to_Fe52) { - rate_mn51pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mn51pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Co55_P) { - rate_fe52ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Co55_P_to_Ni56) { - rate_co55pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_co55pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_N_to_Fe53) { - rate_fe52ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe53_N_to_Fe54) { - rate_fe53ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe53ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe54_P_to_Co55) { - rate_fe54pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe54pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ni56_to_Fe54) { - if (tf.temp >= 1.0e6_rt && rho >= 1.0e-9_rt) { + if (state.tf.temp >= 1.0e6_rt && state.rho >= 1.0e-9_rt) { Real tmp; - langanke(tf.temp, rho, y(Ni56), y_e, fr, tmp); + langanke(state.tf.temp, state.rho, state.y(Ni56), state.y_e, rates.fr, tmp); // We scaled the number of nucleons per reaction // in defining the reaction so that the number of reactants @@ -1428,19 +1426,15 @@ namespace RHS { // so to compensate, we scale the rate by 1/54. This serves // the equivalent role to the c54 constant in the original aprox19. - fr *= 1.0_rt / 54.0_rt; + rates.fr *= 1.0_rt / 54.0_rt; } } } template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void postprocess_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt, - Real& fr1, Real& fr1dt, Real& rr1, Real& rr1dt, - Real& fr2, Real& fr2dt, Real& rr2, Real& rr2dt, - Real& fr3, Real& fr3dt, Real& rr3, Real& rr3dt) + void postprocess_rate (const rhs_state_t& state, rate_t& rates, + rate_t& rates1, rate_t& rates2, rate_t& rates3) { using namespace Species; using namespace Rates; @@ -1452,45 +1446,45 @@ namespace RHS { // Rate 2 == H2_P_to_He3 == irdpg (forward), irhegp (reverse) // Rate 3 == He3_N_to_He4 == irheng (forward), irhegn (reverse) - Real denom = rr2 * rr1 + y(N) * fr3 * rr1 + y(N) * y(P) * fr3 * fr2; + Real denom = rates2.rr * rates1.rr + state.y(N) * rates3.fr * rates1.rr + state.y(N) * state.y(P) * rates3.fr * rates2.fr; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr2dt * rr1 + rr2 * rr1dt + - y(N) * (fr3dt * rr1 + fr3 * rr1dt) + - y(N) * y(P) * (fr3dt * fr2 + fr3 * fr2dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates2.rrdt * rates1.rr + rates2.rr * rates1.rrdt + + state.y(N) * (rates3.frdt * rates1.rr + rates3.fr * rates1.rrdt) + + state.y(N) * state.y(P) * (rates3.frdt * rates2.fr + rates3.fr * rates2.frdt); Real zz = 1.0_rt / denom; // iralf1 in the original aprox19 - rr = rr3 * rr2 * rr1 * zz; - rrdt = rr3dt * rr2 * rr1 * zz + - rr3 * rr2dt * rr1 * zz + - rr3 * rr2 * rr1dt * zz - - rr * zz * denomdt; + rates.rr = rates3.rr * rates2.rr * rates1.rr * zz; + rates.rrdt = rates3.rrdt * rates2.rr * rates1.rr * zz + + rates3.rr * rates2.rrdt * rates1.rr * zz + + rates3.rr * rates2.rr * rates1.rrdt * zz - + rates.rr * zz * denomdt; // iralf2 in the original aprox19 - fr = fr3 * fr2 * fr1 * zz; - frdt = fr3dt * fr2 * fr1 * zz + - fr3 * fr2dt * fr1 * zz + - fr3 * fr2 * fr1dt * zz - - fr * zz * denomdt; + rates.fr = rates3.fr * rates2.fr * rates1.fr * zz; + rates.frdt = rates3.frdt * rates2.fr * rates1.fr * zz + + rates3.fr * rates2.frdt * rates1.fr * zz + + rates3.fr * rates2.fr * rates1.frdt * zz - + rates.fr * zz * denomdt; } } // Set this equal to the p+e->n rate. if constexpr (rate == H1_H1_H1_to_He3) { - fr = fr1; + rates.fr = rates1.fr; } // Beta limit He3 + He4 by the B8 decay half-life if constexpr (rate == He3_He4_H1_to_2He4) { - Real xx = 0.896_rt / y(He4); - fr = amrex::min(fr, xx); - if (fr == xx) { - frdt = 0.0_rt; + Real xx = 0.896_rt / state.y(He4); + rates.fr = amrex::min(rates.fr, xx); + if (rates.fr == xx) { + rates.frdt = 0.0_rt; } } @@ -1514,20 +1508,20 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr += fr1 * (1.0_rt - rate_ir); - frdt += fr1dt * (1.0_rt - rate_ir) - fr1 * dratedt_ir; + rates.fr += rates1.fr * (1.0_rt - rate_ir); + rates.frdt += rates1.frdt * (1.0_rt - rate_ir) - rates1.fr * dratedt_ir; - rr += rr2 * rate_ir; - rrdt += rr2dt * rate_ir + rr2 * dratedt_ir; + rates.rr += rates2.rr * rate_ir; + rates.rrdt += rates2.rrdt * rate_ir + rates2.rr * dratedt_ir; } // The O16+O16->Si28+He4 reaction has an additional contribution from the O16+O16->P31+P @@ -1536,17 +1530,17 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr = 0.56_rt * fr3 + 0.34_rt * fr3 * rate_ir; - frdt = 0.56_rt * fr3dt + 0.34_rt * fr3dt * rate_ir + 0.34_rt * fr3 * dratedt_ir; + rates.fr = 0.56_rt * rates3.fr + 0.34_rt * rates3.fr * rate_ir; + rates.frdt = 0.56_rt * rates3.frdt + 0.34_rt * rates3.frdt * rate_ir + 0.34_rt * rates3.fr * dratedt_ir; } // The O16+O16->S32 reaction has an additional contribution from the O16+O16->P31+P @@ -1555,31 +1549,31 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr = 0.1_rt * fr3 + 0.34_rt * fr3 * (1.0_rt - rate_ir); - frdt = 0.1_rt * fr3dt + 0.34_rt * fr3dt * (1.0_rt - rate_ir) - 0.34_rt * fr3 * dratedt_ir; + rates.fr = 0.1_rt * rates3.fr + 0.34_rt * rates3.fr * (1.0_rt - rate_ir); + rates.frdt = 0.1_rt * rates3.frdt + 0.34_rt * rates3.frdt * (1.0_rt - rate_ir) - 0.34_rt * rates3.fr * dratedt_ir; } if constexpr (rate == C12_2H1_to_N14) { - fr = fr1; - frdt = fr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; } if constexpr (rate == C12_O16_to_Si28) { - fr = fr1; - frdt = fr1dt; - rr = rr1; - rrdt = rr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; + rates.rr = rates1.rr; + rates.rrdt = rates1.rrdt; } if constexpr (rate == N14_2H1_to_O16) @@ -1588,25 +1582,25 @@ namespace RHS { // Beta limit - if (y(H1) > 1.0e-30_rt) { - Real xx = 5.68e-3_rt / (y(H1) * 1.57_rt); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + if (state.y(H1) > 1.0e-30_rt) { + Real xx = 5.68e-3_rt / (state.y(H1) * 1.57_rt); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } } // Rate 2: O15_H1_to_O16 // Rate 3: O15_H1_to_C12_He4 - Real tot = fr2 + fr3; + Real tot = rates2.fr + rates3.fr; Real invtot = 1.0_rt / tot; - fr = fr1 * (fr2 * invtot); + rates.fr = rates1.fr * (rates2.fr * invtot); - Real dtotdt = fr2dt + fr3dt; + Real dtotdt = rates2.frdt + rates3.frdt; - frdt = fr1 * (fr2dt * invtot - fr2 * invtot * invtot * dtotdt); + rates.frdt = rates1.fr * (rates2.frdt * invtot - rates2.fr * invtot * invtot * dtotdt); } if constexpr (rate == N14_2H1_to_C12_He4) @@ -1615,25 +1609,25 @@ namespace RHS { // Beta limit - if (y(H1) > 1.0e-30_rt) { - Real xx = 5.68e-3_rt / (y(H1) * 1.57_rt); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + if (state.y(H1) > 1.0e-30_rt) { + Real xx = 5.68e-3_rt / (state.y(H1) * 1.57_rt); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } } // Rate 2: O15_H1_to_O16 // Rate 3: O15_H1_to_C12_He4 - Real tot = fr2 + fr3; + Real tot = rates2.fr + rates3.fr; Real invtot = 1.0_rt / tot; - fr = fr1 * (fr3 * invtot); + rates.fr = rates1.fr * (rates3.fr * invtot); - Real dtotdt = fr2dt + fr3dt; + Real dtotdt = rates2.frdt + rates3.frdt; - frdt = fr1 * (fr3dt * invtot - fr3 * invtot * invtot * dtotdt); + rates.frdt = rates1.fr * (rates3.frdt * invtot - rates3.fr * invtot * invtot * dtotdt); } if constexpr (rate == N14_He4_to_Ne20) @@ -1643,22 +1637,22 @@ namespace RHS { // on each side would be integers, so to compensate, // we'll halve the rate. - fr = fr1 * 0.5_rt; - frdt = fr1dt * 0.5_rt; + rates.fr = rates1.fr * 0.5_rt; + rates.frdt = rates1.frdt * 0.5_rt; } if constexpr (rate == O16_2H1_to_N14_He4) { // Beta limit - Real xx = 0.0105_rt / y(H1); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + Real xx = 0.0105_rt / state.y(H1); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } - fr = fr1; - frdt = fr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; } if constexpr (rate == Fe52_2N_to_Fe54) @@ -1667,19 +1661,19 @@ namespace RHS { // Rate 1 == Fe52_N_to_Fe53 == ir52ng (forward), ir53gn (reverse) // Rate 2 == Fe53_N_to_Fe54 == ir53ng (forward), ir54gn (reverse) - Real denom = rr1 + y(N) * fr2; - Real denomdt = rr1dt + y(N) * fr2dt; + Real denom = rates1.rr + state.y(N) * rates2.fr; + Real denomdt = rates1.rrdt + state.y(N) * rates2.frdt; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { Real zz = 1.0_rt / denom; // ir2f54 in the original aprox19 - fr = fr1 * fr2 * zz; - frdt = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.fr * zz; + rates.frdt = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - rates.fr * zz * denomdt; // ir1f54 in the original aprox19 - rr = rr2 * rr1 * zz; - rrdt = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.rr * rates1.rr * zz; + rates.rrdt = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1690,20 +1684,20 @@ namespace RHS { // Rate 2 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) // Rate 3 == Fe52_He4_to_Co55_P == irfeap (forward), ircopa (reverse) - Real denom = rr1 + y(P) * (fr2 + rr3); + Real denom = rates1.rr + state.y(P) * (rates2.fr + rates3.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr1dt + y(P) * (fr2dt + rr3dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates1.rrdt + state.y(P) * (rates2.frdt + rates3.rrdt); Real zz = 1.0_rt / denom; // ir3f54 in the original aprox19 - fr = fr1 * fr2 * zz; - frdt = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.fr * zz; + rates.frdt = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - rates.fr * zz * denomdt; // ir4f54 in the original aprox19 - rr = rr2 * rr1 * zz; - rrdt = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.rr * rates1.rr * zz; + rates.rrdt = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1714,20 +1708,20 @@ namespace RHS { // Rate 2 == Fe54_P_to_Co55 == irfepg (forward), ircogp (reverse) // Rate 3 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) - Real denom = rr2 + y(P) * (fr3 + rr1); + Real denom = rates2.rr + state.y(P) * (rates3.fr + rates1.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr2dt + y(P) * (fr3dt + rr1dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates2.rrdt + state.y(P) * (rates3.frdt + rates1.rrdt); Real zz = 1.0_rt / denom; // ir6f54 in the original aprox19 - fr = fr1 * rr2 * zz; - frdt = fr1dt * rr2 * zz + fr1 * rr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.rr * zz; + rates.frdt = rates1.frdt * rates2.rr * zz + rates1.fr * rates2.rrdt * zz - rates.fr * zz * denomdt; // ir5f54 in the original aprox19 - rr = fr2 * rr1 * zz; - rrdt = fr2dt * rr1 * zz + fr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.fr * rates1.rr * zz; + rates.rrdt = rates2.frdt * rates1.rr * zz + rates2.fr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1738,10 +1732,10 @@ namespace RHS { // Rate 2 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) // Rate 3 == Fe54_P_to_Co55 == irfepg (forward), ircogp (reverse) - Real denom = rr3 + y(P) * (fr2 + rr1); + Real denom = rates3.rr + state.y(P) * (rates2.fr + rates1.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr3dt + y(P) * (fr2dt + rr1dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates3.rrdt + state.y(P) * (rates2.frdt + rates1.rrdt); Real zz = 1.0_rt / denom; @@ -1750,16 +1744,16 @@ namespace RHS { // direct alpha capture sequence Fe52(a,g)Ni56. // ir7f54 in the original aprox19 - Real fr_add = fr1 * fr2 * zz; - Real frdt_add = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr_add * zz * denomdt; - fr += y(P) * fr_add; - frdt += y(P) * frdt_add; + Real fr_add = rates1.fr * rates2.fr * zz; + Real frdt_add = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - fr_add * zz * denomdt; + rates.fr += state.y(P) * fr_add; + rates.frdt += state.y(P) * frdt_add; // ir8f54 in the original aprox19 - Real rr_add = rr2 * rr1 * zz; - Real rrdt_add = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr_add * zz * denomdt; - rr += y(P) * rr_add; - rrdt += y(P) * rrdt_add; + Real rr_add = rates2.rr * rates1.rr * zz; + Real rrdt_add = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rr_add * zz * denomdt; + rates.rr += state.y(P) * rr_add; + rates.rrdt += state.y(P) * rrdt_add; } } } diff --git a/networks/aprox21/actual_network.H b/networks/aprox21/actual_network.H index 39a2d88610..62cee10173 100644 --- a/networks/aprox21/actual_network.H +++ b/networks/aprox21/actual_network.H @@ -1396,188 +1396,182 @@ namespace RHS { template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void evaluate_analytical_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt) + void evaluate_analytical_rate (const rhs_state_t& state, rate_t& rates) { using namespace Species; using namespace Rates; if constexpr (rate == H1_H1_to_He3) { Real tmp1, tmp2; - rate_pp(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_pp(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == P_to_N) { - if (tf.temp >= 1.0e6_rt && rho >= 1.0e-9_rt) { + if (state.tf.temp >= 1.0e6_rt && state.rho >= 1.0e-9_rt) { Real tmp1, tmp2; - ecapnuc(eta, tf.temp, fr, rr, tmp1, tmp2); + ecapnuc(state.eta, state.tf.temp, rates.fr, rates.rr, tmp1, tmp2); } } else if constexpr (rate == P_N_to_H2) { - rate_png(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_png(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == H2_P_to_He3) { - rate_dpg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_dpg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He3_N_to_He4) { - rate_he3ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_he3ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He3_He3_to_He4_2H1) { Real tmp1, tmp2; - rate_he3he3(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_he3he3(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == He3_He4_H1_to_2He4) { Real tmp1, tmp2; - rate_he3he4(tf, 1.0_rt, fr, frdt, tmp1, tmp2); + rate_he3he4(state.tf, 1.0_rt, rates.fr, rates.frdt, tmp1, tmp2); } else if constexpr (rate == C12_He4_to_O16) { - rate_c12ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He4_He4_He4_to_C12) { - rate_triplealf(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_triplealf(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_H1_to_N13) { - rate_c12pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_C12_to_Ne20_He4) { - rate_c12c12(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12c12(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Mg24_He4) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N14_H1_to_O15) { - rate_n14pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n14pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N15_H1_to_O16) { - rate_n15pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n15pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N15_H1_to_C12_He4) { - rate_n15pa(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n15pa(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == N14_He4_to_F18) { - rate_n14ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_n14ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_H1_to_F17) { - rate_o16pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_P31_P) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_He4_to_Ne20) { - rate_o16ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ne20_He4_to_Mg24) { - rate_ne20ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ne20ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Si28) { - rate_mg24ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_S32) { - rate_si28ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Ar36) { - rate_s32ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_Ca40) { - rate_ar36ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Ti44) { - rate_ca40ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_Cr48) { - rate_ti44ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Fe52) { - rate_cr48ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Ni56) { - rate_fe52ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Al27_P) { - rate_mg24ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Al27_P_to_Si28) { - rate_al27pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_al27pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Si28_He4_to_P31_P) { - rate_si28ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_si28ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == P31_P_to_S32) { - rate_p31pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_p31pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == S32_He4_to_Cl35_P) { - rate_s32ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_s32ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cl35_P_to_Ar36) { - rate_cl35pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cl35pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ar36_He4_to_K39_P) { - rate_ar36ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ar36ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == K39_P_to_Ca40) { - rate_k39pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_k39pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Sc43_P) { - rate_ca40ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Sc43_P_to_Ti44) { - rate_sc43pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_sc43pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ti44_He4_to_V47_P) { - rate_ti44ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ti44ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == V47_P_to_Cr48) { - rate_v47pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_v47pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Cr48_He4_to_Mn51_P) { - rate_cr48ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_cr48ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mn51_P_to_Fe52) { - rate_mn51pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mn51pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_He4_to_Co55_P) { - rate_fe52ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Co55_P_to_Ni56) { - rate_co55pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_co55pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe52_N_to_Fe53) { - rate_fe52ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe52ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe53_N_to_Fe54) { - rate_fe53ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe53ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe54_P_to_Co55) { - rate_fe54pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe54pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ni56_to_Fe56) { - if (tf.temp >= 1.0e6_rt && rho >= 1.0e-9_rt) { + if (state.tf.temp >= 1.0e6_rt && state.rho >= 1.0e-9_rt) { Real tmp; - langanke(tf.temp, rho, y(Ni56), y_e, fr, tmp); + langanke(state.tf.temp, state.rho, state.y(Ni56), state.y_e, rates.fr, tmp); } } else if constexpr (rate == Fe54_N_to_Fe55) { - rate_fe54ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe54ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe55_N_to_Fe56) { - rate_fe55ng(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe55ng(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe54_He4_to_Co57_P) { - rate_fe54ap(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe54ap(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Fe56_P_to_Co57) { - rate_fe56pg(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_fe56pg(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } } template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void postprocess_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt, - Real& fr1, Real& fr1dt, Real& rr1, Real& rr1dt, - Real& fr2, Real& fr2dt, Real& rr2, Real& rr2dt, - Real& fr3, Real& fr3dt, Real& rr3, Real& rr3dt) + void postprocess_rate (const rhs_state_t& state, rate_t& rates, + rate_t& rates1, rate_t& rates2, rate_t& rates3) { using namespace Species; using namespace Rates; @@ -1589,45 +1583,45 @@ namespace RHS { // Rate 2 == H2_P_to_He3 == irdpg (forward), irhegp (reverse) // Rate 3 == He3_N_to_He4 == irheng (forward), irhegn (reverse) - Real denom = rr2 * rr1 + y(N) * fr3 * rr1 + y(N) * y(P) * fr3 * fr2; + Real denom = rates2.rr * rates1.rr + state.y(N) * rates3.fr * rates1.rr + state.y(N) * state.y(P) * rates3.fr * rates2.fr; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr2dt * rr1 + rr2 * rr1dt + - y(N) * (fr3dt * rr1 + fr3 * rr1dt) + - y(N) * y(P) * (fr3dt * fr2 + fr3 * fr2dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates2.rrdt * rates1.rr + rates2.rr * rates1.rrdt + + state.y(N) * (rates3.frdt * rates1.rr + rates3.fr * rates1.rrdt) + + state.y(N) * state.y(P) * (rates3.frdt * rates2.fr + rates3.fr * rates2.frdt); Real zz = 1.0_rt / denom; // iralf1 in the original aprox21 - rr = rr3 * rr2 * rr1 * zz; - rrdt = rr3dt * rr2 * rr1 * zz + - rr3 * rr2dt * rr1 * zz + - rr3 * rr2 * rr1dt * zz - - rr * zz * denomdt; + rates.rr = rates3.rr * rates2.rr * rates1.rr * zz; + rates.rrdt = rates3.rrdt * rates2.rr * rates1.rr * zz + + rates3.rr * rates2.rrdt * rates1.rr * zz + + rates3.rr * rates2.rr * rates1.rrdt * zz - + rates.rr * zz * denomdt; // iralf2 in the original aprox21 - fr = fr3 * fr2 * fr1 * zz; - frdt = fr3dt * fr2 * fr1 * zz + - fr3 * fr2dt * fr1 * zz + - fr3 * fr2 * fr1dt * zz - - fr * zz * denomdt; + rates.fr = rates3.fr * rates2.fr * rates1.fr * zz; + rates.frdt = rates3.frdt * rates2.fr * rates1.fr * zz + + rates3.fr * rates2.frdt * rates1.fr * zz + + rates3.fr * rates2.fr * rates1.frdt * zz - + rates.fr * zz * denomdt; } } // Set this equal to the p+e->n rate. if constexpr (rate == H1_H1_H1_to_He3) { - fr = fr1; + rates.fr = rates1.fr; } // Beta limit He3 + He4 by the B8 decay half-life if constexpr (rate == He3_He4_H1_to_2He4) { - Real xx = 0.896_rt / y(He4); - fr = amrex::min(fr, xx); - if (fr == xx) { - frdt = 0.0_rt; + Real xx = 0.896_rt / state.y(He4); + rates.fr = amrex::min(rates.fr, xx); + if (rates.fr == xx) { + rates.frdt = 0.0_rt; } } @@ -1651,20 +1645,20 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr += fr1 * (1.0_rt - rate_ir); - frdt += fr1dt * (1.0_rt - rate_ir) - fr1 * dratedt_ir; + rates.fr += rates1.fr * (1.0_rt - rate_ir); + rates.frdt += rates1.frdt * (1.0_rt - rate_ir) - rates1.fr * dratedt_ir; - rr += rr2 * rate_ir; - rrdt += rr2dt * rate_ir + rr2 * dratedt_ir; + rates.rr += rates2.rr * rate_ir; + rates.rrdt += rates2.rrdt * rate_ir + rates2.rr * dratedt_ir; } // The O16+O16->Si28+He4 reaction has an additional contribution from the O16+O16->P31+P @@ -1673,17 +1667,17 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr = 0.56_rt * fr3 + 0.34_rt * fr3 * rate_ir; - frdt = 0.56_rt * fr3dt + 0.34_rt * fr3dt * rate_ir + 0.34_rt * fr3 * dratedt_ir; + rates.fr = 0.56_rt * rates3.fr + 0.34_rt * rates3.fr * rate_ir; + rates.frdt = 0.56_rt * rates3.frdt + 0.34_rt * rates3.frdt * rate_ir + 0.34_rt * rates3.fr * dratedt_ir; } // The O16+O16->S32 reaction has an additional contribution from the O16+O16->P31+P @@ -1692,31 +1686,31 @@ namespace RHS { { Real rate_ir = 0.0e0_rt; Real dratedt_ir = 0.0e0_rt; - Real denom = rr1 + fr2; - Real denomdt = rr1dt + fr2dt; + Real denom = rates1.rr + rates2.fr; + Real denomdt = rates1.rrdt + rates2.frdt; if (denom > 1.0e-30_rt) { Real zz = 1.0_rt / denom; - rate_ir = rr1 * zz; - dratedt_ir = (rr1dt - rate_ir * denomdt) * zz; + rate_ir = rates1.rr * zz; + dratedt_ir = (rates1.rrdt - rate_ir * denomdt) * zz; } - fr = 0.1_rt * fr3 + 0.34_rt * fr3 * (1.0_rt - rate_ir); - frdt = 0.1_rt * fr3dt + 0.34_rt * fr3dt * (1.0_rt - rate_ir) - 0.34_rt * fr3 * dratedt_ir; + rates.fr = 0.1_rt * rates3.fr + 0.34_rt * rates3.fr * (1.0_rt - rate_ir); + rates.frdt = 0.1_rt * rates3.frdt + 0.34_rt * rates3.frdt * (1.0_rt - rate_ir) - 0.34_rt * rates3.fr * dratedt_ir; } if constexpr (rate == C12_2H1_to_N14) { - fr = fr1; - frdt = fr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; } if constexpr (rate == C12_O16_to_Si28) { - fr = fr1; - frdt = fr1dt; - rr = rr1; - rrdt = rr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; + rates.rr = rates1.rr; + rates.rrdt = rates1.rrdt; } if constexpr (rate == N14_2H1_to_O16) @@ -1725,25 +1719,25 @@ namespace RHS { // Beta limit - if (y(H1) > 1.0e-30_rt) { - Real xx = 5.68e-3_rt / (y(H1) * 1.57_rt); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + if (state.y(H1) > 1.0e-30_rt) { + Real xx = 5.68e-3_rt / (state.y(H1) * 1.57_rt); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } } // Rate 2: O15_H1_to_O16 // Rate 3: O15_H1_to_C12_He4 - Real tot = fr2 + fr3; + Real tot = rates2.fr + rates3.fr; Real invtot = 1.0_rt / tot; - fr = fr1 * (fr2 * invtot); + rates.fr = rates1.fr * (rates2.fr * invtot); - Real dtotdt = fr2dt + fr3dt; + Real dtotdt = rates2.frdt + rates3.frdt; - frdt = fr1 * (fr2dt * invtot - fr2 * invtot * invtot * dtotdt); + rates.frdt = rates1.fr * (rates2.frdt * invtot - rates2.fr * invtot * invtot * dtotdt); } if constexpr (rate == N14_2H1_to_C12_He4) @@ -1752,25 +1746,25 @@ namespace RHS { // Beta limit - if (y(H1) > 1.0e-30_rt) { - Real xx = 5.68e-3_rt / (y(H1) * 1.57_rt); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + if (state.y(H1) > 1.0e-30_rt) { + Real xx = 5.68e-3_rt / (state.y(H1) * 1.57_rt); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } } // Rate 2: O15_H1_to_O16 // Rate 3: O15_H1_to_C12_He4 - Real tot = fr2 + fr3; + Real tot = rates2.fr + rates3.fr; Real invtot = 1.0_rt / tot; - fr = fr1 * (fr3 * invtot); + rates.fr = rates1.fr * (rates3.fr * invtot); - Real dtotdt = fr2dt + fr3dt; + Real dtotdt = rates2.frdt + rates3.frdt; - frdt = fr1 * (fr3dt * invtot - fr3 * invtot * invtot * dtotdt); + rates.frdt = rates1.fr * (rates3.frdt * invtot - rates3.fr * invtot * invtot * dtotdt); } if constexpr (rate == N14_He4_to_Ne20) @@ -1780,22 +1774,22 @@ namespace RHS { // on each side would be integers, so to compensate, // we'll halve the rate. - fr = fr1 * 0.5_rt; - frdt = fr1dt * 0.5_rt; + rates.fr = rates1.fr * 0.5_rt; + rates.frdt = rates1.frdt * 0.5_rt; } if constexpr (rate == O16_2H1_to_N14_He4) { // Beta limit - Real xx = 0.0105_rt / y(H1); - fr1 = amrex::min(fr1, xx); - if (fr1 == xx) { - fr1dt = 0.0_rt; + Real xx = 0.0105_rt / state.y(H1); + rates1.fr = amrex::min(rates1.fr, xx); + if (rates1.fr == xx) { + rates1.frdt = 0.0_rt; } - fr = fr1; - frdt = fr1dt; + rates.fr = rates1.fr; + rates.frdt = rates1.frdt; } if constexpr (rate == Fe52_2N_to_Fe54) @@ -1804,19 +1798,19 @@ namespace RHS { // Rate 1 == Fe52_N_to_Fe53 == ir52ng (forward), ir53gn (reverse) // Rate 2 == Fe53_N_to_Fe54 == ir53ng (forward), ir54gn (reverse) - Real denom = rr1 + y(N) * fr2; - Real denomdt = rr1dt + y(N) * fr2dt; + Real denom = rates1.rr + state.y(N) * rates2.fr; + Real denomdt = rates1.rrdt + state.y(N) * rates2.frdt; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { Real zz = 1.0_rt / denom; // ir2f54 in the original aprox21 - fr = fr1 * fr2 * zz; - frdt = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.fr * zz; + rates.frdt = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - rates.fr * zz * denomdt; // ir1f54 in the original aprox21 - rr = rr2 * rr1 * zz; - rrdt = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.rr * rates1.rr * zz; + rates.rrdt = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1827,20 +1821,20 @@ namespace RHS { // Rate 2 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) // Rate 3 == Fe52_He4_to_Co55_P == irfeap (forward), ircopa (reverse) - Real denom = rr1 + y(P) * (fr2 + rr3); + Real denom = rates1.rr + state.y(P) * (rates2.fr + rates3.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr1dt + y(P) * (fr2dt + rr3dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates1.rrdt + state.y(P) * (rates2.frdt + rates3.rrdt); Real zz = 1.0_rt / denom; // ir3f54 in the original aprox21 - fr = fr1 * fr2 * zz; - frdt = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.fr * zz; + rates.frdt = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - rates.fr * zz * denomdt; // ir4f54 in the original aprox21 - rr = rr2 * rr1 * zz; - rrdt = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.rr * rates1.rr * zz; + rates.rrdt = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1851,20 +1845,20 @@ namespace RHS { // Rate 2 == Fe54_P_to_Co55 == irfepg (forward), ircogp (reverse) // Rate 3 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) - Real denom = rr2 + y(P) * (fr3 + rr1); + Real denom = rates2.rr + state.y(P) * (rates3.fr + rates1.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr2dt + y(P) * (fr3dt + rr1dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates2.rrdt + state.y(P) * (rates3.frdt + rates1.rrdt); Real zz = 1.0_rt / denom; // ir6f54 in the original aprox21 - fr = fr1 * rr2 * zz; - frdt = fr1dt * rr2 * zz + fr1 * rr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.rr * zz; + rates.frdt = rates1.frdt * rates2.rr * zz + rates1.fr * rates2.rrdt * zz - rates.fr * zz * denomdt; // ir5f54 in the original aprox21 - rr = fr2 * rr1 * zz; - rrdt = fr2dt * rr1 * zz + fr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.fr * rates1.rr * zz; + rates.rrdt = rates2.frdt * rates1.rr * zz + rates2.fr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1875,10 +1869,10 @@ namespace RHS { // Rate 2 == Co55_P_to_Ni56 == ircopg (forward), irnigp (reverse) // Rate 3 == Fe54_P_to_Co55 == irfepg (forward), ircogp (reverse) - Real denom = rr3 + y(P) * (fr2 + rr1); + Real denom = rates3.rr + state.y(P) * (rates2.fr + rates1.rr); - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { - Real denomdt = rr3dt + y(P) * (fr2dt + rr1dt); + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { + Real denomdt = rates3.rrdt + state.y(P) * (rates2.frdt + rates1.rrdt); Real zz = 1.0_rt / denom; @@ -1887,23 +1881,23 @@ namespace RHS { // direct alpha capture sequence Fe52(a,g)Ni56. // ir7f54 in the original aprox19 - Real fr_add = fr1 * fr2 * zz; - Real frdt_add = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr_add * zz * denomdt; - fr += y(P) * fr_add; - frdt += y(P) * frdt_add; + Real fr_add = rates1.fr * rates2.fr * zz; + Real frdt_add = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - fr_add * zz * denomdt; + rates.fr += state.y(P) * fr_add; + rates.frdt += state.y(P) * frdt_add; // ir8f54 in the original aprox19 - Real rr_add = rr2 * rr1 * zz; - Real rrdt_add = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr_add * zz * denomdt; - rr += y(P) * rr_add; - rrdt += y(P) * rrdt_add; + Real rr_add = rates2.rr * rates1.rr * zz; + Real rrdt_add = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rr_add * zz * denomdt; + rates.rr += state.y(P) * rr_add; + rates.rrdt += state.y(P) * rrdt_add; } } if constexpr (rate == Fe56_to_Cr56) { // Rate 1 == Ni56_to_Fe56 - fr = 1.0e-4 * fr1; + rates.fr = 1.0e-4 * rates1.fr; } if constexpr (rate == Fe54_2N_to_Fe56) @@ -1912,19 +1906,19 @@ namespace RHS { // Rate 1 == Fe54_N_to_Fe55 == ir54ng (forward), ir55gn (reverse) // Rate 2 == Fe55_N_to_Fe56 == ir55ng (forward), ir56gn (reverse) - Real denom = rr1 + y(N) * fr2; - Real denomdt = rr1dt + y(N) * fr2dt; + Real denom = rates1.rr + state.y(N) * rates2.fr; + Real denomdt = rates1.rrdt + state.y(N) * rates2.frdt; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { Real zz = 1.0_rt / denom; // irfe56_aux2 in the original aprox21 - fr = fr1 * fr2 * zz; - frdt = fr1dt * fr2 * zz + fr1 * fr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.fr * zz; + rates.frdt = rates1.frdt * rates2.fr * zz + rates1.fr * rates2.frdt * zz - rates.fr * zz * denomdt; // irfe56_aux1 in the original aprox21 - rr = rr2 * rr1 * zz; - rrdt = rr2dt * rr1 * zz + rr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.rr * rates1.rr * zz; + rates.rrdt = rates2.rrdt * rates1.rr * zz + rates2.rr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } @@ -1934,19 +1928,19 @@ namespace RHS { // Rate 1 == Fe54_He4_to_Co57_P == ir54ap (forward), irco57pa (reverse) // Rate 2 == Fe56_P_to_Co57 == irfe56pg (forward), irco57gp (reverse) - Real denom = rr2 + y(P) * rr1; - Real denomdt = rr2dt + y(P) * rr1dt; + Real denom = rates2.rr + state.y(P) * rates1.rr; + Real denomdt = rates2.rrdt + state.y(P) * rates1.rrdt; - if (denom > 1.0e-50_rt && tf.t9 > 1.5_rt) { + if (denom > 1.0e-50_rt && state.tf.t9 > 1.5_rt) { Real zz = 1.0_rt / denom; // irfe56_aux4 in the original aprox21 - fr = fr1 * rr2 * zz; - frdt = fr1dt * rr2 * zz + fr1 * rr2dt * zz - fr * zz * denomdt; + rates.fr = rates1.fr * rates2.rr * zz; + rates.frdt = rates1.frdt * rates2.rr * zz + rates1.fr * rates2.rrdt * zz - rates.fr * zz * denomdt; // irfe56_aux3 in the original aprox21 - rr = fr2 * rr1 * zz; - rrdt = fr2dt * rr1 * zz + fr2 * rr1dt * zz - rr * zz * denomdt; + rates.rr = rates2.fr * rates1.rr * zz; + rates.rrdt = rates2.frdt * rates1.rr * zz + rates2.fr * rates1.rrdt * zz - rates.rr * zz * denomdt; } } } diff --git a/networks/iso7/actual_network.H b/networks/iso7/actual_network.H index 14fc223105..8454c8439c 100644 --- a/networks/iso7/actual_network.H +++ b/networks/iso7/actual_network.H @@ -260,53 +260,47 @@ namespace RHS { template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void evaluate_analytical_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt) + void evaluate_analytical_rate (const rhs_state_t& state, rate_t& rates) { using namespace Species; using namespace Rates; if constexpr (rate == C12_He4_to_O16) { - rate_c12ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == He4_He4_He4_to_C12) { - rate_triplealf(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_triplealf(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_C12_to_Ne20_He4) { - rate_c12c12(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12c12(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Mg24_He4) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == C12_O16_to_Si28) { - rate_c12o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_c12o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_O16_to_Si28_He4) { - rate_o16o16(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16o16(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == O16_He4_to_Ne20) { - rate_o16ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_o16ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ne20_He4_to_Mg24) { - rate_ne20ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ne20ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Mg24_He4_to_Si28) { - rate_mg24ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_mg24ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } else if constexpr (rate == Ca40_He4_to_Ti44) { - rate_ca40ag(tf, 1.0_rt, fr, frdt, rr, rrdt); + rate_ca40ag(state.tf, 1.0_rt, rates.fr, rates.frdt, rates.rr, rates.rrdt); } } template AMREX_GPU_HOST_DEVICE AMREX_INLINE - void postprocess_rate (const Real& rho, const tf_t& tf, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt, - Real& fr1, Real& fr1dt, Real& rr1, Real& rr1dt, - Real& fr2, Real& fr2dt, Real& rr2, Real& rr2dt, - Real& fr3, Real& fr3dt, Real& rr3, Real& rr3dt) + void postprocess_rate (const rhs_state_t& state, rate_t& rates, + rate_t& rates1, [[maybe_unused]] rate_t& rates2, [[maybe_unused]] rate_t& rates3) { using namespace Species; using namespace Rates; @@ -314,29 +308,29 @@ namespace RHS { if constexpr (rate == Si28_7He4_to_Ni56) { // first rate corresponds to Ca40_He4_to_Ti44 - if (tf.t9 > 2.5_rt && y(C12) + y(O16) <= 4.0e-3_rt) { + if (state.tf.t9 > 2.5_rt && state.y(C12) + state.y(O16) <= 4.0e-3_rt) { - Real t992 = tf.t972 * tf.t9; + Real t992 = state.tf.t972 * state.tf.t9; Real t9i92 = 1.0_rt / t992; - Real yeff_ca40 = t9i92 * std::exp(239.42_rt * tf.t9i - 74.741_rt); - Real yeff_ca40dt = -yeff_ca40 * (239.42_rt * tf.t9i2 + 4.5_rt * tf.t9i); + Real yeff_ca40 = t9i92 * std::exp(239.42_rt * state.tf.t9i - 74.741_rt); + Real yeff_ca40dt = -yeff_ca40 * (239.42_rt * state.tf.t9i2 + 4.5_rt * state.tf.t9i); - Real yeff_ti44 = t992 * std::exp(-274.12_rt * tf.t9i + 74.914_rt); - Real yeff_ti44dt = yeff_ti44*(274.12_rt * tf.t9i2 + 4.5_rt * tf.t9i); + Real yeff_ti44 = t992 * std::exp(-274.12_rt * state.tf.t9i + 74.914_rt); + Real yeff_ti44dt = yeff_ti44*(274.12_rt * state.tf.t9i2 + 4.5_rt * state.tf.t9i); - Real denom = std::pow(rho * y(He4), 3.0_rt); + Real denom = std::pow(state.rho * state.y(He4), 3.0_rt); - fr = yeff_ca40 * denom * fr1; - frdt = (yeff_ca40dt * fr1 + yeff_ca40 * fr1dt) * denom * 1.0e-9_rt; + rates.fr = yeff_ca40 * denom * rates1.fr; + rates.frdt = (yeff_ca40dt * rates1.fr + yeff_ca40 * rates1.frdt) * denom * 1.0e-9_rt; Real zz = 1.0_rt / denom; - rr = amrex::min(1.0e10_rt, yeff_ti44 * rr1 * zz); + rates.rr = amrex::min(1.0e10_rt, yeff_ti44 * rates1.rr * zz); - if (rr == 1.0e10_rt) { - rrdt = 0.0_rt; + if (rates.rr == 1.0e10_rt) { + rates.rrdt = 0.0_rt; } else { - rrdt = (yeff_ti44dt * rr1 + yeff_ti44 * rr1dt) * zz * 1.0e-9_rt; + rates.rrdt = (yeff_ti44dt * rates1.rr + yeff_ti44 * rates1.rrdt) * zz * 1.0e-9_rt; } } } diff --git a/networks/rhs.H b/networks/rhs.H index a6e3d7d01a..3a5df51c6d 100644 --- a/networks/rhs.H +++ b/networks/rhs.H @@ -26,52 +26,8 @@ namespace RHS { // Rate tabulation data. -constexpr Real tab_tlo = 6.0e0_rt; -constexpr Real tab_thi = 10.0e0_rt; -constexpr int tab_per_decade = 2000; -constexpr int nrattab = static_cast(tab_thi - tab_tlo) * tab_per_decade + 1; -constexpr int tab_imax = static_cast(tab_thi - tab_tlo) * tab_per_decade + 1; -constexpr Real tab_tstp = (tab_thi - tab_tlo) / static_cast(tab_imax - 1); - extern AMREX_GPU_MANAGED Array3D rattab; extern AMREX_GPU_MANAGED Array3D drattabdt; -extern AMREX_GPU_MANAGED Array1D ttab; - -// A struct that contains the terms needed to evaluate a tabulated rate. -struct rate_tab_t -{ - Real alfa, beta, gama, delt; - int iat; - - AMREX_GPU_HOST_DEVICE AMREX_INLINE - void initialize (Real temp) { - // hash locate - constexpr int mp = 4; - iat = static_cast((std::log10(temp) - tab_tlo) / tab_tstp) + 1; - iat = amrex::max(1, amrex::min(iat - 1, tab_imax - mp + 1)); - - // setup the lagrange interpolation coefficients for a cubic - Real x = temp; - Real x1 = ttab(iat); - Real x2 = ttab(iat+1); - Real x3 = ttab(iat+2); - Real x4 = ttab(iat+3); - Real a = x - x1; - Real b = x - x2; - Real c = x - x3; - Real d = x - x4; - Real e = x1 - x2; - Real f = x1 - x3; - Real g = x1 - x4; - Real h = x2 - x3; - Real p = x2 - x4; - Real q = x3 - x4; - alfa = b * c * d / (e * f * g); - beta = -a * c * d / (e * h * p); - gama = a * b * d / (f * h * q); - delt = -a * b * c / (g * p * q); - } -}; // Calculate an integer factorial. AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -393,13 +349,13 @@ void dgefa (RArray2D& a1) constexpr_for<0, INT_NEQS-1>([&] (auto n1) { - constexpr int k = n1; + [[maybe_unused]] constexpr int k = n1; // compute multipliers Real t = -1.0_rt / a(k,k); constexpr_for([&] (auto n2) { - constexpr int j = n2; + [[maybe_unused]] constexpr int j = n2; if (is_jacobian_term_used()) { a(j,k) *= t; @@ -409,12 +365,12 @@ void dgefa (RArray2D& a1) // row elimination with column indexing constexpr_for([&] (auto n2) { - constexpr int j = n2; + [[maybe_unused]] constexpr int j = n2; t = a(k,j); constexpr_for([&] (auto n3) { - constexpr int i = n3; + [[maybe_unused]] constexpr int i = n3; if constexpr (is_jacobian_term_used()) { a(i,j) += t * a(i,k); @@ -481,7 +437,7 @@ constexpr int density_exponent_reverse () // Scale a rate using the density terms. template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void apply_density_scaling (const Real& rho, Real& fr, Real& frdt, Real& rr, Real& rrdt) +void apply_density_scaling (const rhs_state_t& state, rate_t& rates) { constexpr int forward_exponent = density_exponent_forward(); constexpr int reverse_exponent = density_exponent_reverse(); @@ -494,36 +450,36 @@ void apply_density_scaling (const Real& rho, Real& fr, Real& frdt, Real& rr, Rea Real density_term_reverse = 1.0; if constexpr (forward_exponent == 1) { - density_term_forward = rho; + density_term_forward = state.rho; } else if constexpr (forward_exponent == 2) { - density_term_forward = rho * rho; + density_term_forward = state.rho * state.rho; } else if constexpr (forward_exponent == 3) { - density_term_forward = rho * rho * rho; + density_term_forward = state.rho * state.rho * state.rho; } if constexpr (reverse_exponent == 1) { - density_term_reverse = rho; + density_term_reverse = state.rho; } else if constexpr (reverse_exponent == 2) { - density_term_reverse = rho * rho; + density_term_reverse = state.rho * state.rho; } else if constexpr (reverse_exponent == 3) { - density_term_reverse = rho * rho * rho; + density_term_reverse = state.rho * state.rho * state.rho; } - fr *= density_term_forward; - frdt *= density_term_forward; - rr *= density_term_reverse; - rrdt *= density_term_reverse; + rates.fr *= density_term_forward; + rates.frdt *= density_term_forward; + rates.rr *= density_term_reverse; + rates.rrdt *= density_term_reverse; } #ifdef SCREENING // Apply the screening term to a given rate. template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& rr, Real& rrdt) +void apply_screening (const rhs_state_t& state, rate_t& rates) { // The screening behavior depends on the type of reaction. We provide screening // here for the reaction classes we know about, and any other reactions are unscreened. @@ -550,16 +506,16 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& static_assert(scn_fac.z1 == Z1); Real sc, scdt; - actual_screen(pstate, scn_fac, sc, scdt); + actual_screen(state.pstate, scn_fac, sc, scdt); if constexpr (data.screen_forward_reaction == 1) { - frdt = frdt * sc + fr * scdt; - fr = fr * sc; + rates.frdt = rates.frdt * sc + rates.fr * scdt; + rates.fr = rates.fr * sc; } if constexpr (data.screen_reverse_reaction == 1) { - rrdt = rrdt * sc + rr * scdt; - rr = rr * sc; + rates.rrdt = rates.rrdt * sc + rates.rr * scdt; + rates.rr = rates.rr * sc; } } @@ -574,16 +530,16 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& static_assert(scn_fac.z1 == Z1); Real sc, scdt; - actual_screen(pstate, scn_fac, sc, scdt); + actual_screen(state.pstate, scn_fac, sc, scdt); if constexpr (data.screen_forward_reaction == 1) { - frdt = frdt * sc + fr * scdt; - fr = fr * sc; + rates.frdt = rates.frdt * sc + rates.fr * scdt; + rates.fr = rates.fr * sc; } if constexpr (data.screen_reverse_reaction == 1) { - rrdt = rrdt * sc + rr * scdt; - rr = rr * sc; + rates.rrdt = rates.rrdt * sc + rates.rr * scdt; + rates.rr = rates.rr * sc; } } @@ -599,7 +555,7 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& static_assert(scn_fac1.z1 == Z1); Real sc1, sc1dt; - actual_screen(pstate, scn_fac1, sc1, sc1dt); + actual_screen(state.pstate, scn_fac1, sc1, sc1dt); constexpr Real Z2 = 2.0_rt * Z1; constexpr Real A2 = 2.0_rt * A1; @@ -609,7 +565,7 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& static_assert(scn_fac2.z1 == Z1); Real sc2, sc2dt; - actual_screen(pstate, scn_fac2, sc2, sc2dt); + actual_screen(state.pstate, scn_fac2, sc2, sc2dt); // Compute combined screening factor @@ -617,13 +573,13 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& Real scdt = sc1dt * sc2 + sc1 * sc2dt; if constexpr (data.screen_forward_reaction == 1) { - frdt = frdt * sc + fr * scdt; - fr = fr * sc; + rates.frdt = rates.frdt * sc + rates.fr * scdt; + rates.fr = rates.fr * sc; } if constexpr (data.screen_reverse_reaction == 1) { - rrdt = rrdt * sc + rr * scdt; - rr = rr * sc; + rates.rrdt = rates.rrdt * sc + rates.rr * scdt; + rates.rr = rates.rr * sc; } } } @@ -632,24 +588,24 @@ void apply_screening (const plasma_state_t& pstate, Real& fr, Real& frdt, Real& // Apply the branching ratios to a given rate. template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void apply_branching (Real& fr, Real& frdt, Real& rr, Real& rrdt) +void apply_branching (rate_t& rates) { constexpr rhs_t data = rhs_data(rate); if constexpr (data.forward_branching_ratio != 1.0_rt) { - fr *= data.forward_branching_ratio; - frdt *= data.forward_branching_ratio; + rates.fr *= data.forward_branching_ratio; + rates.frdt *= data.forward_branching_ratio; } if constexpr (data.reverse_branching_ratio != 1.0_rt) { - rr *= data.reverse_branching_ratio; - rrdt *= data.reverse_branching_ratio; + rates.rr *= data.reverse_branching_ratio; + rates.rrdt *= data.reverse_branching_ratio; } } // Do the initial tabulation of rates. We loop over various // temperatures and evaluate all the rates at each temperature. -AMREX_INLINE +AMREX_GPU_HOST_DEVICE AMREX_INLINE void tabulate_rates () { using namespace Rates; @@ -660,35 +616,39 @@ void tabulate_rates () ttab(i) = temp; + rhs_state_t state; + // Get the temperature factors - tf_t tf = get_tfactors(temp); + state.tf = get_tfactors(temp); // Arbitrary density, y, and y_e values (should be unused) - const Real rho = 0.0_rt; - const Real y_e = 0.0_rt; - const Real eta = 0.0_rt; - Array1D y = {0.0_rt}; + state.rho = 0.0_rt; + state.y_e = 0.0_rt; + state.eta = 0.0_rt; + state. y = {0.0_rt}; constexpr_for<1, NumRates+1>([&] (auto n) { - constexpr int rate = n; + [[maybe_unused]] constexpr int rate = n; - Real fr = 0.0_rt; - Real frdt = 0.0_rt; - - Real rr = 0.0_rt; - Real rrdt = 0.0_rt; + rate_t rates; constexpr rhs_t data = RHS::rhs_data(rate); if constexpr (data.rate_can_be_tabulated) { - evaluate_analytical_rate(rho, tf, y_e, eta, y, fr, frdt, rr, rrdt); + evaluate_analytical_rate(state, rates); + } + else { + rates.fr = 0.0_rt; + rates.rr = 0.0_rt; + rates.frdt = 0.0_rt; + rates.rrdt = 0.0_rt; } - rattab(rate, 1, i) = fr; - rattab(rate, 2, i) = rr; - drattabdt(rate, 1, i) = frdt; - drattabdt(rate, 2, i) = rrdt; + rattab(rate, 1, i) = rates.fr; + rattab(rate, 2, i) = rates.rr; + drattabdt(rate, 1, i) = rates.frdt; + drattabdt(rate, 2, i) = rates.rrdt; }); } } @@ -696,29 +656,27 @@ void tabulate_rates () // Evaluate a rate using the rate tables. template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void evaluate_tabulated_rate (const Real& rho, const rate_tab_t& tab, const Real& y_e, - const Real& eta, const Array1D& y, - Real& fr, Real& frdt, Real& rr, Real& rrdt) +void evaluate_tabulated_rate (const rhs_state_t& state, rate_t& rates) { - fr = (tab.alfa * rattab(rate, 1, tab.iat ) + - tab.beta * rattab(rate, 1, tab.iat+1) + - tab.gama * rattab(rate, 1, tab.iat+2) + - tab.delt * rattab(rate, 1, tab.iat+3)); - - rr = (tab.alfa * rattab(rate, 2, tab.iat ) + - tab.beta * rattab(rate, 2, tab.iat+1) + - tab.gama * rattab(rate, 2, tab.iat+2) + - tab.delt * rattab(rate, 2, tab.iat+3)); - - frdt = (tab.alfa * drattabdt(rate, 1, tab.iat ) + - tab.beta * drattabdt(rate, 1, tab.iat+1) + - tab.gama * drattabdt(rate, 1, tab.iat+2) + - tab.delt * drattabdt(rate, 1, tab.iat+3)); - - rrdt = (tab.alfa * drattabdt(rate, 2, tab.iat ) + - tab.beta * drattabdt(rate, 2, tab.iat+1) + - tab.gama * drattabdt(rate, 2, tab.iat+2) + - tab.delt * drattabdt(rate, 2, tab.iat+3)); + rates.fr = (state.tab.alfa * rattab(rate, 1, state.tab.iat ) + + state.tab.beta * rattab(rate, 1, state.tab.iat+1) + + state.tab.gama * rattab(rate, 1, state.tab.iat+2) + + state.tab.delt * rattab(rate, 1, state.tab.iat+3)); + + rates.rr = (state.tab.alfa * rattab(rate, 2, state.tab.iat ) + + state.tab.beta * rattab(rate, 2, state.tab.iat+1) + + state.tab.gama * rattab(rate, 2, state.tab.iat+2) + + state.tab.delt * rattab(rate, 2, state.tab.iat+3)); + + rates.frdt = (state.tab.alfa * drattabdt(rate, 1, state.tab.iat ) + + state.tab.beta * drattabdt(rate, 1, state.tab.iat+1) + + state.tab.gama * drattabdt(rate, 1, state.tab.iat+2) + + state.tab.delt * drattabdt(rate, 1, state.tab.iat+3)); + + rates.rrdt = (state.tab.alfa * drattabdt(rate, 2, state.tab.iat ) + + state.tab.beta * drattabdt(rate, 2, state.tab.iat+1) + + state.tab.gama * drattabdt(rate, 2, state.tab.iat+2) + + state.tab.delt * drattabdt(rate, 2, state.tab.iat+3)); } // Calculate the RHS term for a given species and rate. @@ -1215,49 +1173,43 @@ constexpr Real jac_term (const burn_t& state, const Real& fr, const Real& rr) template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void construct_rate (const Real& rho, const tf_t& tf, const rate_tab_t& tab, const Real& y_e, - const Real& eta, const Array1D& y, const plasma_state_t& pstate, - Real& fr, Real& frdt, Real& rr, Real& rrdt) +void construct_rate (const rhs_state_t& state, rate_t& rates) { using namespace Species; using namespace Rates; - fr = 0.0; - frdt = 0.0; - rr = 0.0; - rrdt = 0.0; + rates.fr = 0.0; + rates.frdt = 0.0; + rates.rr = 0.0; + rates.rrdt = 0.0; constexpr rhs_t data = RHS::rhs_data(rate); if (use_tables && data.rate_can_be_tabulated) { - evaluate_tabulated_rate(rho, tab, y_e, eta, y, fr, frdt, rr, rrdt); + evaluate_tabulated_rate(state, rates); } else { - evaluate_analytical_rate(rho, tf, y_e, eta, y, fr, frdt, rr, rrdt); + evaluate_analytical_rate(state, rates); } // Set the density dependence - apply_density_scaling(rho, fr, frdt, rr, rrdt); + apply_density_scaling(state, rates); #ifdef SCREENING // Screen - apply_screening(pstate, fr, frdt, rr, rrdt); + apply_screening(state, rates); #endif // Branching ratios - apply_branching(fr, frdt, rr, rrdt); + apply_branching(rates); } template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void fill_additional_rates (Arr intermediate_fr, Arr intermediate_frdt, - Arr intermediate_rr, Arr intermediate_rrdt, - Real& fr1, Real& fr1dt, Real& rr1, Real& rr1dt, - Real& fr2, Real& fr2dt, Real& rr2, Real& rr2dt, - Real& fr3, Real& fr3dt, Real& rr3, Real& rr3dt) +void fill_additional_rates (const Arr& intermediate_rates, rate_t& rates1, rate_t& rates2, rate_t& rates3) { constexpr rhs_t data = RHS::rhs_data(rate); @@ -1270,10 +1222,7 @@ void fill_additional_rates (Arr intermediate_fr, Arr intermediate_frdt, static_assert(index >= 1 && index <= num_intermediate_reactions()); - fr1 = intermediate_fr(index); - fr1dt = intermediate_frdt(index); - rr1 = intermediate_rr(index); - rr1dt = intermediate_rrdt(index); + rates1 = intermediate_rates(index); } if constexpr (rate2 >= 0) { @@ -1281,10 +1230,7 @@ void fill_additional_rates (Arr intermediate_fr, Arr intermediate_frdt, static_assert(index >= 1 && index <= num_intermediate_reactions()); - fr2 = intermediate_fr(index); - fr2dt = intermediate_frdt(index); - rr2 = intermediate_rr(index); - rr2dt = intermediate_rrdt(index); + rates2 = intermediate_rates(index); } if constexpr (rate3 >= 0) { @@ -1292,10 +1238,7 @@ void fill_additional_rates (Arr intermediate_fr, Arr intermediate_frdt, static_assert(index >= 1 && index <= num_intermediate_reactions()); - fr3 = intermediate_fr(index); - fr3dt = intermediate_frdt(index); - rr3 = intermediate_rr(index); - rr3dt = intermediate_rrdt(index); + rates3 = intermediate_rates(index); } } @@ -1317,25 +1260,28 @@ void rhs_init () // for each term in ydot). template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void rhs (burn_t& state, Array1D& ydot) +void rhs (burn_t& burn_state, Array1D& ydot) { static_assert(nrhs == neqs || nrhs == 2 * neqs); + rhs_state_t rhs_state; + + rhs_state.rho = burn_state.rho; + rhs_state.eta = burn_state.eta; + rhs_state.y_e = burn_state.y_e; + // Convert X to Y. - Array1D y; for (int n = 1; n <= NumSpec; ++n) { - y(n) = state.xn[n-1] * aion_inv[n-1]; + rhs_state.y(n) = burn_state.xn[n-1] * aion_inv[n-1]; } // Set up the state data, which is the same for all screening factors. - plasma_state_t pstate; - fill_plasma_state(pstate, state.T, state.rho, y); + fill_plasma_state(rhs_state.pstate, burn_state.T, burn_state.rho, rhs_state.y); // Initialize the rate temperature term. - tf_t tf = get_tfactors(state.T); - rate_tab_t tab; + rhs_state.tf = get_tfactors(burn_state.T); if (use_tables) { - tab.initialize(state.T); + rhs_state.tab.initialize(burn_state.T); } // Initialize the RHS terms. @@ -1350,10 +1296,7 @@ void rhs (burn_t& state, Array1D& ydot) constexpr int intermediate_array_size = num_intermediate > 0 ? num_intermediate : 1; // Define forward and reverse (and d/dT) rate arrays. - Array1D intermediate_fr; - Array1D intermediate_frdt; - Array1D intermediate_rr; - Array1D intermediate_rrdt; + Array1D intermediate_rates; // Fill all intermediate rates first. constexpr_for<1, Rates::NumRates+1>([&] (auto n) @@ -1363,11 +1306,7 @@ void rhs (burn_t& state, Array1D& ydot) constexpr int index = locate_intermediate_rate_index(rate); if constexpr (index >= 1) { - construct_rate(state.rho, tf, tab, state.y_e, state.eta, y, pstate, - intermediate_fr(index), - intermediate_frdt(index), - intermediate_rr(index), - intermediate_rrdt(index)); + construct_rate(rhs_state, intermediate_rates(index)); } }); @@ -1377,51 +1316,38 @@ void rhs (burn_t& state, Array1D& ydot) { constexpr int rate = n1; - Real fr, frdt, rr, rrdt; + rate_t rates; // We only need to compute the rate at this point if it's not intermediate. If it // is intermediate, retrieve it from the cached array. constexpr int index = locate_intermediate_rate_index(rate); if constexpr (index < 0) { - construct_rate(state.rho, tf, tab, state.y_e, state.eta, y, pstate, fr, frdt, rr, rrdt); + construct_rate(rhs_state, rates); } else { - fr = intermediate_fr(index); - frdt = intermediate_frdt(index); - rr = intermediate_rr(index); - rrdt = intermediate_rrdt(index); + rates = intermediate_rates(index); } // Locate all intermediate rates needed to augment this reaction. // To keep the problem bounded we assume that there are no more than // three intermediate reactions needed. - Real fr1, fr1dt, rr1, rr1dt; - Real fr2, fr2dt, rr2, rr2dt; - Real fr3, fr3dt, rr3, rr3dt; + rate_t rates1, rates2, rates3; - fill_additional_rates(intermediate_fr, intermediate_frdt, - intermediate_rr, intermediate_rrdt, - fr1, fr1dt, rr1, rr1dt, - fr2, fr2dt, rr2, rr2dt, - fr3, fr3dt, rr3, rr3dt); + fill_additional_rates(intermediate_rates, rates1, rates2, rates3); // Perform rate postprocessing, using additional reactions as inputs. // If there is no postprocessing for this rate, this will be a no-op. - postprocess_rate(state.rho, tf, state.y_e, state.eta, y, - fr, frdt, rr, rrdt, - fr1, fr1dt, rr1, rr1dt, - fr2, fr2dt, rr2, rr2dt, - fr3, fr3dt, rr3, rr3dt); + postprocess_rate(rhs_state, rates, rates1, rates2, rates3); constexpr_for<1, NumSpec+1>([&] (auto n2) { constexpr int species = n2; if constexpr (is_rate_used()) { - auto [forward_term, reverse_term] = rhs_term(state, fr, rr); + auto [forward_term, reverse_term] = rhs_term(burn_state, rates.fr, rates.rr); if constexpr (nrhs == 2 * neqs) { if (forward_term >= 0.0_rt) { @@ -1443,7 +1369,8 @@ void rhs (burn_t& state, Array1D& ydot) // Evaluate the neutrino cooling. #ifdef NEUTRINOS Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, + sneut, dsneutdt, dsneutdd, snuda, snudz); #else Real sneut = 0.0; #endif @@ -1473,23 +1400,26 @@ void rhs (burn_t& state, Array1D& ydot) // Analytical Jacobian AMREX_GPU_HOST_DEVICE AMREX_INLINE -void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) +void jac (burn_t& burn_state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) { + rhs_state_t rhs_state; + + rhs_state.rho = burn_state.rho; + rhs_state.eta = burn_state.eta; + rhs_state.y_e = burn_state.y_e; + // Convert X to Y. - Array1D y; for (int n = 1; n <= NumSpec; ++n) { - y(n) = state.xn[n-1] * aion_inv[n-1]; + rhs_state.y(n) = burn_state.xn[n-1] * aion_inv[n-1]; } // Set up the state data, which is the same for all screening factors. - plasma_state_t pstate; - fill_plasma_state(pstate, state.T, state.rho, y); + fill_plasma_state(rhs_state.pstate, burn_state.T, burn_state.rho, rhs_state.y); // Initialize the rate temperature term. - tf_t tf = get_tfactors(state.T); - rate_tab_t tab; + rhs_state.tf = get_tfactors(burn_state.T); if (use_tables) { - tab.initialize(state.T); + rhs_state.tab.initialize(burn_state.T); } // Initialize the Jacobian terms. @@ -1506,14 +1436,9 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) constexpr int intermediate_array_size = num_intermediate > 0 ? num_intermediate : 1; // Define forward and reverse (and d/dT) rate arrays. - Array1D intermediate_fr; - Array1D intermediate_frdt; - Array1D intermediate_rr; - Array1D intermediate_rrdt; + Array1D intermediate_rates; - Real fr1, fr1dt, rr1, rr1dt; - Real fr2, fr2dt, rr2, rr2dt; - Real fr3, fr3dt, rr3, rr3dt; + rate_t rates1, rates2, rates3; // Fill all intermediate rates first. constexpr_for<1, Rates::NumRates+1>([&] (auto n) @@ -1523,11 +1448,7 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) constexpr int index = locate_intermediate_rate_index(rate); if constexpr (index >= 1) { - construct_rate(state.rho, tf, tab, state.y_e, state.eta, y, pstate, - intermediate_fr(index), - intermediate_frdt(index), - intermediate_rr(index), - intermediate_rrdt(index)); + construct_rate(rhs_state, intermediate_rates(index)); } }); @@ -1536,51 +1457,41 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) { constexpr int rate = n1; - Real fr, frdt, rr, rrdt; + rate_t rates; // We only need to compute the rate at this point if it's not intermediate. If it // is intermediate, retrieve it from the cached array. constexpr int index = locate_intermediate_rate_index(rate); if constexpr (index < 0) { - construct_rate(state.rho, tf, tab, state.y_e, state.eta, y, pstate, fr, frdt, rr, rrdt); + construct_rate(rhs_state, rates); } else { - fr = intermediate_fr(index); - frdt = intermediate_frdt(index); - rr = intermediate_rr(index); - rrdt = intermediate_rrdt(index); + rates = intermediate_rates(index); } // Locate all intermediate rates needed to augment this reaction. // To keep the problem bounded we assume that there are no more than // three intermediate reactions needed. - fill_additional_rates(intermediate_fr, intermediate_frdt, - intermediate_rr, intermediate_rrdt, - fr1, fr1dt, rr1, rr1dt, - fr2, fr2dt, rr2, rr2dt, - fr3, fr3dt, rr3, rr3dt); + fill_additional_rates(intermediate_rates, rates1, rates2, rates3); // Perform rate postprocessing, using additional reactions as inputs. // If there is no postprocessing for this rate, this will be a no-op. - postprocess_rate(state.rho, tf, state.y_e, state.eta, y, - fr, frdt, rr, rrdt, - fr1, fr1dt, rr1, rr1dt, - fr2, fr2dt, rr2, rr2dt, - fr3, fr3dt, rr3, rr3dt); + postprocess_rate(rhs_state, rates, rates1, rates2, rates3); // Species Jacobian elements with respect to other species. constexpr_for<1, NumSpec+1>([&] (auto n2) { - constexpr int spec1 = n2; + [[maybe_unused]] constexpr int spec1 = n2; constexpr_for<1, NumSpec+1>([&] (auto n3) { + [[maybe_unused]] constexpr int spec2 = n3; + if constexpr (is_rate_used()) { - constexpr int spec2 = n3; - jac(spec1, spec2) += jac_term(state, fr, rr); + jac(spec1, spec2) += jac_term(burn_state, rates.fr, rates.rr); } }); }); @@ -1589,10 +1500,10 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) // We'll convert them from d/dT to d/de later. constexpr_for<1, NumSpec+1>([&] (auto n2) { - constexpr int species = n2; + [[maybe_unused]] constexpr int species = n2; if constexpr (is_rate_used()) { - auto [forward_term, reverse_term] = rhs_term(state, frdt, rrdt); + auto [forward_term, reverse_term] = rhs_term(burn_state, rates.frdt, rates.rrdt); jac(species, net_ienuc) += forward_term + reverse_term; } }); @@ -1601,20 +1512,20 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) // Evaluate the neutrino cooling. #ifdef NEUTRINOS Real sneut, dsneutdt, dsneutdd, snuda, snudz; - sneut5(state.T, state.rho, state.abar, state.zbar, + sneut5(burn_state.T, burn_state.rho, burn_state.abar, burn_state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz); #else Real sneut = 0.0, dsneutdt = 0.0, dsneutdd = 0.0, snuda = 0.0, snudz = 0.0; #endif - jac(net_ienuc, net_ienuc) = -temperature_to_energy_jacobian(state, dsneutdt); + jac(net_ienuc, net_ienuc) = -temperature_to_energy_jacobian(burn_state, dsneutdt); constexpr_for<1, NumSpec+1>([&] (auto j) { constexpr int species = j; // Energy generation rate Jacobian elements with respect to species. - Real b1 = (-state.abar * state.abar * snuda + (NetworkProperties::zion(species) - state.zbar) * state.abar * snudz); + Real b1 = (-burn_state.abar * burn_state.abar * snuda + (NetworkProperties::zion(species) - burn_state.zbar) * burn_state.abar * snudz); jac(net_ienuc, species) = -b1; constexpr_for<1, NumSpec+1>([&] (auto i) @@ -1625,7 +1536,7 @@ void jac (burn_t& state, ArrayUtil::MathArray2D<1, neqs, 1, neqs>& jac) }); // Convert previously computed terms from d/dT to d/de. - jac(species, net_ienuc) = temperature_to_energy_jacobian(state, jac(species, net_ienuc)); + jac(species, net_ienuc) = temperature_to_energy_jacobian(burn_state, jac(species, net_ienuc)); // Compute df(e) / de term. jac(net_ienuc, net_ienuc) += ener_gener_rate(jac(species, net_ienuc)); diff --git a/rates/Make.package b/rates/Make.package index 593b2f5ac3..2c78534c31 100644 --- a/rates/Make.package +++ b/rates/Make.package @@ -1,4 +1,3 @@ CEXE_headers += aprox_rates_data.H CEXE_sources += aprox_rates_data.cpp CEXE_headers += aprox_rates.H -CEXE_headers += tfactors.H