Skip to content

Commit

Permalink
update ECSN to latest pynucastro
Browse files Browse the repository at this point in the history
  • Loading branch information
zingale committed Oct 23, 2023
1 parent 8ec4f08 commit 89a1b8d
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 14 deletions.
19 changes: 10 additions & 9 deletions networks/ECSN/actual_rhs.H
Original file line number Diff line number Diff line change
Expand Up @@ -268,37 +268,39 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {

[[maybe_unused]] Real rate, drate_dt, edot_nu, edot_gamma;

rate_eval.enuc_weak = 0.0;

tabular_evaluate(j_f20_o20_meta, j_f20_o20_rhoy, j_f20_o20_temp, j_f20_o20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_f20_to_o20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_f20_to_o20) = drate_dt;
}
rate_eval.add_energy_rate(k_f20_to_o20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);

tabular_evaluate(j_ne20_f20_meta, j_ne20_f20_rhoy, j_ne20_f20_temp, j_ne20_f20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_ne20_to_f20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_ne20_to_f20) = drate_dt;
}
rate_eval.add_energy_rate(k_ne20_to_f20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne20) * (edot_nu + edot_gamma);

tabular_evaluate(j_o20_f20_meta, j_o20_f20_rhoy, j_o20_f20_temp, j_o20_f20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_o20_to_f20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_o20_to_f20) = drate_dt;
}
rate_eval.add_energy_rate(k_o20_to_f20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(O20) * (edot_nu + edot_gamma);

tabular_evaluate(j_f20_ne20_meta, j_f20_ne20_rhoy, j_f20_ne20_temp, j_f20_ne20_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_f20_to_ne20) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_f20_to_ne20) = drate_dt;
}
rate_eval.add_energy_rate(k_f20_to_ne20) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);


}
Expand Down Expand Up @@ -394,6 +396,7 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
rate_t rate_eval;

constexpr int do_T_derivatives = 0;

evaluate_rates<do_T_derivatives, rate_t>(state, rate_eval);

rhs_nuc(state, ydot, Y, rate_eval.screened_rates);
Expand All @@ -403,11 +406,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
Real enuc;
ener_gener_rate(ydot, enuc);

// include reaction neutrino losses (non-thermal) and gamma heating
enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_o20);
enuc += C::Legacy::n_A * Y(Ne20) * rate_eval.add_energy_rate(k_ne20_to_f20);
enuc += C::Legacy::n_A * Y(O20) * rate_eval.add_energy_rate(k_o20_to_f20);
enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_ne20);
// include any weak rate neutrino losses
enuc += rate_eval.enuc_weak;

// Get the thermal neutrino losses

Expand Down Expand Up @@ -613,6 +613,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
rate_derivs_t rate_eval;

constexpr int do_T_derivatives = 1;

evaluate_rates<do_T_derivatives, rate_derivs_t>(state, rate_eval);

// Species Jacobian elements with respect to other species
Expand Down
4 changes: 2 additions & 2 deletions networks/ECSN/reaclib_rates.H
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,13 @@ using namespace Species;

struct rate_t {
Array1D<Real, 1, NumRates> screened_rates;
Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
Real enuc_weak;
};

struct rate_derivs_t {
Array1D<Real, 1, NumRates> screened_rates;
Array1D<Real, 1, NumRates> dscreened_rates_dT;
Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
Real enuc_weak;
};


Expand Down
13 changes: 11 additions & 2 deletions networks/ECSN/table_rates.H
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,11 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
std::ifstream table;
table.open(file);

if (!table.is_open()) {
// the table was not present or we could not open it; abort
amrex::Error("table could not be opened");
}

std::string line;

// read and skip over the header
Expand All @@ -99,6 +104,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table
for (int j = 1; j <= tf.nrhoy; ++j) {
for (int i = 1; i <= tf.ntemp; ++i) {
std::getline(table, line);
if (line.empty()) {
amrex::Error("Error reading table data");
}

std::istringstream sdata(line);

sdata >> log_rhoy_table(j) >> log_temp_table(i);
Expand Down Expand Up @@ -265,8 +274,8 @@ evaluate_dr_dtemp(const table_t& table_meta, const R& log_rhoy_table, const T& l
// Finally, we perform a 1d-linear interpolation between dlogr_dlogt_ip1
// and dlogr_dlogt_i to compute dlogr_dlogt

Real log_rhoy_lo = log_temp_table(irhoy_lo);
Real log_rhoy_hi = log_temp_table(irhoy_hi);
Real log_rhoy_lo = log_rhoy_table(irhoy_lo);
Real log_rhoy_hi = log_rhoy_table(irhoy_hi);

Real log_temp_lo = log_temp_table(jtemp_lo);
Real log_temp_hi = log_temp_table(jtemp_hi);
Expand Down
2 changes: 1 addition & 1 deletion unit_test/test_rhs/ci-benchmarks/ecsn.out
Original file line number Diff line number Diff line change
Expand Up @@ -161,7 +161,7 @@
J_hydrogen-1_E -0.12580793235 1.2089803642e+12
J_helium-4_E -0.0078065301852 1.4905283338e+12
J_oxygen-16_E -5.3990173959e+12 1.2739583515e-06
J_oxygen-20_E -2.079874023e-18 3.5544574816e-18
J_oxygen-20_E -2.8880340474e-19 3.5544574816e-18
J_fluorine-20_E -2.0716943987e-18 7.6531304664e-18
J_neon-20_E -8.8836179369e-05 1.9862184663e-05
J_magnesium-24_E -0.0040169125096 0.018563573855
Expand Down

0 comments on commit 89a1b8d

Please sign in to comment.