Skip to content

Commit

Permalink
update He-C-Fe-group with new pynucastro (#1364)
Browse files Browse the repository at this point in the history 
we now need to remove duplicates before the RateCollection
  • Loading branch information
@zingale
zingale authored Oct 22, 2023
1 parent b6838d2 commit 8ec4f08
Show file tree
Hide file tree
Showing 6 changed files with 231 additions and 227 deletions.
Binary file modified networks/He-C-Fe-group/He-C-Fe-group.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
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
Loading

0 comments on commit 8ec4f08

Please sign in to comment.