Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update Urca network with the p -> n e-capture rate from Langanke #1359

Merged
merged 8 commits into from
Oct 24, 2023
Merged

update Urca network with the p -> n e-capture rate from Langanke #1359

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
b6ae2c3
update the Urca network with the p -> n electron capture rate from La…
zingale Oct 18, 2023
9915744
add missing files
zingale Oct 18, 2023
7d185b4
Merge branch 'development' into update_urca
zingale Oct 19, 2023
ee046b1
Merge branch 'development' into update_urca
zingale Oct 20, 2023
1983da1
Merge branch 'development' into update_urca
zingale Oct 22, 2023
43779f1
another udpate
zingale Oct 22, 2023
9a52d1e
update network
zingale Oct 24, 2023
8eaf6ad
Merge branch 'development' into update_urca
zingale Oct 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 11 additions & 9 deletions networks/ignition_reaclib/URCA-simple/actual_network.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,19 +28,20 @@ namespace Rates
k_c12_c12_to_n_mg23 = 2,
k_c12_c12_to_p_na23 = 3,
k_he4_c12_to_o16 = 4,
k_n_to_p_weak_wc12 = 5,
k_na23_to_ne23 = 6,
k_ne23_to_na23 = 7,
NumRates = k_ne23_to_na23
k_na23_to_ne23 = 5,
k_ne23_to_na23 = 6,
k_n_to_p = 7,
k_p_to_n = 8,
NumRates = k_p_to_n
};

// number of reaclib rates

const int NrateReaclib = 5;
const int NrateReaclib = 4;

// number of tabular rates

const int NrateTabular = 2;
const int NrateTabular = 4;

// rate names -- note: the rates are 1-based, not zero-based, so we pad
// this vector with rate_names[0] = "" so the indices line up with the
Expand All @@ -52,9 +53,10 @@ namespace Rates
"c12_c12_to_n_mg23", // 2,
"c12_c12_to_p_na23", // 3,
"he4_c12_to_o16", // 4,
"n_to_p_weak_wc12", // 5,
"na23_to_ne23", // 6,
"ne23_to_na23" // 7,
"na23_to_ne23", // 5,
"ne23_to_na23", // 6,
"n_to_p", // 7,
"p_to_n" // 8,
};

}
Expand Down
5 changes: 3 additions & 2 deletions networks/ignition_reaclib/URCA-simple/actual_network_data.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,10 @@ namespace NSE_INDEX
-1, 3, 3, -1, 0, 8, -1,
-1, 3, 3, -1, 1, 7, -1,
-1, 2, 3, -1, -1, 4, -1,
-1, -1, 0, -1, -1, 1, -1,
-1, -1, 7, -1, -1, 6, -1,
-1, -1, 6, -1, -1, 7, 6
-1, -1, 6, -1, -1, 7, 5,
-1, -1, 0, -1, -1, 1, -1,
-1, -1, 1, -1, -1, 0, -1
};
}
#endif
Expand Down
45 changes: 36 additions & 9 deletions networks/ignition_reaclib/URCA-simple/actual_rhs.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,21 +133,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_na23_ne23_meta, j_na23_ne23_rhoy, j_na23_ne23_temp, j_na23_ne23_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_na23_to_ne23) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_na23_to_ne23) = drate_dt;
}
rate_eval.add_energy_rate(k_na23_to_ne23) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Na23) * (edot_nu + edot_gamma);

tabular_evaluate(j_ne23_na23_meta, j_ne23_na23_rhoy, j_ne23_na23_temp, j_ne23_na23_data,
rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
rate_eval.screened_rates(k_ne23_to_na23) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_ne23_to_na23) = drate_dt;
}
rate_eval.add_energy_rate(k_ne23_to_na23) = edot_nu + edot_gamma;
rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne23) * (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.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.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);


}
Expand All @@ -161,12 +179,14 @@ void rhs_nuc(const burn_t& state,
using namespace Rates;

ydot_nuc(N) =
-screened_rates(k_n_to_p_weak_wc12)*Y(N) +
-screened_rates(k_n_to_p)*Y(N) +
screened_rates(k_p_to_n)*Y(H1) +
0.5*screened_rates(k_c12_c12_to_n_mg23)*std::pow(Y(C12), 2)*state.rho;

ydot_nuc(H1) =
0.5*screened_rates(k_c12_c12_to_p_na23)*std::pow(Y(C12), 2)*state.rho +
screened_rates(k_n_to_p_weak_wc12)*Y(N);
screened_rates(k_n_to_p)*Y(N) +
-screened_rates(k_p_to_n)*Y(H1);

ydot_nuc(He4) =
0.5*screened_rates(k_c12_c12_to_he4_ne20)*std::pow(Y(C12), 2)*state.rho +
Expand Down Expand Up @@ -216,6 +236,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 @@ -225,9 +246,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(Na23) * rate_eval.add_energy_rate(k_na23_to_ne23);
enuc += C::Legacy::n_A * Y(Ne23) * rate_eval.add_energy_rate(k_ne23_to_na23);
// include any weak rate neutrino losses
enuc += rate_eval.enuc_weak;

// Get the thermal neutrino losses

Expand All @@ -252,15 +272,21 @@ void jac_nuc(const burn_t& state,

Real scratch;

scratch = -screened_rates(k_n_to_p_weak_wc12);
scratch = -screened_rates(k_n_to_p);
jac.set(N, N, scratch);

scratch = screened_rates(k_p_to_n);
jac.set(N, H1, scratch);

scratch = 1.0*screened_rates(k_c12_c12_to_n_mg23)*Y(C12)*state.rho;
jac.set(N, C12, scratch);

scratch = screened_rates(k_n_to_p_weak_wc12);
scratch = screened_rates(k_n_to_p);
jac.set(H1, N, scratch);

scratch = -screened_rates(k_p_to_n);
jac.set(H1, H1, scratch);

scratch = 1.0*screened_rates(k_c12_c12_to_p_na23)*Y(C12)*state.rho;
jac.set(H1, C12, scratch);

Expand Down Expand Up @@ -325,6 +351,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