reorder some loops over Jacobians (#1639)

the MathArray2D is Fortran-ordered.  Some of the loops should be more
efficient if their order is swapped.

Author: zingale

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);
         }
     }