[WIP] regenerate the networks with the latest pynucastro

networks/3alpha_pyna/partition_functions_data.cpp

@@ -0,0 +1,14 @@
+#include <AMReX_Array.H>
+#include <string>
+#include <table_rates.H>
+#include <AMReX_Print.H>
+
+#include <partition_functions.H>
+
+using namespace amrex;
+
+namespace part_fun {
+
+
+}
+

networks/CNO_He_burn/actual_rhs.H

@@ -1584,7 +1584,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
 
     // Fill approximate rates
 
-    fill_approx_rates<do_T_derivatives, T>(tfactors, rate_eval);
+    fill_approx_rates<do_T_derivatives, T>(tfactors, state.rho, Y, rate_eval);
 
     // Calculate tabular rates
 
@@ -2005,9 +2005,12 @@ void actual_rhs (burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 
     // Get the thermal neutrino losses
 
-    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real sneut{};
+#ifdef NEUTRINOS
     constexpr int do_derivatives{0};
+    amrex::Real dsneutdt{}, dsneutdd{}, dsnuda{}, dsnudz{};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+#endif
 
     // Append the energy equation (this is erg/g/s)
 
@@ -2829,6 +2832,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
 
     // Account for the thermal neutrino losses
 
+#ifdef NEUTRINOS
     amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
     constexpr int do_derivatives{1};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
@@ -2837,7 +2841,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
        amrex::Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);
        jac.add(net_ienuc, j, -b1);
     }
-
+#endif
 
     // Evaluate the Jacobian elements with respect to energy by
     // calling the RHS using d(rate) / dT and then transform them

networks/CNO_He_burn/reaclib_rates.H

@@ -8342,7 +8342,10 @@ 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, [[maybe_unused]] T& rate_eval)
+fill_approx_rates([[maybe_unused]] const tf_t& tfactors,
+                  [[maybe_unused]] const amrex::Real rho,
+                  [[maybe_unused]] const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
+                  [[maybe_unused]] T& rate_eval)
 {
 
     [[maybe_unused]] amrex::Real rate{};

networks/CNO_extras/Make.package

@@ -5,6 +5,7 @@ ifeq ($(USE_REACT),TRUE)
   CEXE_headers += actual_network.H
   CEXE_headers += tfactors.H
   CEXE_headers += partition_functions.H
+  CEXE_sources += partition_functions_data.cpp
   CEXE_headers += actual_rhs.H
   CEXE_headers += reaclib_rates.H
   CEXE_headers += table_rates.H

networks/CNO_extras/actual_network.H

@@ -30,64 +30,64 @@ namespace network
             return 0.0_rt;
         }
         else if constexpr (spec == He4) {
-            return 28.29566_rt;
+            return 28.295662457999697_rt;
         }
         else if constexpr (spec == C12) {
-            return 92.16172800000001_rt;
+            return 92.16173498399803_rt;
         }
         else if constexpr (spec == C13) {
-            return 97.108037_rt;
+            return 97.10804378399916_rt;
         }
         else if constexpr (spec == N13) {
-            return 94.105219_rt;
+            return 94.10522604799917_rt;
         }
         else if constexpr (spec == N14) {
-            return 104.65859599999999_rt;
+            return 104.65860734799753_rt;
         }
         else if constexpr (spec == N15) {
-            return 115.4919_rt;
+            return 115.49190414799887_rt;
         }
         else if constexpr (spec == O14) {
-            return 98.731892_rt;
+            return 98.73189611199996_rt;
         }
         else if constexpr (spec == O15) {
-            return 111.95538_rt;
+            return 111.95539521199862_rt;
         }
         else if constexpr (spec == O16) {
-            return 127.619296_rt;
+            return 127.6193154119992_rt;
         }
         else if constexpr (spec == O17) {
-            return 131.76237600000002_rt;
+            return 131.76239561199873_rt;
         }
         else if constexpr (spec == O18) {
-            return 139.807746_rt;
+            return 139.8077658120019_rt;
         }
         else if constexpr (spec == F17) {
-            return 128.21957600000002_rt;
+            return 128.21958437599824_rt;
         }
         else if constexpr (spec == F18) {
-            return 137.369484_rt;
+            return 137.36950247599816_rt;
         }
         else if constexpr (spec == F19) {
-            return 147.801342_rt;
+            return 147.80136567599766_rt;
         }
         else if constexpr (spec == Ne18) {
-            return 132.142626_rt;
+            return 132.14265544000227_rt;
         }
         else if constexpr (spec == Ne19) {
-            return 143.779517_rt;
+            return 143.7795235400008_rt;
         }
         else if constexpr (spec == Ne20) {
-            return 160.6448_rt;
+            return 160.64482384000075_rt;
         }
         else if constexpr (spec == Mg22) {
-            return 168.58074200000001_rt;
+            return 168.58082376800303_rt;
         }
         else if constexpr (spec == Mg24) {
-            return 198.25701600000002_rt;
+            return 198.2570479679962_rt;
         }
         else if constexpr (spec == Fe56) {
-            return 492.2598239999999_rt;
+            return 492.2599506639962_rt;
         }
 
 

