Merge branch 'development' into README

CHANGES.md

@@ -1,3 +1,54 @@
+# 24.08
+
+  * autodiff is now used with the templated reaction networks (#1614)
+    + some autodiff clean-ups and derivative fixes (#1604, #1612,
+    #1613, #1616, #1619, #1633)
+
+  * we can now output warnings from GPUs if you compile with
+    `USE_GPU_PRINTF=TRUE` (#1629, #1635)
+
+  * documentation improvements (#1570, #1628)
+
+  * a new jacobian unit (`jac_cell`) test was added that compares the
+    numerical and analytic Jacobians (#1618)
+
+  * support for Strang + NSE has been removed.  NSE only works with
+    SDC now (#1549, #1621)
+
+  * the network `CNO_He_burn` was added for explosive H/He burning
+    (#1622)
+
+  * code clean-ups (#1582, #1602, #1609, #1623, #1624, #1625, #1626,
+    #1627, #1631, #1639)
+
+  * `test_nse_net` now also tests the NSE EOS interface (#1621)
+
+  * the self-consistent NSE + SDC update has been synced with the
+    tabular NSE implementation (#1569, #1607, #1617)
+
+  * `test_jac` was not correctly evaluating the numerical Jacobian
+    (#1615)
+
+  * the `fast_atan` function is now more accurate (#1611)
+
+  * `test_ase` was renamed `test_nse_net` and the old `test_nse` was
+    removed (#1610)
+
+  * the old `test_screening` unit test was removed (#1608)
+
+  * the RKC integrator now supports NSE bailout (#1544)
+
+  * a second temperature check for tabular NSE was added -- above this
+    temperature we don't consider composition (#1547)
+
+  * a SDC+NSE unit test was added (#1548)
+
+  * a fast log and fast pow approximation was added (#1591)
+
+  * the primordial_chem network now uses the fast math routines (#1605)
+
+  * fix potential Inf in constexpr linear algebra (#1603)
+
 # 24.07
 
    * added an autodiff library and converted all of the screening

integration/BackwardEuler/be_integrator.H

@@ -76,8 +76,8 @@ int single_step (BurnT& state, BeT& be, const amrex::Real dt)
         // construct the matrix for the linear system
         // (I - dt J) dy^{n+1} = rhs
 
-        for (int m  = 1; m <= int_neqs; m++) {
-            for (int n = 1; n <= int_neqs; n++) {
+        for (int n = 1; n <= int_neqs; n++) {
+            for (int m  = 1; m <= int_neqs; m++) {
                 be.jac(m, n) *= -dt;
                 if (m == n) {
                     be.jac(m, n) = 1.0_rt + be.jac(m, n);

integration/VODE/vode_dvjac.H

@@ -79,9 +79,6 @@ void dvjac (int& IERPJ, BurnT& state, DvodeT& vstate)
             // Indicate that the Jacobian is current for this solve.
             vstate.JCUR = 1;
 
-            // Initialize the Jacobian to zero
-            vstate.jac.zero();
-
             jac(vstate.tn, state, vstate, vstate.jac);
 
 #ifdef ALLOW_JACOBIAN_CACHING

integration/integrator_rhs_sdc.H

@@ -105,15 +105,8 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)
 
 
     if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) {
-
-        for (int j = 1; j <= INT_NEQS; ++j) {
-            for (int i = 1; i <= INT_NEQS; ++i) {
-                pd(i,j) = 0.0_rt;
-            }
-        }
-
+        pd.zero();
         return;
-
     }
 
     // Call the specific network routine to get the Jacobian.
@@ -122,15 +115,16 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)
 
     // The Jacobian from the nets is in terms of dYdot/dY, but we want
     // it was dXdot/dX, so convert here.
-    for (int n = 1; n <= NumSpec; n++) {
-        for (int m = 1; m <= neqs; m++) {
-            pd(n,m) = pd(n,m) * aion[n-1];
+
+    for (int jcol = 1; jcol <= neqs; ++jcol) {
+        for (int irow = 1; irow <= NumSpec; irow++) {
+            pd(irow, jcol) = pd(irow, jcol) * aion[irow-1];
         }
     }
 
-    for (int m = 1; m <= neqs; m++) {
-        for (int n = 1; n <= NumSpec; n++) {
-            pd(m,n) = pd(m,n) * aion_inv[n-1];
+    for (int jcol = 1; jcol <= NumSpec; ++jcol) {
+        for (int irow = 1; irow <= neqs; ++irow) {
+            pd(irow, jcol) = pd(irow, jcol) * aion_inv[jcol-1];
         }
     }
 
@@ -170,9 +164,9 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)
 
     eos_xderivs_t eos_xderivs = composition_derivatives(eos_state);
 
-    for (int m = 1; m <= neqs; m++) {
-        for (int n = 1; n <= NumSpec; n++) {
-            pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc);
+    for (int jcol = 1; jcol <= NumSpec;++jcol) {
+        for (int irow = 1; irow <= neqs; ++irow) {
+            pd(irow, jcol) -= eos_xderivs.dedX[jcol-1] * pd(irow, net_ienuc);
         }
     }
 

integration/integrator_rhs_strang.H

@@ -119,15 +119,8 @@ void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, M
     // bounds. Otherwise set the Jacobian to zero and return.
 
     if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) {
-
-        for (int j = 1; j <= INT_NEQS; ++j) {
-            for (int i = 1; i <= INT_NEQS; ++i) {
-                pd(i,j) = 0.0_rt;
-            }
-        }
-
+        pd.zero();
         return;
-
     }
 
     // Call the specific network routine to get the Jacobian.

integration/utils/numerical_jacobian.H

@@ -136,15 +136,8 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac)
     state_delp.T += dy;
 
     if (state_delp.T <= EOSData::mintemp || state_delp.T >= integrator_rp::MAX_TEMP) {
-
-        for (int i = 1; i <= int_neqs; i++) {
-            for (int j = 1; j <= int_neqs; j++) {
-                jac(i,j) = 0.0_rt;
-            }
-        }
-
+        jac.zero();
         return;
-
     }
 
 
@@ -195,8 +188,8 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac)
     // now correct the species derivatives
     // this constructs dy/dX_k |_e = dy/dX_k |_T - e_{X_k} |_T dy/dT / c_v
 
-    for (int m = 1; m <= int_neqs; m++) {
-        for (int n = 1; n <= NumSpec; n++) {
+    for (int n = 1; n <= NumSpec; n++) {
+        for (int m = 1; m <= int_neqs; m++) {
             jac(m, n) -= eos_xderivs.dedX[n-1] * jac(m, net_ienuc);
         }
     }

sphinx_docs/source/integrators.rst

@@ -276,6 +276,11 @@ The form looks like:
    A network is not required to provide a Jacobian if a numerical
    Jacobian is used.
 
+.. important::
+
+   The integrator does not zero the Jacobian elements.  It is the responsibility
+   of the Jacobian implementation to zero the Jacobian array if necessary.
+
 
 Jacobian wrapper
 ^^^^^^^^^^^^^^^^