From 89a1b8ddb4afdfcc4c0d8191f160b9b7f2be87bd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 23 Oct 2023 13:56:03 -0400 Subject: [PATCH 1/2] update ECSN to latest pynucastro --- networks/ECSN/actual_rhs.H | 19 ++++++++++--------- networks/ECSN/reaclib_rates.H | 4 ++-- networks/ECSN/table_rates.H | 13 +++++++++++-- unit_test/test_rhs/ci-benchmarks/ecsn.out | 2 +- 4 files changed, 24 insertions(+), 14 deletions(-) diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H index 640781ea89..99a2fc4ae8 100644 --- a/networks/ECSN/actual_rhs.H +++ b/networks/ECSN/actual_rhs.H @@ -268,13 +268,15 @@ 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_f20_o20_meta, j_f20_o20_rhoy, j_f20_o20_temp, j_f20_o20_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); rate_eval.screened_rates(k_f20_to_o20) = rate; if constexpr (std::is_same::value) { rate_eval.dscreened_rates_dT(k_f20_to_o20) = drate_dt; } - rate_eval.add_energy_rate(k_f20_to_o20) = edot_nu + edot_gamma; + rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma); tabular_evaluate(j_ne20_f20_meta, j_ne20_f20_rhoy, j_ne20_f20_temp, j_ne20_f20_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -282,7 +284,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { if constexpr (std::is_same::value) { rate_eval.dscreened_rates_dT(k_ne20_to_f20) = drate_dt; } - rate_eval.add_energy_rate(k_ne20_to_f20) = edot_nu + edot_gamma; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne20) * (edot_nu + edot_gamma); tabular_evaluate(j_o20_f20_meta, j_o20_f20_rhoy, j_o20_f20_temp, j_o20_f20_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -290,7 +292,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { if constexpr (std::is_same::value) { rate_eval.dscreened_rates_dT(k_o20_to_f20) = drate_dt; } - rate_eval.add_energy_rate(k_o20_to_f20) = edot_nu + edot_gamma; + rate_eval.enuc_weak += C::Legacy::n_A * Y(O20) * (edot_nu + edot_gamma); tabular_evaluate(j_f20_ne20_meta, j_f20_ne20_rhoy, j_f20_ne20_temp, j_f20_ne20_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -298,7 +300,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { if constexpr (std::is_same::value) { rate_eval.dscreened_rates_dT(k_f20_to_ne20) = drate_dt; } - rate_eval.add_energy_rate(k_f20_to_ne20) = edot_nu + edot_gamma; + rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma); } @@ -394,6 +396,7 @@ void actual_rhs (burn_t& state, Array1D& ydot) rate_t rate_eval; constexpr int do_T_derivatives = 0; + evaluate_rates(state, rate_eval); rhs_nuc(state, ydot, Y, rate_eval.screened_rates); @@ -403,11 +406,8 @@ void actual_rhs (burn_t& state, Array1D& ydot) Real enuc; ener_gener_rate(ydot, enuc); - // include reaction neutrino losses (non-thermal) and gamma heating - enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_o20); - enuc += C::Legacy::n_A * Y(Ne20) * rate_eval.add_energy_rate(k_ne20_to_f20); - enuc += C::Legacy::n_A * Y(O20) * rate_eval.add_energy_rate(k_o20_to_f20); - enuc += C::Legacy::n_A * Y(F20) * rate_eval.add_energy_rate(k_f20_to_ne20); + // include any weak rate neutrino losses + enuc += rate_eval.enuc_weak; // Get the thermal neutrino losses @@ -613,6 +613,7 @@ void actual_jac(const burn_t& state, MatrixType& jac) rate_derivs_t rate_eval; constexpr int do_T_derivatives = 1; + evaluate_rates(state, rate_eval); // Species Jacobian elements with respect to other species diff --git a/networks/ECSN/reaclib_rates.H b/networks/ECSN/reaclib_rates.H index efbfe0e818..723ae25a37 100644 --- a/networks/ECSN/reaclib_rates.H +++ b/networks/ECSN/reaclib_rates.H @@ -13,13 +13,13 @@ using namespace Species; struct rate_t { Array1D screened_rates; - Array1D add_energy_rate; + Real enuc_weak; }; struct rate_derivs_t { Array1D screened_rates; Array1D dscreened_rates_dT; - Array1D add_energy_rate; + Real enuc_weak; }; diff --git a/networks/ECSN/table_rates.H b/networks/ECSN/table_rates.H index 9c5946e9fe..64d1413cc3 100644 --- a/networks/ECSN/table_rates.H +++ b/networks/ECSN/table_rates.H @@ -86,6 +86,11 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table std::ifstream table; table.open(file); + if (!table.is_open()) { + // the table was not present or we could not open it; abort + amrex::Error("table could not be opened"); + } + std::string line; // read and skip over the header @@ -99,6 +104,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& log_rhoy_table for (int j = 1; j <= tf.nrhoy; ++j) { for (int i = 1; i <= tf.ntemp; ++i) { std::getline(table, line); + if (line.empty()) { + amrex::Error("Error reading table data"); + } + std::istringstream sdata(line); sdata >> log_rhoy_table(j) >> log_temp_table(i); @@ -265,8 +274,8 @@ evaluate_dr_dtemp(const table_t& table_meta, const R& log_rhoy_table, const T& l // Finally, we perform a 1d-linear interpolation between dlogr_dlogt_ip1 // and dlogr_dlogt_i to compute dlogr_dlogt - Real log_rhoy_lo = log_temp_table(irhoy_lo); - Real log_rhoy_hi = log_temp_table(irhoy_hi); + Real log_rhoy_lo = log_rhoy_table(irhoy_lo); + Real log_rhoy_hi = log_rhoy_table(irhoy_hi); Real log_temp_lo = log_temp_table(jtemp_lo); Real log_temp_hi = log_temp_table(jtemp_hi); diff --git a/unit_test/test_rhs/ci-benchmarks/ecsn.out b/unit_test/test_rhs/ci-benchmarks/ecsn.out index 4b7b5e98d0..9489ba5c4d 100644 --- a/unit_test/test_rhs/ci-benchmarks/ecsn.out +++ b/unit_test/test_rhs/ci-benchmarks/ecsn.out @@ -161,7 +161,7 @@ J_hydrogen-1_E -0.12580793235 1.2089803642e+12 J_helium-4_E -0.0078065301852 1.4905283338e+12 J_oxygen-16_E -5.3990173959e+12 1.2739583515e-06 - J_oxygen-20_E -2.079874023e-18 3.5544574816e-18 + J_oxygen-20_E -2.8880340474e-19 3.5544574816e-18 J_fluorine-20_E -2.0716943987e-18 7.6531304664e-18 J_neon-20_E -8.8836179369e-05 1.9862184663e-05 J_magnesium-24_E -0.0040169125096 0.018563573855 From e9d4be4f5d261c5f8efe767105ab175a363bc324 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 23 Oct 2023 14:04:46 -0400 Subject: [PATCH 2/2] update --- .../ci-benchmarks/ecsn_unit_test.out | 50 +++++++++---------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out index 1ae3f2bc19..bc26ca103b 100644 --- a/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out +++ b/unit_test/burn_cell/ci-benchmarks/ecsn_unit_test.out @@ -1,4 +1,4 @@ -AMReX (23.07-106-gde7c6189623b-dirty) initialized +AMReX (23.10-13-gecaa14456e3f) initialized starting the single zone burn... reading in network electron-capture / beta-decay tables... Maximum Time (s): 0.01 @@ -29,37 +29,37 @@ RHS at t = 0 s32 1323519.174 ------------------------------------ successful? 1 - - Hnuc = 4.716248406e+19 - - added e = 4.716248406e+17 - - final T = 6804058463 + - Hnuc = 4.4691188e+19 + - added e = 4.4691188e+17 + - final T = 6711517341 ------------------------------------ e initial = 5.995956082e+17 -e final = 1.071220449e+18 +e final = 1.046507488e+18 ------------------------------------ new mass fractions: -h1 0.009032124504 -he4 4.003293528e-13 -o16 5.888714278e-07 -o20 0.00985678636 -f20 0.009855432543 -ne20 7.501479981e-14 -mg24 1.21710929e-09 -al27 1.951815457e-21 -si28 0.26794899 -p31 8.245884489e-14 -s32 0.7033060765 +h1 0.007036728874 +he4 4.569221711e-13 +o16 6.11671909e-07 +o20 0.009363850335 +f20 0.009362566168 +ne20 8.595667728e-14 +mg24 8.576669062e-06 +al27 9.999999996e-31 +si28 0.2286832772 +p31 9.999999996e-31 +s32 0.745544389 ------------------------------------ species creation rates: -omegadot(h1): -0.0967875496 +omegadot(h1): -0.2963271126 omegadot(he4): -3 -omegadot(o16): -49.99994111 -omegadot(o20): -0.01432136399 -omegadot(f20): -0.0144567457 +omegadot(o16): -49.99993883 +omegadot(o20): -0.0636149665 +omegadot(f20): -0.06374338321 omegadot(ne20): -30 -omegadot(mg24): -9.999999878 +omegadot(mg24): -9.999142333 omegadot(al27): -1 -omegadot(si28): 25.794899 +omegadot(si28): 21.86832772 omegadot(p31): -1 -omegadot(s32): 69.33060765 -number of steps taken: 31502 -AMReX (23.07-106-gde7c6189623b-dirty) finalized +omegadot(s32): 73.5544389 +number of steps taken: 14714 +AMReX (23.10-13-gecaa14456e3f) finalized