networks/CNO_extras/actual_rhs.H

@@ -687,7 +687,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
 
     // Fill approximate rates
 
-    fill_approx_rates<do_T_derivatives, T>(tfactors, rate_eval);
+    fill_approx_rates<do_T_derivatives, T>(tfactors, state.rho, Y, rate_eval);
 
     // Calculate tabular rates
 
@@ -973,9 +973,12 @@ void actual_rhs (burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 
     // Get the thermal neutrino losses
 
-    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real sneut{};
+#ifdef NEUTRINOS
     constexpr int do_derivatives{0};
+    amrex::Real dsneutdt{}, dsneutdd{}, dsnuda{}, dsnudz{};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+#endif
 
     // Append the energy equation (this is erg/g/s)
 
@@ -1497,6 +1500,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
 
     // Account for the thermal neutrino losses
 
+#ifdef NEUTRINOS
     amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
     constexpr int do_derivatives{1};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
@@ -1505,7 +1509,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
        amrex::Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);
        jac.add(net_ienuc, j, -b1);
     }
-
+#endif
 
     // Evaluate the Jacobian elements with respect to energy by
     // calling the RHS using d(rate) / dT and then transform them

networks/CNO_extras/partition_functions.H

@@ -17,21 +17,21 @@ namespace part_fun {
 
     // interpolation routine
 
-    template <int npts>
+    template <typename T>
     AMREX_GPU_HOST_DEVICE AMREX_INLINE
-    void interpolate_pf(const amrex::Real t9, const amrex::Real (&temp_array)[npts], const amrex::Real (&pf_array)[npts],
+    void interpolate_pf(const amrex::Real t9, const T& temp_array, const T& pf_array,
                         amrex::Real& pf, amrex::Real& dpf_dT) {
 
-        if (t9 >= temp_array[0] && t9 < temp_array[npts-1]) {
+        if (t9 >= temp_array.lo() && t9 < temp_array.hi()) {
 
             // find the largest temperature element <= t9 using a binary search
 
-            int left = 0;
-            int right = npts;
+            int left = temp_array.lo();
+            int right = temp_array.hi();
 
             while (left < right) {
                 int mid = (left + right) / 2;
-                if (temp_array[mid] > t9) {
+                if (temp_array(mid) > t9) {
                     right = mid;
                 } else {
                     left = mid + 1;
@@ -44,11 +44,12 @@ namespace part_fun {
 
             // construct the slope -- this is (log10(pf_{i+1}) - log10(pf_i)) / (T_{i+1} - T_i)
 
-            amrex::Real slope = (pf_array[idx+1] - pf_array[idx]) / (temp_array[idx+1] - temp_array[idx]);
+            amrex::Real slope = (pf_array(idx+1) - pf_array(idx)) /
+                                (temp_array(idx+1) - temp_array(idx));
 
             // find the PF
 
-            amrex::Real log10_pf = pf_array[idx] + slope * (t9 - temp_array[idx]);
+            amrex::Real log10_pf = pf_array(idx) + slope * (t9 - temp_array(idx));
             pf = std::pow(10.0_rt, log10_pf);
 
             // find the derivative (with respect to T, not T9)

networks/CNO_extras/partition_functions_data.cpp

@@ -0,0 +1,14 @@
+#include <AMReX_Array.H>
+#include <string>
+#include <table_rates.H>
+#include <AMReX_Print.H>
+
+#include <partition_functions.H>
+
+using namespace amrex;
+
+namespace part_fun {
+
+
+}
+

networks/CNO_extras/reaclib_rates.H

@@ -4483,7 +4483,10 @@ 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, [[maybe_unused]] T& rate_eval)
+fill_approx_rates([[maybe_unused]] const tf_t& tfactors,
+                  [[maybe_unused]] const amrex::Real rho,
+                  [[maybe_unused]] const amrex::Array1D<amrex::Real, 1, NumSpec>& Y,
+                  [[maybe_unused]] T& rate_eval)
 {
 
     [[maybe_unused]] amrex::Real rate{};

networks/ECSN/Make.package

@@ -5,6 +5,7 @@ ifeq ($(USE_REACT),TRUE)
   CEXE_headers += actual_network.H
   CEXE_headers += tfactors.H
   CEXE_headers += partition_functions.H
+  CEXE_sources += partition_functions_data.cpp
   CEXE_headers += actual_rhs.H
   CEXE_headers += reaclib_rates.H
   CEXE_headers += table_rates.H

networks/ECSN/actual_network.H

@@ -30,34 +30,34 @@ namespace network
             return 0.0_rt;
         }
         else if constexpr (spec == He4) {
-            return 28.29566_rt;
+            return 28.295662457999697_rt;
         }
         else if constexpr (spec == O16) {
-            return 127.619296_rt;
+            return 127.6193154119992_rt;
         }
         else if constexpr (spec == O20) {
-            return 151.3714_rt;
+            return 151.37138571200194_rt;
         }
         else if constexpr (spec == F20) {
-            return 154.40268_rt;
+            return 154.40270167600102_rt;
         }
         else if constexpr (spec == Ne20) {
-            return 160.6448_rt;
+            return 160.64482384000075_rt;
         }
         else if constexpr (spec == Mg24) {
-            return 198.25701600000002_rt;
+            return 198.2570479679962_rt;
         }
         else if constexpr (spec == Al27) {
-            return 224.951931_rt;
+            return 224.95193723199915_rt;
         }
         else if constexpr (spec == Si28) {
-            return 236.536832_rt;
+            return 236.53684539599638_rt;
         }
         else if constexpr (spec == P31) {
-            return 262.91617699999995_rt;
+            return 262.9161999600037_rt;
         }
         else if constexpr (spec == S32) {
-            return 271.78012800000005_rt;
+            return 271.78016372399725_rt;
         }
 
 

networks/ECSN/actual_rhs.H

@@ -269,7 +269,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
 
     // Fill approximate rates
 
-    fill_approx_rates<do_T_derivatives, T>(tfactors, rate_eval);
+    fill_approx_rates<do_T_derivatives, T>(tfactors, state.rho, Y, rate_eval);
 
     // Calculate tabular rates
 
@@ -499,9 +499,12 @@ void actual_rhs (burn_t& state, amrex::Array1D<amrex::Real, 1, neqs>& ydot)
 
     // Get the thermal neutrino losses
 
-    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real sneut{};
+#ifdef NEUTRINOS
     constexpr int do_derivatives{0};
+    amrex::Real dsneutdt{}, dsneutdd{}, dsnuda{}, dsnudz{};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
+#endif
 
     // Append the energy equation (this is erg/g/s)
 
@@ -717,6 +720,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
 
     // Account for the thermal neutrino losses
 
+#ifdef NEUTRINOS
     amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
     constexpr int do_derivatives{1};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
@@ -725,7 +729,7 @@ void actual_jac(const burn_t& state, MatrixType& jac)
        amrex::Real b1 = (-state.abar * state.abar * dsnuda + (zion[j-1] - state.zbar) * state.abar * dsnudz);
        jac.add(net_ienuc, j, -b1);
     }
-
+#endif
 
     // Evaluate the Jacobian elements with respect to energy by
     // calling the RHS using d(rate) / dT and then transform them