another udpate

AMReX-Astro · Oct 22, 2023 · 43779f1 · 43779f1
1 parent 1983da1
commit 43779f1
Show file tree

Hide file tree

Showing 3 changed files with 23 additions and 13 deletions.
diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H
@@ -133,37 +133,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_na23_ne23_meta, j_na23_ne23_rhoy, j_na23_ne23_temp, j_na23_ne23_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
     rate_eval.screened_rates(k_na23_to_ne23) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_na23_to_ne23) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_na23_to_ne23) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(Na23) * (edot_nu + edot_gamma);
 
     tabular_evaluate(j_ne23_na23_meta, j_ne23_na23_rhoy, j_ne23_na23_temp, j_ne23_na23_data,
                      rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma);
     rate_eval.screened_rates(k_ne23_to_na23) = rate;
     if constexpr (std::is_same<T, rate_derivs_t>::value) {
         rate_eval.dscreened_rates_dT(k_ne23_to_na23) = drate_dt;
     }
-    rate_eval.add_energy_rate(k_ne23_to_na23) = edot_nu + edot_gamma;
+    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.add_energy_rate(k_n_to_p) = edot_nu + edot_gamma;
+    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.add_energy_rate(k_p_to_n) = edot_nu + edot_gamma;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma);
 
 
 }
@@ -234,6 +236,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);
@@ -243,11 +246,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(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);
+    // include any weak rate neutrino losses
+    enuc += rate_eval.enuc_weak;
 
     // Get the thermal neutrino losses
 
@@ -351,6 +351,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/ignition_reaclib/URCA-simple/reaclib_rates.H b/networks/ignition_reaclib/URCA-simple/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/ignition_reaclib/URCA-simple/table_rates.H b/networks/ignition_reaclib/URCA-simple/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);