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