[DON'T MERGE] test the constexpr ener_gener_rate

networks/CNO_extras/Make.package

@@ -4,6 +4,7 @@ ifeq ($(USE_REACT),TRUE)
   CEXE_sources += actual_network_data.cpp
   CEXE_headers += actual_network.H
   CEXE_headers += tfactors.H
+  CEXE_headers += interp_tools.H
   CEXE_headers += partition_functions.H
   CEXE_sources += partition_functions_data.cpp
   CEXE_headers += actual_rhs.H

networks/CNO_extras/actual_network.H

@@ -285,14 +285,6 @@ namespace Rates
         NumRates = k_He4_He4_He4_to_C12
     };
 
-    // number of reaclib rates
-
-    const int NrateReaclib = 67;
-
-    // number of tabular rates
-
-    const int NrateTabular = 0;
-
     // 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
     // NetworkRates enum
@@ -381,7 +373,7 @@ namespace NSE_INDEX
     // First 3 row indices for reactants, followed by 3 product indices
     // last index is the corresponding reverse rate index.
 
-    extern AMREX_GPU_MANAGED amrex::Array2D<int, 1, Rates::NumRates, 1, 7, amrex::Order::C> rate_indices;
+    extern AMREX_GPU_MANAGED amrex::Array2D<std::int8_t, 1, Rates::NumRates, 1, 7, amrex::Order::C> rate_indices;
 }
 #endif
 

networks/CNO_extras/actual_network_data.cpp

