update ECSN to latest pynucastro (#1368)

AMReX-Astro · Oct 23, 2023 · e8749f9 · e8749f9
1 parent 8ec4f08
commit e8749f9
Show file tree

Hide file tree

Showing 5 changed files with 49 additions and 39 deletions.
diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H
@@ -268,37 +268,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_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<T, rate_derivs_t>::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);
     rate_eval.screened_rates(k_ne20_to_f20) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::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);
     rate_eval.screened_rates(k_o20_to_f20) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::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);
     rate_eval.screened_rates(k_f20_to_ne20) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::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<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);
@@ -403,11 +406,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(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<do_T_derivatives, rate_derivs_t>(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
@@ -13,13 +13,13 @@ using namespace Species;
 
 struct rate_t {
     Array1D<Real, 1, NumRates>  screened_rates;
-    Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
+    Real enuc_weak;
 };
 
 struct rate_derivs_t {
     Array1D<Real, 1, NumRates>  screened_rates;
     Array1D<Real, 1, NumRates>  dscreened_rates_dT;
-    Array1D<Real, NrateReaclib+1, NrateReaclib+NrateTabular> add_energy_rate;
+    Real enuc_weak;
 };
 
 

diff --git 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/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
diff --git 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