Merge branch 'development' into joss_paper

Author: zingale

CHANGES.md

@@ -1,6 +1,25 @@
 # Changelog
 
-## 23.03
+## 25.05
+
+  * some clang-tidy cleaning (#1784)
+
+  * with HIP we were disabling inlining due to ROCm issues.  This is
+    now an option, with the default enabling inlining, since it works
+    as expected for ROCm >= 6.3.1 (#1780)
+
+  * clean up the pynucastro scripts that derived reverse rates (#1781)
+
+## 25.04
+
+  * the "he-burn" networks now will look for weak rates for
+    all nuclei, not just the Fe-group (#1763)
+
+  * clean up the ECSN network generation script (#1776)
+
+  * Allow for single-step backward Euler (#1773)
+
+## 25.03
 
   * the `nova2` net was renamed `nova`, and the old `nova`
     net was removed (#1746, #1768)

Docs/source/rp_intro.rst

@@ -6,19 +6,31 @@ The behavior of the network and EOS are controlled by many runtime
 parameters. These parameters are defined in plain-text files
 ``_parameters`` located in the different directories that hold the
 microphysics code. At compile time, a script in the AMReX build
-system, ``findparams.py``, locates all of the ``_parameters`` files that
-are needed for the given choice of network, integrator, and EOS, and
-assembles all of the runtime parameters into a set of header files
-(using the ``write_probin.py`` script).
+system, ``findparams.py``, locates all of the ``_parameters`` files
+that are needed for the given choice of network, integrator, and EOS,
+and other modules, and assembles all of the runtime parameters into a
+set of header files (using the ``write_probin.py`` script).
 
+A ``_parameter`` file starts with a namespace declaration like:
 
-Parameter definitions take the form::
+.. code::
+
+   @namespace: integrator
+
+Followed by a list of parameters.  Each parameter defined after the namespace
+declaration would be prefixed by the namespace in an inputs file or commandline,
+like ``integrator.parameter``.
+
+The individual parameter definitions take the form::
 
     # comment describing the parameter
     name              data-type       default-value      priority
 
-Here, the priority is simply an integer. When two directories
-define the same parameter, but with different defaults, the version of
+where "data-type" can be: ``real``, ``bool``, ``int``, ``string``.
+
+The "priority" is simply an integer that is used to determine
+what happens when two different ``_paramerter`` files define the same
+parameter but with different defaults.  In this case, the version of
 the parameter with the highest priority takes precedence. This allows
 specific implementations to override the general parameter defaults.
 

Make.Microphysics_extern

@@ -1,9 +1,12 @@
 # This is the main include makefile for applications that want to use Microphysics
 # You should set NETWORK_OUTPUT_PATH before including this file
 
-# For HIP, inlining produces large kernels that seem to cause memory issues
+# For some (older?) versions of HIP, inlining produces large kernels
+# that seem to cause memory issues
 ifeq ($(USE_HIP),TRUE)
+  ifeq ($(DISABLE_HIP_INLINE),TRUE)
     EXTRACXXFLAGS += -mllvm -amdgpu-function-calls=true
+  endif
 endif
 
 # for common code paths between simplified-SDC and true-SDC, we'll

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;
         }
 
 
@@ -129,11 +129,11 @@ namespace Rates
         k_p_Al27_to_He4_Mg24 = 12,
         k_He4_Si28_to_p_P31 = 13,
         k_p_P31_to_He4_Si28 = 14,
-        k_F20_to_O20 = 15,
-        k_Ne20_to_F20 = 16,
-        k_O20_to_F20 = 17,
-        k_F20_to_Ne20 = 18,
-        NumRates = k_F20_to_Ne20
+        k_F20_to_Ne20 = 15,
+        k_F20_to_O20 = 16,
+        k_Ne20_to_F20 = 17,
+        k_O20_to_F20 = 18,
+        NumRates = k_O20_to_F20
     };
 
     // number of reaclib rates
@@ -164,10 +164,10 @@ namespace Rates
         "p_Al27_to_He4_Mg24",  // 12,
         "He4_Si28_to_p_P31",  // 13,
         "p_P31_to_He4_Si28",  // 14,
-        "F20_to_O20",  // 15,
-        "Ne20_to_F20",  // 16,
-        "O20_to_F20",  // 17,
-        "F20_to_Ne20"  // 18,
+        "F20_to_Ne20",  // 15,
+        "F20_to_O20",  // 16,
+        "Ne20_to_F20",  // 17,
+        "O20_to_F20"  // 18,
     };
 
 }

networks/ECSN/actual_network_data.cpp

@@ -19,10 +19,10 @@ namespace NSE_INDEX
         -1, 0, 7, -1, 1, 6, 11,
         -1, 1, 8, -1, 0, 9, -1,
         -1, 0, 9, -1, 1, 8, 13,
+        -1, -1, 4, -1, -1, 5, 17,
         -1, -1, 4, -1, -1, 3, -1,
         -1, -1, 5, -1, -1, 4, -1,
-        -1, -1, 3, -1, -1, 4, 15,
-        -1, -1, 4, -1, -1, 5, 16
+        -1, -1, 3, -1, -1, 4, 16
     };
 }
 #endif

networks/ECSN/actual_rhs.H

@@ -10,7 +10,9 @@
 #include <jacobian_utilities.H>
 #include <screen.H>
 #include <microphysics_autodiff.H>
+#ifdef NEUTRINOS
 #include <sneut5.H>
+#endif
 #include <reaclib_rates.H>
 #include <table_rates.H>
 
@@ -269,14 +271,22 @@ 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
 
     [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma;
 
     rate_eval.enuc_weak = 0.0_rt;
 
+    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_v<T, rate_derivs_t>) {
+        rate_eval.dscreened_rates_dT(k_F20_to_Ne20) = drate_dt;
+    }
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);
+
     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;
