Skip to content

Commit

Permalink
update the Urca network with the p -> n electron capture rate from La…
Browse files Browse the repository at this point in the history 
…nganke
  • Loading branch information
@zingale
zingale committed Oct 18, 2023
1 parent 8b2c485 commit b6ae2c3
Show file tree
Hide file tree
Showing 7 changed files with 94 additions and 62 deletions.
20 changes: 11 additions & 9 deletions networks/ignition_reaclib/URCA-simple/actual_network.H
View file
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
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
34 changes: 30 additions & 4 deletions networks/ignition_reaclib/URCA-simple/actual_rhs.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,22 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
}
rate_eval.add_energy_rate(k_ne23_to_na23) = 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;

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;


}

Expand All @@ -161,12 +177,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 @@ -228,6 +246,8 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
// 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);
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);

// 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
37 changes: 0 additions & 37 deletions networks/ignition_reaclib/URCA-simple/reaclib_rates.H
View file
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
View file
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
View file
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
View file
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()

0 comments on commit b6ae2c3

Please sign in to comment.