From 8924bd8d2c17abf930fa1e0301cb05266957d424 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 31 Jul 2024 14:06:17 -0400 Subject: [PATCH 1/4] reorder some loops over Jacobians (#1639) 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..dead4a654 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) = 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); } } 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 14b8b0e3173041968943d4bbac2c4803a33abceb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 1 Aug 2024 07:55:46 -0400 Subject: [PATCH 2/4] update CHANGES.md for 24.08 (#1632) --- CHANGES.md | 51 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index 5f91931e1..0b17ad85d 100644 --- a/CHANGES.md +++ b/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 From 115aec8846801361d790190c74e1c8c2b5addcaa Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 5 Aug 2024 18:19:27 -0400 Subject: [PATCH 3/4] fix the loop order over Jac in BackwardEuler (#1640) the MathArray is Fortran-ordered --- integration/BackwardEuler/be_integrator.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H index 1a8e5c954..f82cab471 100644 --- a/integration/BackwardEuler/be_integrator.H +++ b/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); From 8ce3375afa5ecd900a359c584fbcf78fddb0234e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 5 Aug 2024 18:19:52 -0400 Subject: [PATCH 4/4] only zero the Jacobian once (#1636) now it is the role of the jac() wrapper to zero the analytic Jacobian. The networks can assume that it is zero. Previously both VODE and the network were doing the zeroing, which was redundant. --- integration/VODE/vode_dvjac.H | 3 --- sphinx_docs/source/integrators.rst | 5 +++++ 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/integration/VODE/vode_dvjac.H b/integration/VODE/vode_dvjac.H index 848a31621..8c420e6f1 100644 --- a/integration/VODE/vode_dvjac.H +++ b/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 diff --git a/sphinx_docs/source/integrators.rst b/sphinx_docs/source/integrators.rst index 1e7b5fe54..cf6c87565 100644 --- a/sphinx_docs/source/integrators.rst +++ b/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 ^^^^^^^^^^^^^^^^