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.

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
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
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
32 changes: 28 additions & 4 deletions networks/ignition_reaclib/URCA-simple/actual_rhs.H
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,22 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
}
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 @@ -163,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 @@ -254,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
148 changes: 148 additions & 0 deletions networks/ignition_reaclib/URCA-simple/n-p_betadecay.dat

Large diffs are not rendered by default.

148 changes: 148 additions & 0 deletions networks/ignition_reaclib/URCA-simple/p-n_electroncapture.dat

Large diffs are not rendered by default.

37 changes: 0 additions & 37 deletions networks/ignition_reaclib/URCA-simple/reaclib_rates.H
Original file line number Diff line number Diff line change
Expand Up @@ -168,37 +168,6 @@ void rate_He4_C12_to_O16(const tf_t& tfactors, Real& rate, Real& drate_dT) {

}

template <int do_T_derivatives>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void rate_n_to_p_weak_wc12(const tf_t& tfactors, Real& rate, Real& drate_dT) {

// n --> p

rate = 0.0;
drate_dT = 0.0;

Real ln_set_rate{0.0};
Real dln_set_rate_dT9{0.0};
Real set_rate{0.0};

// wc12w
ln_set_rate = -6.78161;
amrex::ignore_unused(tfactors);

if constexpr (do_T_derivatives) {
dln_set_rate_dT9 = 0.0;
}

// avoid underflows by zeroing rates in [0.0, 1.e-100]
ln_set_rate = std::max(ln_set_rate, -230.0);
set_rate = std::exp(ln_set_rate);
rate += set_rate;
if constexpr (do_T_derivatives) {
drate_dT += set_rate * dln_set_rate_dT9 / 1.0e9;
}

}



template <int do_T_derivatives, typename T>
Expand Down Expand Up @@ -234,12 +203,6 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval)
rate_eval.dscreened_rates_dT(k_He4_C12_to_O16) = drate_dT;

}
rate_n_to_p_weak_wc12<do_T_derivatives>(tfactors, rate, drate_dT);
rate_eval.screened_rates(k_n_to_p_weak_wc12) = rate;
if constexpr (std::is_same<T, rate_derivs_t>::value) {
rate_eval.dscreened_rates_dT(k_n_to_p_weak_wc12) = drate_dT;

}

}

Expand Down
12 changes: 11 additions & 1 deletion networks/ignition_reaclib/URCA-simple/table_rates.H
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ void init_tabular();
// Log(g/cm^3) Log(K) erg erg erg Log(1/s) Log(erg/s) Log(erg/s)
//

const int num_tables = 2;
const int num_tables = 4;

enum TableVars
{
Expand Down Expand Up @@ -65,6 +65,16 @@ namespace rate_tables
extern AMREX_GPU_MANAGED Array1D<Real, 1, 152> j_Ne23_Na23_rhoy;
extern AMREX_GPU_MANAGED Array1D<Real, 1, 39> j_Ne23_Na23_temp;

extern AMREX_GPU_MANAGED table_t j_n_p_meta;
extern AMREX_GPU_MANAGED Array3D<Real, 1, 13, 1, 11, 1, 6> j_n_p_data;
extern AMREX_GPU_MANAGED Array1D<Real, 1, 11> j_n_p_rhoy;
extern AMREX_GPU_MANAGED Array1D<Real, 1, 13> j_n_p_temp;

extern AMREX_GPU_MANAGED table_t j_p_n_meta;
extern AMREX_GPU_MANAGED Array3D<Real, 1, 13, 1, 11, 1, 6> j_p_n_data;
extern AMREX_GPU_MANAGED Array1D<Real, 1, 11> j_p_n_rhoy;
extern AMREX_GPU_MANAGED Array1D<Real, 1, 13> j_p_n_temp;

}

template <typename R, typename T, typename D>
Expand Down
26 changes: 26 additions & 0 deletions networks/ignition_reaclib/URCA-simple/table_rates_data.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,16 @@ namespace rate_tables
AMREX_GPU_MANAGED Array1D<Real, 1, 152> j_Ne23_Na23_rhoy;
AMREX_GPU_MANAGED Array1D<Real, 1, 39> j_Ne23_Na23_temp;

AMREX_GPU_MANAGED table_t j_n_p_meta;
AMREX_GPU_MANAGED Array3D<Real, 1, 13, 1, 11, 1, 6> j_n_p_data;
AMREX_GPU_MANAGED Array1D<Real, 1, 11> j_n_p_rhoy;
AMREX_GPU_MANAGED Array1D<Real, 1, 13> j_n_p_temp;

AMREX_GPU_MANAGED table_t j_p_n_meta;
AMREX_GPU_MANAGED Array3D<Real, 1, 13, 1, 11, 1, 6> j_p_n_data;
AMREX_GPU_MANAGED Array1D<Real, 1, 11> j_p_n_rhoy;
AMREX_GPU_MANAGED Array1D<Real, 1, 13> j_p_n_temp;


}

Expand Down Expand Up @@ -45,5 +55,21 @@ void init_tabular()
init_tab_info(j_Ne23_Na23_meta, "23ne-23na_betadecay.dat", j_Ne23_Na23_rhoy, j_Ne23_Na23_temp, j_Ne23_Na23_data);


j_n_p_meta.ntemp = 13;
j_n_p_meta.nrhoy = 11;
j_n_p_meta.nvars = 6;
j_n_p_meta.nheader = 5;

init_tab_info(j_n_p_meta, "n-p_betadecay.dat", j_n_p_rhoy, j_n_p_temp, j_n_p_data);


j_p_n_meta.ntemp = 13;
j_p_n_meta.nrhoy = 11;
j_p_n_meta.nvars = 6;
j_p_n_meta.nheader = 5;

init_tab_info(j_p_n_meta, "p-n_electroncapture.dat", j_p_n_rhoy, j_p_n_temp, j_p_n_data);



}
22 changes: 13 additions & 9 deletions networks/ignition_reaclib/URCA-simple/urca.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
# C-burning with A=23 URCA rate module generator

from pynucastro.networks import AmrexAstroCxxNetwork
import pynucastro as pyna

files = ["c12-c12a-ne20-cf88",
"c12-c12n-mg23-cf88",
"c12-c12p-na23-cf88",
"c12-ag-o16-nac2",
"na23--ne23-toki",
"ne23--na23-toki",
"n--p-wc12"]
rl = pyna.ReacLibLibrary()
rl_rates = rl.get_rate_by_name(["c12(c12,a)ne20",
"c12(c12,n)mg23",
"c12(c12,p)na23",
"c12(a,g)o16"])

urca_net = AmrexAstroCxxNetwork(files)
tl = pyna.TabularLibrary()
tl_rates = tl.get_rate_by_name(["na23(,)ne23",
"ne23(,)na23",
"n(,)p",
"p(,)n"])

urca_net = pyna.AmrexAstroCxxNetwork(rates=rl_rates+tl_rates)
urca_net.write_network()