@@ -4,74 +4,74 @@
 #ifdef NSE_NET
 namespace NSE_INDEX
 {
-    AMREX_GPU_MANAGED amrex::Array2D<int, 1, Rates::NumRates, 1, 7, amrex::Order::C> rate_indices {
-        -1, -1, 4, -1, -1, 3, -1,
-        -1, -1, 7, -1, -1, 5, -1,
-        -1, -1, 8, -1, -1, 6, -1,
-        -1, -1, 12, -1, -1, 10, -1,
-        -1, -1, 13, -1, -1, 11, -1,
-        -1, -1, 15, -1, -1, 13, -1,
-        -1, -1, 16, -1, -1, 14, -1,
-        -1, -1, 4, -1, 0, 2, -1,
-        -1, -1, 5, -1, 0, 3, -1,
-        -1, -1, 7, -1, 0, 4, -1,
-        -1, -1, 8, -1, 0, 5, -1,
-        -1, -1, 9, -1, 0, 6, -1,
-        -1, -1, 9, -1, 1, 2, -1,
-        -1, -1, 12, -1, 0, 9, -1,
-        -1, -1, 13, -1, 0, 10, -1,
-        -1, -1, 13, -1, 1, 5, -1,
-        -1, -1, 14, -1, 0, 11, -1,
-        -1, -1, 14, -1, 1, 6, -1,
-        -1, -1, 15, -1, 0, 12, -1,
-        -1, -1, 15, -1, 1, 7, -1,
-        -1, -1, 16, -1, 0, 13, -1,
-        -1, -1, 16, -1, 1, 8, -1,
-        -1, -1, 17, -1, 0, 14, -1,
-        -1, -1, 17, -1, 1, 9, -1,
-        -1, -1, 18, -1, 1, 15, -1,
-        -1, -1, 19, -1, 1, 17, -1,
-        -1, -1, 2, 1, 1, 1, -1,
-        -1, 0, 2, -1, -1, 4, 8,
-        -1, 1, 2, -1, -1, 9, 13,
-        -1, 0, 3, -1, -1, 5, 9,
-        -1, 0, 4, -1, -1, 7, 10,
-        -1, 0, 5, -1, -1, 8, 11,
-        -1, 1, 5, -1, -1, 13, 16,
-        -1, 0, 6, -1, -1, 9, 12,
-        -1, 1, 6, -1, -1, 14, 18,
-        -1, 1, 7, -1, -1, 15, 20,
-        -1, 1, 8, -1, -1, 16, 22,
-        -1, 0, 9, -1, -1, 12, 14,
-        -1, 1, 9, -1, -1, 17, 24,
-        -1, 0, 10, -1, -1, 13, 15,
-        -1, 0, 11, -1, -1, 14, 17,
-        -1, 0, 12, -1, -1, 15, 19,
-        -1, 0, 13, -1, -1, 16, 21,
-        -1, 0, 14, -1, -1, 17, 23,
-        -1, 1, 15, -1, -1, 18, 25,
-        -1, 1, 17, -1, -1, 19, 26,
-        -1, 1, 2, -1, 0, 6, -1,
-        -1, 2, 2, -1, 1, 17, 65,
-        -1, 1, 4, -1, 0, 9, 55,
-        -1, 1, 5, -1, 0, 10, -1,
-        -1, 0, 6, -1, 1, 2, 47,
-        -1, 1, 6, -1, 0, 11, -1,
-        -1, 1, 7, -1, 0, 12, 60,
-        -1, 1, 8, -1, 0, 13, -1,
-        -1, 0, 9, -1, 1, 4, -1,
-        -1, 1, 9, -1, 0, 14, -1,
-        -1, 2, 9, -1, 1, 19, 66,
-        -1, 0, 10, -1, 1, 5, 50,
-        -1, 0, 11, -1, 1, 6, 52,
-        -1, 0, 12, -1, 1, 7, -1,
-        -1, 1, 12, -1, 0, 17, 64,
-        -1, 0, 13, -1, 1, 8, 54,
-        -1, 0, 14, -1, 1, 9, 56,
-        -1, 0, 17, -1, 1, 12, -1,
-        -1, 1, 17, -1, 2, 2, -1,
-        -1, 1, 19, -1, 2, 9, -1,
-        1, 1, 1, -1, -1, 2, 27
+    AMREX_GPU_MANAGED amrex::Array2D<std::int8_t, 1, Rates::NumRates, 1, 7, amrex::Order::C> rate_indices {
+        -1, -1, 4, -1, -1, 3, -1,  // N13_to_C13_weak_wc12
+        -1, -1, 7, -1, -1, 5, -1,  // O14_to_N14_weak_wc12
+        -1, -1, 8, -1, -1, 6, -1,  // O15_to_N15_weak_wc12
+        -1, -1, 12, -1, -1, 10, -1,  // F17_to_O17_weak_wc12
+        -1, -1, 13, -1, -1, 11, -1,  // F18_to_O18_weak_wc12
+        -1, -1, 15, -1, -1, 13, -1,  // Ne18_to_F18_weak_wc12
+        -1, -1, 16, -1, -1, 14, -1,  // Ne19_to_F19_weak_wc12
+        -1, -1, 4, -1, 0, 2, -1,  // N13_to_p_C12
+        -1, -1, 5, -1, 0, 3, -1,  // N14_to_p_C13
+        -1, -1, 7, -1, 0, 4, -1,  // O14_to_p_N13
+        -1, -1, 8, -1, 0, 5, -1,  // O15_to_p_N14
+        -1, -1, 9, -1, 0, 6, -1,  // O16_to_p_N15
+        -1, -1, 9, -1, 1, 2, -1,  // O16_to_He4_C12
+        -1, -1, 12, -1, 0, 9, -1,  // F17_to_p_O16
+        -1, -1, 13, -1, 0, 10, -1,  // F18_to_p_O17
+        -1, -1, 13, -1, 1, 5, -1,  // F18_to_He4_N14
+        -1, -1, 14, -1, 0, 11, -1,  // F19_to_p_O18
+        -1, -1, 14, -1, 1, 6, -1,  // F19_to_He4_N15
+        -1, -1, 15, -1, 0, 12, -1,  // Ne18_to_p_F17
+        -1, -1, 15, -1, 1, 7, -1,  // Ne18_to_He4_O14
+        -1, -1, 16, -1, 0, 13, -1,  // Ne19_to_p_F18
+        -1, -1, 16, -1, 1, 8, -1,  // Ne19_to_He4_O15
+        -1, -1, 17, -1, 0, 14, -1,  // Ne20_to_p_F19
+        -1, -1, 17, -1, 1, 9, -1,  // Ne20_to_He4_O16
+        -1, -1, 18, -1, 1, 15, -1,  // Mg22_to_He4_Ne18
+        -1, -1, 19, -1, 1, 17, -1,  // Mg24_to_He4_Ne20
+        -1, -1, 2, 1, 1, 1, -1,  // C12_to_He4_He4_He4
+        -1, 0, 2, -1, -1, 4, 8,  // p_C12_to_N13
+        -1, 1, 2, -1, -1, 9, 13,  // He4_C12_to_O16
+        -1, 0, 3, -1, -1, 5, 9,  // p_C13_to_N14
+        -1, 0, 4, -1, -1, 7, 10,  // p_N13_to_O14
+        -1, 0, 5, -1, -1, 8, 11,  // p_N14_to_O15
+        -1, 1, 5, -1, -1, 13, 16,  // He4_N14_to_F18
+        -1, 0, 6, -1, -1, 9, 12,  // p_N15_to_O16
+        -1, 1, 6, -1, -1, 14, 18,  // He4_N15_to_F19
+        -1, 1, 7, -1, -1, 15, 20,  // He4_O14_to_Ne18
+        -1, 1, 8, -1, -1, 16, 22,  // He4_O15_to_Ne19
+        -1, 0, 9, -1, -1, 12, 14,  // p_O16_to_F17
+        -1, 1, 9, -1, -1, 17, 24,  // He4_O16_to_Ne20
+        -1, 0, 10, -1, -1, 13, 15,  // p_O17_to_F18
+        -1, 0, 11, -1, -1, 14, 17,  // p_O18_to_F19
+        -1, 0, 12, -1, -1, 15, 19,  // p_F17_to_Ne18
+        -1, 0, 13, -1, -1, 16, 21,  // p_F18_to_Ne19
+        -1, 0, 14, -1, -1, 17, 23,  // p_F19_to_Ne20
+        -1, 1, 15, -1, -1, 18, 25,  // He4_Ne18_to_Mg22
+        -1, 1, 17, -1, -1, 19, 26,  // He4_Ne20_to_Mg24
+        -1, 1, 2, -1, 0, 6, -1,  // He4_C12_to_p_N15
+        -1, 2, 2, -1, 1, 17, 65,  // C12_C12_to_He4_Ne20
+        -1, 1, 4, -1, 0, 9, 55,  // He4_N13_to_p_O16
+        -1, 1, 5, -1, 0, 10, -1,  // He4_N14_to_p_O17
+        -1, 0, 6, -1, 1, 2, 47,  // p_N15_to_He4_C12
+        -1, 1, 6, -1, 0, 11, -1,  // He4_N15_to_p_O18
+        -1, 1, 7, -1, 0, 12, 60,  // He4_O14_to_p_F17
+        -1, 1, 8, -1, 0, 13, -1,  // He4_O15_to_p_F18
+        -1, 0, 9, -1, 1, 4, -1,  // p_O16_to_He4_N13
+        -1, 1, 9, -1, 0, 14, -1,  // He4_O16_to_p_F19
+        -1, 2, 9, -1, 1, 19, 66,  // C12_O16_to_He4_Mg24
+        -1, 0, 10, -1, 1, 5, 50,  // p_O17_to_He4_N14
+        -1, 0, 11, -1, 1, 6, 52,  // p_O18_to_He4_N15
+        -1, 0, 12, -1, 1, 7, -1,  // p_F17_to_He4_O14
+        -1, 1, 12, -1, 0, 17, 64,  // He4_F17_to_p_Ne20
+        -1, 0, 13, -1, 1, 8, 54,  // p_F18_to_He4_O15
+        -1, 0, 14, -1, 1, 9, 56,  // p_F19_to_He4_O16
+        -1, 0, 17, -1, 1, 12, -1,  // p_Ne20_to_He4_F17
+        -1, 1, 17, -1, 2, 2, -1,  // He4_Ne20_to_C12_C12
+        -1, 1, 19, -1, 2, 9, -1,  // He4_Mg24_to_C12_O16
+        1, 1, 1, -1, -1, 2, 27  // He4_He4_He4_to_C12
     };
 }
 #endif

networks/CNO_extras/actual_rhs.H

@@ -38,9 +38,10 @@ void ener_gener_rate(T const& dydt, amrex::Real& enuc)
 
     enuc = 0.0_rt;
 
-    for (int n = 1; n <= NumSpec; ++n) {
-        enuc += dydt(n) * network::mion(n);
-    }
+    amrex::constexpr_for<1, NumSpec+1>([&] (auto n)
+    {
+        enuc += dydt(n) * network::mion<n>();
+    });
 
     enuc *= C::enuc_conv2;
 }