@@ -301,14 +311,6 @@ void evaluate_rates(const burn_t& state, T& rate_eval) {
     }
     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_v<T, rate_derivs_t>) {
-        rate_eval.dscreened_rates_dT(k_F20_to_Ne20) = drate_dt;
-    }
-    rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);
-
 
 }
 
@@ -340,6 +342,11 @@ void get_ydot_weak(const burn_t& state,
 
     // Calculate tabular rates and get ydot_weak
 
+    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;
+    rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);
+
     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;
@@ -355,11 +362,6 @@ void get_ydot_weak(const burn_t& state,
     rate_eval.screened_rates(k_O20_to_F20) = rate;
     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;
-    rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma);
-
     auto screened_rates = rate_eval.screened_rates;
 
     ydot_nuc(H1) = 0.0_rt;
@@ -372,8 +374,8 @@ void get_ydot_weak(const burn_t& state,
         (-screened_rates(k_O20_to_F20)*Y(O20) + screened_rates(k_F20_to_O20)*Y(F20));
 
     ydot_nuc(F20) =
-        (screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20)) +
-        (-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20));
+        (-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20)) +
+        (screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20));
 
     ydot_nuc(Ne20) =
         (screened_rates(k_F20_to_Ne20)*Y(F20) + -screened_rates(k_Ne20_to_F20)*Y(Ne20));
@@ -427,13 +429,13 @@ void rhs_nuc(const burn_t& state,
         (-screened_rates(k_O20_to_F20)*Y(O20) + screened_rates(k_F20_to_O20)*Y(F20));
 
     ydot_nuc(F20) =
-        (screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20)) +
-        (-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20));
+        (-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20)) +
+        (screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20));
 
     ydot_nuc(Ne20) =
-        (screened_rates(k_F20_to_Ne20)*Y(F20) + -screened_rates(k_Ne20_to_F20)*Y(Ne20)) +
         (screened_rates(k_He4_O16_to_Ne20)*Y(He4)*Y(O16)*state.rho + -screened_rates(k_Ne20_to_He4_O16)*Y(Ne20)) +
-        -screened_rates(k_He4_Ne20_to_Mg24)*Y(He4)*Y(Ne20)*state.rho;
+        -screened_rates(k_He4_Ne20_to_Mg24)*Y(He4)*Y(Ne20)*state.rho +
+        (screened_rates(k_F20_to_Ne20)*Y(F20) + -screened_rates(k_Ne20_to_F20)*Y(Ne20));
 
     ydot_nuc(Mg24) =
         screened_rates(k_He4_Ne20_to_Mg24)*Y(He4)*Y(Ne20)*state.rho +
@@ -499,9 +501,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,15 +722,17 @@ void actual_jac(const burn_t& state, MatrixType& jac)
 
     // Account for the thermal neutrino losses
 
-    amrex::Real sneut, dsneutdt, dsneutdd, dsnuda, dsnudz;
+    amrex::Real dsneutdt{};
+#ifdef NEUTRINOS
+    amrex::Real sneut, dsneutdd, dsnuda, dsnudz;
     constexpr int do_derivatives{1};
     sneut5<do_derivatives>(state.T, state.rho, state.abar, state.zbar, sneut, dsneutdt, dsneutdd, dsnuda, dsnudz);
 
     for (int j = 1; j <= NumSpec; ++j) {
        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/ECSN/ecsn_network_generation.py

@@ -5,10 +5,14 @@
 
 data_list = mylibrary.get_rates()
 
-all_nuclei = ["p","he4","ne20","o20","f20","mg24","al27","o16","si28","s32","p31"]
+all_nuclei = ["p", "he4", "ne20", "o20", "f20",
+              "mg24", "al27", "o16", "si28",
+              "s32", "p31"]
 
-escn_library = mylibrary.linking_nuclei(all_nuclei,with_reverse=True)
-escn_tabular = ["f20--o20-toki","ne20--f20-toki","o20--f20-toki","f20--ne20-toki"]
+escn_library = mylibrary.linking_nuclei(all_nuclei, with_reverse=True)
+
+sl = pyna.SuzukiLibrary()
+tabular_rates = sl.linking_nuclei(["f20", "o20", "ne20"])
 
 rc = pyna.RateCollection(libraries=[escn_library])
 
@@ -36,7 +40,8 @@
     if ydots[rate] >= 1.e-20 and rate.weak == False:
         new_rate_list.append(rate)
 
-wd_net = AmrexAstroCxxNetwork(rates=new_rate_list, rate_files=escn_tabular)
+
+wd_net = AmrexAstroCxxNetwork(rates=new_rate_list, libraries=[tabular_rates])
 wd_net.write_network()
 
 wd_net.plot(rho, T, comp, outfile="ECSN.png",

networks/ECSN/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/ECSN/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 {
+
+
+}
+