Skip to content

Commit

Permalink
Construct some containers for data in rhs.H (#1054)
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@maxpkatz
maxpkatz authored Nov 27, 2022
1 parent bee474f commit ca6d256
Show file tree
Hide file tree
Showing 9 changed files with 634 additions and 674 deletions.
1 change: 1 addition & 0 deletions interfaces/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
71 changes: 71 additions & 0 deletions interfaces/rhs_type.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,11 @@
#define rhs_type_H

#include <AMReX_REAL.H>
#include <network_properties.H>
#ifdef SCREENING
#include <screen.H>
#endif
#include <tfactors.H>

namespace RHS
{
Expand Down Expand Up @@ -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<int>(tab_thi - tab_tlo) * tab_per_decade + 1;
constexpr int tab_imax = static_cast<int>(tab_thi - tab_tlo) * tab_per_decade + 1;
constexpr amrex::Real tab_tstp = (tab_thi - tab_tlo) / static_cast<Real>(tab_imax - 1);

extern AMREX_GPU_MANAGED Array1D<amrex::Real, 1, nrattab> 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<int>((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<amrex::Real, 1, NumSpec> y;
};

struct rate_t
{
amrex::Real fr;
amrex::Real rr;
amrex::Real frdt;
amrex::Real rrdt;
};

} // namespace RHS

#endif
2 changes: 2 additions & 0 deletions rates/tfactors.H → interfaces/tfactors.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
#include <AMReX.H>
#include <cmath>

using namespace amrex;

struct tf_t {
amrex::Real temp;
amrex::Real t9;
Expand Down
120 changes: 57 additions & 63 deletions networks/aprox13/actual_network.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -712,125 +712,119 @@ namespace RHS {

template<int rate>
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<Real, 1, NumSpec>& 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<int rate>
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<Real, 1, NumSpec>& 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;
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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;
}
}

Expand Down
Loading

0 comments on commit ca6d256

Please sign in to comment.