Skip to content

Commit

Permalink
Move pivot into the VODE state (#1456)
Browse files Browse the repository at this point in the history
This addresses #1451, which was caused by the pivot array being reused for the first time in a new step when the Jacobian had been cached without keeping the corresponding values of the pivot array in scope.
  • Loading branch information
maxpkatz authored Jan 27, 2024
1 parent 0c58dd5 commit 0e28d42
Show file tree
Hide file tree
Showing 4 changed files with 10 additions and 11 deletions.
7 changes: 3 additions & 4 deletions integration/VODE/vode_dvjac.H
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,9 @@
#include <integrator_rhs_simplified_sdc.H>
#endif

template <typename IArray1D, typename BurnT, typename DvodeT>
template <typename BurnT, typename DvodeT>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void dvjac (IArray1D& pivot, int& IERPJ, BurnT& state, DvodeT& vstate)
void dvjac (int& IERPJ, BurnT& state, DvodeT& vstate)
{
// dvjac is called by dvnlsd to compute and process the matrix
// P = I - h*rl1*J , where J is an approximation to the Jacobian
Expand Down Expand Up @@ -169,11 +169,10 @@ void dvjac (IArray1D& pivot, int& IERPJ, BurnT& state, DvodeT& vstate)
int IER;

#ifdef NEW_NETWORK_IMPLEMENTATION
amrex::ignore_unused(pivot);
RHS::dgefa(vstate.jac);
IER = 0;
#else
dgefa<int_neqs>(vstate.jac, pivot, IER);
dgefa<int_neqs>(vstate.jac, vstate.pivot, IER);
#endif

if (IER != 0) {
Expand Down
8 changes: 4 additions & 4 deletions integration/VODE/vode_dvnlsd.H
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@
#endif
#include <vode_dvjac.H>

template <typename IArray1D, typename BurnT, typename DvodeT>
template <typename BurnT, typename DvodeT>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Real dvnlsd (IArray1D& pivot, int& NFLAG, BurnT& state, DvodeT& vstate)
Real dvnlsd (int& NFLAG, BurnT& state, DvodeT& vstate)
{
constexpr int int_neqs = integrator_neqs<BurnT>();

Expand Down Expand Up @@ -76,7 +76,7 @@ Real dvnlsd (IArray1D& pivot, int& NFLAG, BurnT& state, DvodeT& vstate)
// to 0 as an indicator that this has been done.

int IERPJ;
dvjac(pivot, IERPJ, state, vstate);
dvjac(IERPJ, state, vstate);
vstate.IPUP = 0;
vstate.RC = 1.0_rt;
vstate.DRC = 0.0_rt;
Expand Down Expand Up @@ -114,7 +114,7 @@ Real dvnlsd (IArray1D& pivot, int& NFLAG, BurnT& state, DvodeT& vstate)
#ifdef NEW_NETWORK_IMPLEMENTATION
RHS::dgesl(vstate.jac, vstate.y);
#else
dgesl<int_neqs>(vstate.jac, pivot, vstate.y);
dgesl<int_neqs>(vstate.jac, vstate.pivot, vstate.y);
#endif

if (vstate.RC != 1.0_rt) {
Expand Down
4 changes: 1 addition & 3 deletions integration/VODE/vode_dvstep.H
Original file line number Diff line number Diff line change
Expand Up @@ -132,8 +132,6 @@ int dvstep (BurnT& state, DvodeT& vstate)

}

Array1D<short, 1, int_neqs> pivot;

// Compute the predicted values by effectively
// multiplying the yh array by the Pascal triangle matrix.
// dvset is called to calculate all integration coefficients.
Expand Down Expand Up @@ -174,7 +172,7 @@ int dvstep (BurnT& state, DvodeT& vstate)

// Call the nonlinear system solver.

Real ACNRM = dvnlsd(pivot, NFLAG, state, vstate);
Real ACNRM = dvnlsd(NFLAG, state, vstate);

if (NFLAG != 0) {

Expand Down
2 changes: 2 additions & 0 deletions integration/VODE/vode_type.H
Original file line number Diff line number Diff line change
Expand Up @@ -159,6 +159,8 @@ struct dvode_t

Array1D<Real, 1, int_neqs> ewt, savf;

Array1D<short, 1, int_neqs> pivot;

// Array of size NEQ used for the accumulated corrections on each
// step, scaled in the output to represent the estimated local
// error in Y on the last step. This is the vector e in the
Expand Down

0 comments on commit 0e28d42

Please sign in to comment.