Skip to content

Commit

Permalink
Merge branch 'development' into update_nuc_capitalize
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhichen3
zhichen3 committed Oct 23, 2023
2 parents 5a4e367 + 8ec4f08 commit 546435b
Show file tree
Hide file tree
Showing 5 changed files with 42 additions and 39 deletions.
16 changes: 9 additions & 7 deletions networks/He-C-Fe-group/He-C-Fe-group.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,23 +95,25 @@
d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True)
all_lib.add_rate(d)

# combine all three libraries into a single network

net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib],
symmetric_screening=False)

# we will have duplicate rates -- we want to remove any ReacLib rates
# that we have tabular rates for

dupes = net.find_duplicate_links()
dupes = all_lib.find_duplicate_links()

rates_to_remove = []
for d in dupes:
for r in d:
if isinstance(r, ReacLibRate):
rates_to_remove.append(r)

net.remove_rates(rates_to_remove)
for r in rates_to_remove:
all_lib.remove_rate(r)

# combine all three libraries into a single network

net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib],
symmetric_screening=False)


# now we approximate some (alpha, p)(p, gamma) links

Expand Down
43 changes: 18 additions & 25 deletions networks/He-C-Fe-group/actual_rhs.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1269,101 +1269,103 @@ 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_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co55_to_Fe55) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Co55_to_Fe55) = drate_dt;
}
rate_eval.add_energy_rate(k_Co55_to_Fe55) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co56_Fe56_meta, j_Co56_Fe56_rhoy, j_Co56_Fe56_temp, j_Co56_Fe56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co56_to_Fe56) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Co56_to_Fe56) = drate_dt;
}
rate_eval.add_energy_rate(k_Co56_to_Fe56) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co56_Ni56_meta, j_Co56_Ni56_rhoy, j_Co56_Ni56_temp, j_Co56_Ni56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co56_to_Ni56) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Co56_to_Ni56) = drate_dt;
}
rate_eval.add_energy_rate(k_Co56_to_Ni56) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Co57_Ni57_meta, j_Co57_Ni57_rhoy, j_Co57_Ni57_temp, j_Co57_Ni57_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Co57_to_Ni57) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Co57_to_Ni57) = drate_dt;
}
rate_eval.add_energy_rate(k_Co57_to_Ni57) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Co57) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe55_Co55_meta, j_Fe55_Co55_rhoy, j_Fe55_Co55_temp, j_Fe55_Co55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe55_to_Co55) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Fe55_to_Co55) = drate_dt;
}
rate_eval.add_energy_rate(k_Fe55_to_Co55) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe55_Mn55_meta, j_Fe55_Mn55_rhoy, j_Fe55_Mn55_temp, j_Fe55_Mn55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe55_to_Mn55) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Fe55_to_Mn55) = drate_dt;
}
rate_eval.add_energy_rate(k_Fe55_to_Mn55) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Fe56_Co56_meta, j_Fe56_Co56_rhoy, j_Fe56_Co56_temp, j_Fe56_Co56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Fe56_to_Co56) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Fe56_to_Co56) = drate_dt;
}
rate_eval.add_energy_rate(k_Fe56_to_Co56) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Mn55_Fe55_meta, j_Mn55_Fe55_rhoy, j_Mn55_Fe55_temp, j_Mn55_Fe55_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Mn55_to_Fe55) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Mn55_to_Fe55) = drate_dt;
}
rate_eval.add_energy_rate(k_Mn55_to_Fe55) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Mn55) * (edot_nu + edot_gamma);

tabular_evaluate(j_Ni56_Co56_meta, j_Ni56_Co56_rhoy, j_Ni56_Co56_temp, j_Ni56_Co56_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Ni56_to_Co56) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Ni56_to_Co56) = drate_dt;
}
rate_eval.add_energy_rate(k_Ni56_to_Co56) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni56) * (edot_nu + edot_gamma);

tabular_evaluate(j_Ni57_Co57_meta, j_Ni57_Co57_rhoy, j_Ni57_Co57_temp, j_Ni57_Co57_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_Ni57_to_Co57) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_Ni57_to_Co57) = drate_dt;
}
rate_eval.add_energy_rate(k_Ni57_to_Co57) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni57) * (edot_nu + edot_gamma);

tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_n_to_p) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_n_to_p) = drate_dt;
}
rate_eval.add_energy_rate(k_n_to_p) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma);

tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_p_to_n) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_p_to_n) = drate_dt;
}
rate_eval.add_energy_rate(k_p_to_n) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);


}
Expand Down Expand Up @@ -1724,6 +1726,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 @@ -1733,19 +1736,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(Co55) * rate_eval.add_energy_rate(k_Co55_to_Fe55);
enuc += C::Legacy::n_A * Y(Co56) * rate_eval.add_energy_rate(k_Co56_to_Fe56);
enuc += C::Legacy::n_A * Y(Co56) * rate_eval.add_energy_rate(k_Co56_to_Ni56);
enuc += C::Legacy::n_A * Y(Co57) * rate_eval.add_energy_rate(k_Co57_to_Ni57);
enuc += C::Legacy::n_A * Y(Fe55) * rate_eval.add_energy_rate(k_Fe55_to_Co55);
enuc += C::Legacy::n_A * Y(Fe55) * rate_eval.add_energy_rate(k_Fe55_to_Mn55);
enuc += C::Legacy::n_A * Y(Fe56) * rate_eval.add_energy_rate(k_Fe56_to_Co56);
enuc += C::Legacy::n_A * Y(Mn55) * rate_eval.add_energy_rate(k_Mn55_to_Fe55);
enuc += C::Legacy::n_A * Y(Ni56) * rate_eval.add_energy_rate(k_Ni56_to_Co56);
enuc += C::Legacy::n_A * Y(Ni57) * rate_eval.add_energy_rate(k_Ni57_to_Co57);
enuc += C::Legacy::n_A * Y(N) * rate_eval.add_energy_rate(k_n_to_p);
enuc += C::Legacy::n_A * Y(H1) * rate_eval.add_energy_rate(k_p_to_n);
// include any weak rate neutrino losses
enuc += rate_eval.enuc_weak;

// Get the thermal neutrino losses

Expand Down Expand Up @@ -2770,6 +2762,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/He-C-Fe-group/reaclib_rates.H
View file
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/He-C-Fe-group/table_rates.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,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 @@ -139,6 +144,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 @@ -305,8 +314,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
5 changes: 2 additions & 3 deletions unit_test/react_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,8 +138,8 @@ void get_xn(const int k, const init_t cd, Real *xn_zone, int uniform_composition

for (int n = 0; n < NumSpec; n++) {
if (n == cd.is1 ||
(n == cd.is2 && cd.nprim >= 2) ||
(n == cd.is3 && cd.nprim >= 3)) {
(cd.nprim >= 2 && n == cd.is2) ||
(cd.nprim >= 3 && n == cd.is3)) {
continue;
}

Expand Down Expand Up @@ -286,4 +286,3 @@ Real get_xn(const int index) {
}

#endif

0 comments on commit 546435b

Please sign in to comment.