From 1af7a784d8579a97bf5d60fe3ec25c84fe1f3b6d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 31 Jul 2024 11:43:43 -0400 Subject: [PATCH 1/3] reorder some loops over Jacobians the MathArray2D is Fortran-ordered. Some of the loops should be more efficient if their order is swapped. --- integration/integrator_rhs_sdc.H | 28 ++++++++++---------------- integration/integrator_rhs_strang.H | 9 +-------- integration/utils/numerical_jacobian.H | 13 +++--------- 3 files changed, 15 insertions(+), 35 deletions(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index a953d253b..ab3fcddfd 100644 --- a/integration/integrator_rhs_sdc.H +++ b/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) *= [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) *= 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); } } diff --git a/integration/integrator_rhs_strang.H b/integration/integrator_rhs_strang.H index 8c9d6f4a9..90b35a167 100644 --- a/integration/integrator_rhs_strang.H +++ b/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. diff --git a/integration/utils/numerical_jacobian.H b/integration/utils/numerical_jacobian.H index 7c367fcee..982238041 100644 --- a/integration/utils/numerical_jacobian.H +++ b/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); } } From b72771fdd79e813e638af8a3036b3921467e202f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 31 Jul 2024 11:46:32 -0400 Subject: [PATCH 2/3] fix compilation --- integration/integrator_rhs_sdc.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index ab3fcddfd..800d1c20f 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -105,7 +105,7 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) { - pd.zero() + pd.zero(); return; } From 6b8b917a1c50cd3702fcf360e361eadfc0ba4c56 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 31 Jul 2024 11:49:49 -0400 Subject: [PATCH 3/3] fix compilation --- integration/integrator_rhs_sdc.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/integration/integrator_rhs_sdc.H b/integration/integrator_rhs_sdc.H index 800d1c20f..dead4a654 100644 --- a/integration/integrator_rhs_sdc.H +++ b/integration/integrator_rhs_sdc.H @@ -118,13 +118,13 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd) for (int jcol = 1; jcol <= neqs; ++jcol) { for (int irow = 1; irow <= NumSpec; irow++) { - pd(irow, jcol) *= [irow-1]; + pd(irow, jcol) = pd(irow, jcol) * aion[irow-1]; } } for (int jcol = 1; jcol <= NumSpec; ++jcol) { for (int irow = 1; irow <= neqs; ++irow) { - pd(irow, jcol) *= aion_inv[jcol-1]; + pd(irow, jcol) = pd(irow, jcol) * aion_inv[jcol-1]; } }