networks/CNO_extras/interp_tools.H

@@ -0,0 +1,66 @@
+#ifndef INTERP_TOOLS_H
+#define INTERP_TOOLS_H
+
+#include <AMReX_REAL.H>
+
+namespace interp_net {
+
+    // index locator
+    // return the index i such that x_array[i] <= x0 <= x_array[i+1]
+
+    template <typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    int
+    find_index(const amrex::Real x0, const T& x_array) {
+
+        int left = x_array.lo();
+        int right = x_array.hi();
+
+        int idx = -1;
+
+        if (x0 >= x_array(left) && x0 < x_array(right)) {
+
+            // find the largest x element <= x0 using a
+            // binary search
+
+            while (left < right) {
+                int mid = (left + right) / 2;
+                if (x_array(mid) > x0) {
+                    right = mid;
+                } else {
+                    left = mid + 1;
+                }
+            }
+
+            idx = right - 1;
+        }
+
+        return idx;
+    }
+
+    // an index locator that extrapolates at the boundaries instead
+    // of returning -1
+
+    template <typename T>
+    AMREX_GPU_HOST_DEVICE AMREX_INLINE
+    int
+    find_index_extrap(const amrex::Real x0, const T& x_array) {
+
+        int left = x_array.lo();
+        int right = x_array.hi();
+
+        int idx;
+
+        if (x0 < x_array(left)) {
+            idx = left;
+        } else if (x0 > x_array(right)) {
+            idx = right - 1;
+        } else {
+            idx = find_index(x0, x_array);
+        }
+
+        return idx;
+    }
+
+}
+#endif

