From b6ae2c3256d050f40c8e64001c6fda0eab2793ea Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 18 Oct 2023 14:36:56 -0400 Subject: [PATCH] update the Urca network with the p -> n electron capture rate from Langanke --- .../URCA-simple/actual_network.H | 20 +++++----- .../URCA-simple/actual_network_data.cpp | 5 ++- .../ignition_reaclib/URCA-simple/actual_rhs.H | 34 +++++++++++++++-- .../URCA-simple/reaclib_rates.H | 37 ------------------- .../URCA-simple/table_rates.H | 12 +++++- .../URCA-simple/table_rates_data.cpp | 26 +++++++++++++ networks/ignition_reaclib/URCA-simple/urca.py | 22 ++++++----- 7 files changed, 94 insertions(+), 62 deletions(-) diff --git a/networks/ignition_reaclib/URCA-simple/actual_network.H b/networks/ignition_reaclib/URCA-simple/actual_network.H index 51cbc00c79..15b00fb144 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_network.H +++ b/networks/ignition_reaclib/URCA-simple/actual_network.H @@ -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 @@ -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, }; } diff --git a/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp b/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp index 4a7e7368ae..39fcc09e09 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp +++ b/networks/ignition_reaclib/URCA-simple/actual_network_data.cpp @@ -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 diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H index d020349fb6..21af576863 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H @@ -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::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::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; + } @@ -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 + @@ -228,6 +246,8 @@ void actual_rhs (burn_t& state, Array1D& 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 @@ -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); diff --git a/networks/ignition_reaclib/URCA-simple/reaclib_rates.H b/networks/ignition_reaclib/URCA-simple/reaclib_rates.H index 8b6ddc7e39..00ceb89957 100644 --- a/networks/ignition_reaclib/URCA-simple/reaclib_rates.H +++ b/networks/ignition_reaclib/URCA-simple/reaclib_rates.H @@ -168,37 +168,6 @@ void rate_he4_c12_to_o16(const tf_t& tfactors, Real& rate, Real& drate_dT) { } -template -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 @@ -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(tfactors, rate, drate_dT); - rate_eval.screened_rates(k_n_to_p_weak_wc12) = rate; - if constexpr (std::is_same::value) { - rate_eval.dscreened_rates_dT(k_n_to_p_weak_wc12) = drate_dT; - - } } diff --git a/networks/ignition_reaclib/URCA-simple/table_rates.H b/networks/ignition_reaclib/URCA-simple/table_rates.H index 6f4af69558..b2936192db 100644 --- a/networks/ignition_reaclib/URCA-simple/table_rates.H +++ b/networks/ignition_reaclib/URCA-simple/table_rates.H @@ -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 { @@ -65,6 +65,16 @@ namespace rate_tables extern AMREX_GPU_MANAGED Array1D j_ne23_na23_rhoy; extern AMREX_GPU_MANAGED Array1D j_ne23_na23_temp; + extern AMREX_GPU_MANAGED table_t j_n_p_meta; + extern AMREX_GPU_MANAGED Array3D j_n_p_data; + extern AMREX_GPU_MANAGED Array1D j_n_p_rhoy; + extern AMREX_GPU_MANAGED Array1D j_n_p_temp; + + extern AMREX_GPU_MANAGED table_t j_p_n_meta; + extern AMREX_GPU_MANAGED Array3D j_p_n_data; + extern AMREX_GPU_MANAGED Array1D j_p_n_rhoy; + extern AMREX_GPU_MANAGED Array1D j_p_n_temp; + } template diff --git a/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp b/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp index 61b399cf4b..4719bd77c9 100644 --- a/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp +++ b/networks/ignition_reaclib/URCA-simple/table_rates_data.cpp @@ -18,6 +18,16 @@ namespace rate_tables AMREX_GPU_MANAGED Array1D j_ne23_na23_rhoy; AMREX_GPU_MANAGED Array1D j_ne23_na23_temp; + AMREX_GPU_MANAGED table_t j_n_p_meta; + AMREX_GPU_MANAGED Array3D j_n_p_data; + AMREX_GPU_MANAGED Array1D j_n_p_rhoy; + AMREX_GPU_MANAGED Array1D j_n_p_temp; + + AMREX_GPU_MANAGED table_t j_p_n_meta; + AMREX_GPU_MANAGED Array3D j_p_n_data; + AMREX_GPU_MANAGED Array1D j_p_n_rhoy; + AMREX_GPU_MANAGED Array1D j_p_n_temp; + } @@ -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); + + } diff --git a/networks/ignition_reaclib/URCA-simple/urca.py b/networks/ignition_reaclib/URCA-simple/urca.py index 5db0666b34..4c8e1c8288 100644 --- a/networks/ignition_reaclib/URCA-simple/urca.py +++ b/networks/ignition_reaclib/URCA-simple/urca.py @@ -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()