update ECSN with latest pynucastro (#1320)

this will help us with the diffs in the pynucastro C++ CI

networks/ECSN/actual_network.H

@@ -93,7 +93,6 @@ namespace NSE_INDEX
     // last index is the corresponding reverse rate index.
 
     extern AMREX_GPU_MANAGED amrex::Array2D<int, 1, Rates::NumRates, 1, 7, Order::C> rate_indices;
-    extern AMREX_GPU_MANAGED bool initialized;
 }
 #endif
 

networks/ECSN/actual_network_data.cpp

@@ -30,7 +30,6 @@ namespace NSE_INDEX
         -1, -1, 3, -1, -1, 4, 15,
         -1, -1, 4, -1, -1, 5, 16
     };
-    AMREX_GPU_MANAGED bool initialized = false;
 }
 #endif
 

networks/ECSN/actual_rhs.H

@@ -70,7 +70,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
 
     Real ratraw, dratraw_dT;
     Real scor, dscor_dt;
-    Real scor2, dscor2_dt;
+    [[maybe_unused]] Real scor2, dscor2_dt;
 
 
     {
@@ -411,9 +411,9 @@ void actual_rhs (burn_t& state, Array1D<Real, 1, neqs>& ydot)
 
     // Get the thermal neutrino losses
 
-    Real sneut, dsneutdt, dsneutdd, snuda, snudz;
+    Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
 
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     // Append the energy equation (this is erg/g/s)
 
@@ -628,11 +628,11 @@ void actual_jac(const burn_t& state, MatrixType& jac)
 
     // Account for the thermal neutrino losses
 
-    Real sneut, dsneutdt, dsneutdd, snuda, snudz;
-    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, snuda, snudz);
+    Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    sneut5(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
-       Real b1 = (-state.abar * state.abar * snuda + (zion[j-1] - state.zbar) * state.abar * snudz);
+       Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);
        jac.add(net_ienuc, j, -b1);
     }
 

networks/ECSN/partition_functions.H

@@ -22,18 +22,25 @@ namespace part_fun {
     void interpolate_pf(const Real t9, const Real (&temp_array)[npts], const Real (&pf_array)[npts],
                         Real& pf, Real& dpf_dT) {
 
-        // find the index of the first temperature element we are larger than
+        if (t9 >= temp_array[0] && t9 < temp_array[npts-1]) {
 
-        int idx = -1;
+            // find the largest temperature element <= t9 using a binary search
 
-        for (int i = 0; i < npts-1; ++i) {
-            if (t9 >= temp_array[i] && t9 < temp_array[i+1]) {
-                idx = i;
-                break;
+            int left = 0;
+            int right = npts;
+
+            while (left < right) {
+                int mid = (left + right) / 2;
+                if (temp_array[mid] > t9) {
+                    right = mid;
+                } else {
+                    left = mid + 1;
+                }
             }
-        }
 
-        if (idx >= 0) {
+            const int idx = right - 1;
+
+            // now we have temp_array[idx] <= t9 < temp_array[idx+1]
 
             // construct the slope -- this is (log10(pf_{i+1}) - log10(pf_i)) / (T_{i+1} - T_i)
 
@@ -46,12 +53,12 @@ namespace part_fun {
 
             // find the derivative (with respect to T, not T9)
 
-            Real dpf_dT9 = pf * std::log(10.0_rt) * slope;
+            Real dpf_dT9 = pf * M_LN10 * slope;
             dpf_dT = dpf_dT9 / 1.e9_rt;
 
         } else {
 
-            // T < the smallest T in the partition function table
+            // T < the smallest T or >= the largest T in the partition function table
             pf = 1.0;
             dpf_dT = 0.0;
 
@@ -90,48 +97,27 @@ constexpr Real get_spin_state(const int inuc) {
 
     switch (inuc) {
 
-    case H1:
-        spin = 2;
-        break;
-
     case He4:
-        spin = 1;
-        break;
-
     case O16:
-        spin = 1;
-        break;
-
     case O20:
-        spin = 1;
-        break;
-
-    case F20:
-        spin = 5;
-        break;
-
     case Ne20:
-        spin = 1;
-        break;
-
     case Mg24:
-        spin = 1;
-        break;
-
-    case Al27:
-        spin = 6;
-        break;
-
     case Si28:
+    case S32:
         spin = 1;
         break;
 
+    case H1:
     case P31:
         spin = 2;
         break;
 
-    case S32:
-        spin = 1;
+    case F20:
+        spin = 5;
+        break;
+
+    case Al27:
+        spin = 6;
         break;
 
 

networks/ECSN/reaclib_rates.H

@@ -910,11 +910,11 @@ fill_reaclib_rates(const tf_t& tfactors, T& rate_eval)
 template <int do_T_derivatives, typename T>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void
-fill_approx_rates([[maybe_unused]] const tf_t& tfactors, T& rate_eval)
+fill_approx_rates([[maybe_unused]] const tf_t& tfactors, [[maybe_unused]] T& rate_eval)
 {
 
-    Real rate;
-    Real drate_dT;
+    [[maybe_unused]] Real rate{};
+    [[maybe_unused]] Real drate_dT{};
 
 
 }

networks/ECSN/table_rates.H

@@ -112,6 +112,10 @@ void init_tab_info(const table_t& tf, const std::string& file, R& rhoy, T& temp,
             }
 
             if (1 == j) {
+                if (temp(i) > 33.0_rt){
+                    amrex::Error("init_tab_info: Table data needs to be updated");
+                }
+
                 temp(i) = std::pow(10.0_rt, temp(i));
             }
 
@@ -151,11 +155,11 @@ int vector_index_lu(const int vlen, const V& vector, const Real fvar)
             } else {
                 ndn = j;
             }
-          if (((nup - ndn) == 1)) {
-              index = ndn;
-              return index;
-          }
+            if ((nup - ndn) == 1) {
+                break;
+            }
         }
+        index = ndn;
     }
 
     return index;
@@ -256,6 +260,7 @@ get_entries(const table_t& table_meta,
                                data(itemp_hi, irhoy_hi, ivar),
                                temp);
 
+        // NOLINTNEXTLINE(readability-suspicious-call-argument)
         entries(ivar) = bl_extrap(rhoy_lo, rhoy_hi, f_i, f_ip1, rhoy);
     }