networks/CNO_extras/partition_functions.H

@@ -7,46 +7,14 @@
 #include <tfactors.H>
 #include <fundamental_constants.H>
 #include <network_properties.H>
+#include <interp_tools.H>
 
 using namespace amrex::literals;
 using namespace Species;
 
 namespace part_fun {
 
 
-
-    // index locator
-
-    template <typename T>
-    AMREX_GPU_HOST_DEVICE AMREX_INLINE
-    int
-    index_pf(const amrex::Real t9, const T& temp_array) {
-
-        int left = temp_array.lo();
-        int right = temp_array.hi();
-
-        int idx = -1;
-
-        if (t9 >= temp_array(left) && t9 < temp_array(right)) {
-
-            // find the largest temperature element <= t9 using a
-            // binary search
-
-            while (left < right) {
-                int mid = (left + right) / 2;
-                if (temp_array(mid) > t9) {
-                    right = mid;
-                } else {
-                    left = mid + 1;
-                }
-            }
-
-            idx = right - 1;
-        }
-
-        return idx;
-    }
-
     // interpolation routine
 
     template <typename T>
@@ -58,7 +26,7 @@ namespace part_fun {
         // if the index is -1, we are either outside of the bounds or
         // we never computed it.  try recomputing here
         if (idx == -1) {
-            idx = index_pf(t9, temp_array);
+            idx = interp_net::find_index(t9, temp_array);
         }
 
         // if the index is still -1, then we are out-of-bounds, so

networks/CNO_extras/pynucastro-info.txt

@@ -1 +1 @@
-pynucastro version: 2.8.0-6-g7d01fa58a
+pynucastro version: 2.8.0-39-gcd1ab9edf