Consistently use ycur for user RHS routines

CHANGELOG.md

@@ -6,10 +6,20 @@
 
 ### New Features and Enhancements
 
-  The default numbers of stages for the SSP Runge--Kutta methods `ARKODE_LSRK_SSP_S_2`
-  and `ARKODE_LSRK_SSP_S_3` in LSRKStep were changed from 10 and 9, respectively, to
-  their minimum allowable values of 2 and 4.  Users may revert to the previous values
-  by calling `LSRKStepSetNumSSPStages`.
+Multiple minor updates were made to the ARKODE package.  We removed an extraneous
+copy of the output vector when using ARKODE in `ARK_ONE_STEP` mode.  We
+standardized calls to the user-supplied right-hand-side functions so that these
+are provided the user-supplied solution vector passed to `ARKodeEvolve` whenever
+possible -- the notable exceptions are the Hermite temporal interpolation module,
+the provided preconditioners ARKBANDPRE and ARKBBDPRE, banded or dense linear
+solvers with automatically-approximated Jacobian matrices, iterative linear solvers
+with automatically-approximated Jacobian-times-vector product, temporal root-finding,
+discrete adjoint modules in ARKStep or ERKStep, the SPRKStep stepper, and LSRKStep's
+use of the automated dominant eigenvalue estimation module.  The default numbers of
+stages for the SSP Runge--Kutta methods `ARKODE_LSRK_SSP_S_2` and `ARKODE_LSRK_SSP_S_3`
+in LSRKStep were changed from 10 and 9, respectively, to their minimum allowable
+values of 2 and 4.  Users may revert to the previous values by calling
+`LSRKStepSetNumSSPStages`.
 
 ### Bug Fixes
 
@@ -41,8 +51,6 @@ failing the step.
 The functions `CVodeGetUserDataB` and `IDAGetUserDataB` were added to CVODES and
 IDAS, respectively.
 
-Removed extraneous copy of output vector when using ARKODE in `ARK_ONE_STEP` mode.
-
 ### Bug Fixes
 
 Fixed a bug in the CVODE(S) inequality constraint handling where the predicted

doc/shared/RecentChanges.rst

@@ -5,10 +5,20 @@
 
 **New Features and Enhancements**
 
-  The default numbers of stages for the SSP Runge--Kutta methods :c:enumerator:`ARKODE_LSRK_SSP_S_2`
-  and :c:enumerator:`ARKODE_LSRK_SSP_S_3` in LSRKStep were changed from 10 and 9, respectively, to
-  their minimum allowable values of 2 and 4.  Users may revert to the previous values by calling
-  :c:func:`LSRKStepSetNumSSPStages`.
+Multiple minor updates were made to the ARKODE package.  We removed an extraneous
+copy of the output vector when using ARKODE in ``ARK_ONE_STEP`` mode.  We
+standardized calls to the user-supplied right-hand-side functions so that these
+are provided the user-supplied solution vector passed to :c:func:`ARKodeEvolve` whenever
+possible -- the notable exceptions are the Hermite temporal interpolation module,
+the provided preconditioners ARKBANDPRE and ARKBBDPRE, banded or dense linear
+solvers with automatically-approximated Jacobian matrices, iterative linear solvers
+with automatically-approximated Jacobian-times-vector product, temporal root-finding,
+discrete adjoint modules in ARKStep or ERKStep, the SPRKStep stepper, and LSRKStep's
+use of the automated dominant eigenvalue estimation module.  The default numbers of
+stages for the SSP Runge--Kutta methods :c:enumerator:`ARKODE_LSRK_SSP_S_2` and
+:c:enumerator:`ARKODE_LSRK_SSP_S_3` in LSRKStep were changed from 10 and 9, respectively,
+to their minimum allowable values of 2 and 4.  Users may revert to the previous values by
+calling :c:func:`LSRKStepSetNumSSPStages`.
 
 **Bug Fixes**
 

src/arkode/arkode.c

@@ -2534,7 +2534,8 @@ int arkHin(ARKodeMem ark_mem, sunrealtype tout)
        mass matrix solve when computing the full RHS. Before calling arkHin, h
        is set to |tout - tcur| or 1 and so we do not need to guard against
        h == 0 here before calling the full RHS. */
-    retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->yn,
+    N_VScale(ONE, ark_mem->yn, ark_mem->ycur);
+    retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->ycur,
                                    ark_mem->fn, ARK_FULLRHS_START);
     if (retval) { return ARK_RHSFUNC_FAIL; }
     ark_mem->fn_is_current = SUNTRUE;

src/arkode/arkode_arkstep.c

@@ -1845,7 +1845,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
          of the just completed step (tn, yn) and potentially reuse the
          evaluation (FSAL method) or save the value for later use. */
       mode   = (ark_mem->initsetup) ? ARK_FULLRHS_START : ARK_FULLRHS_END;
-      retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->yn,
+      retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->ycur,
                                      ark_mem->fn, mode);
       if (retval)
       {
@@ -1873,7 +1873,7 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
                                           ark_mem->user_data);
           if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
         }
-        retval = step_mem->fi(ark_mem->tn, ark_mem->yn, step_mem->Fi[0],
+        retval = step_mem->fi(ark_mem->tn, ark_mem->ycur, step_mem->Fi[0],
                               ark_mem->user_data);
         step_mem->nfi++;
 

src/arkode/arkode_erkstep.c

@@ -778,15 +778,15 @@ int erkStep_TakeStep(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   SUNLogExtraDebugVec(ARK_LOGGER, "stage", ark_mem->yn, "z_0(:) =");
 
   /* Call the full RHS if needed. If this is the first step then we may need to
-     evaluate or copy the RHS values from an  earlier evaluation (e.g., to
+     evaluate or copy the RHS values from an earlier evaluation (e.g., to
      compute h0). For subsequent steps treat this RHS evaluation as an
      evaluation at the end of the just completed step to potentially reuse
      (FSAL methods) RHS evaluations from the end of the last step. */
 
   if (!(ark_mem->fn_is_current))
   {
     mode   = (ark_mem->initsetup) ? ARK_FULLRHS_START : ARK_FULLRHS_END;
-    retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->yn,
+    retval = ark_mem->step_fullrhs(ark_mem, ark_mem->tn, ark_mem->ycur,
                                    ark_mem->fn, mode);
     if (retval)
     {

src/arkode/arkode_mristep.c

@@ -1320,9 +1320,13 @@ int mriStep_Init(ARKodeMem ark_mem, sunrealtype tout, int init_type)
     }
     if (ark_mem->hin == ZERO)
     {
+      /*   initialize (tcur,ycur) to (t0,y0) */
+      ark_mem->tcur = ark_mem->tn;
+      N_VScale(ONE, ark_mem->yn, ark_mem->ycur);
+
       /*   tempv1 = fslow(t0, y0) */
-      if (mriStep_SlowRHS(ark_mem, ark_mem->tcur, ark_mem->yn, ark_mem->tempv1,
-                          ARK_FULLRHS_START) != ARK_SUCCESS)
+      if (mriStep_SlowRHS(ark_mem, ark_mem->tcur, ark_mem->ycur,
+                          ark_mem->tempv1, ARK_FULLRHS_START) != ARK_SUCCESS)
       {
         arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__, __FILE__,
                         "error calling slow RHS function(s)");
@@ -1791,7 +1795,6 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
   ARKodeMRIStepMem step_mem;          /* outer stepper memory       */
   int is;                             /* current stage index        */
   int retval;                         /* reusable return flag       */
-  N_Vector tmp;                       /* N_Vector pointer           */
   SUNAdaptController_Type adapt_type; /* timestep adaptivity type   */
   sunrealtype t0, tf;                 /* start/end of each stage    */
   sunbooleantype calc_fslow;
@@ -1843,7 +1846,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
   if (!ark_mem->fixedstep)
   {
     retval = mriStepInnerStepper_Reset(step_mem->stepper, ark_mem->tn,
-                                       ark_mem->yn);
+                                       ark_mem->ycur);
     if (retval != ARK_SUCCESS)
     {
       arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__,
@@ -1878,7 +1881,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
   nested_mri = step_mem->expforcing || step_mem->impforcing;
   if (ark_mem->fn == NULL || nested_mri)
   {
-    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->yn,
+    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->ycur,
                               ARK_FULLRHS_START);
     if (retval)
     {
@@ -1901,7 +1904,7 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
   }
   else if (ark_mem->fn != NULL && !ark_mem->fn_is_current)
   {
-    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->yn, ark_mem->fn,
+    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->ycur, ark_mem->fn,
                              ARK_FULLRHS_START);
     if (retval)
     {
@@ -2119,15 +2122,9 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
   {
     is = step_mem->stages;
 
-    /* Temporarily swap ark_mem->ycur and ark_mem->tempv4 pointers, copying
-       data so that both hold the current ark_mem->ycur value.  This ensures
-       that during this embedding "stage":
-         - ark_mem->ycur will be the correct initial condition for the final stage.
-         - ark_mem->tempv4 will hold the embedded solution vector. */
+    /* Copy ark_mem->ycur into ark_mem->tempv4, since it serves as the initial condition
+       for both this embedding stage, and the subsequent final stage. */
     N_VScale(ONE, ark_mem->ycur, ark_mem->tempv4);
-    tmp             = ark_mem->ycur;
-    ark_mem->ycur   = ark_mem->tempv4;
-    ark_mem->tempv4 = tmp;
 
     /* Reset ark_mem->tcur as the time value corresponding with the end of the step */
     /* Set relevant stage times (including desired stage time for implicit solves) */
@@ -2181,10 +2178,11 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
     SUNLogExtraDebugVec(ARK_LOGGER, "embedded solution", ark_mem->ycur,
                         "y_embedded(:) =");
 
-    /* Swap back ark_mem->ycur with ark_mem->tempv4, and reset the inner integrator */
-    tmp             = ark_mem->ycur;
-    ark_mem->ycur   = ark_mem->tempv4;
-    ark_mem->tempv4 = tmp;
+    /* Copy embedding solution into ark_mem->tempv1 for later error estimation,
+       reset ark_mem->ycur to ark_mem->tempv4 in preparation for the final stage,
+       and reset the inner integrator */
+    N_VScale(ONE, ark_mem->ycur, ark_mem->tempv1);
+    N_VScale(ONE, ark_mem->tempv4, ark_mem->ycur);
     retval = mriStepInnerStepper_Reset(step_mem->stepper, t0, ark_mem->ycur);
     if (retval != ARK_SUCCESS)
     {
@@ -2286,10 +2284,11 @@ int mriStep_TakeStepMRIGARK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPt
     }
 
     /* Compute temporal error estimate via difference between step
-       solution and embedding, store in ark_mem->tempv1, and take norm. */
+       solution and embedding (embedding is already stored in ark_mem->tempv1
+       so we need only subtract off ark_mem->ycur) and take norm. */
     if (do_embedding)
     {
-      N_VLinearSum(ONE, ark_mem->tempv4, -ONE, ark_mem->ycur, ark_mem->tempv1);
+      N_VLinearSum(ONE, ark_mem->tempv1, -ONE, ark_mem->ycur, ark_mem->tempv1);
       *dsmPtr = N_VWrmsNorm(ark_mem->tempv1, ark_mem->ewt);
     }
 
@@ -2399,7 +2398,7 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   if (!ark_mem->fixedstep)
   {
     retval = mriStepInnerStepper_Reset(step_mem->stepper, ark_mem->tn,
-                                       ark_mem->yn);
+                                       ark_mem->ycur);
     if (retval != ARK_SUCCESS)
     {
       arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__,
@@ -2422,7 +2421,7 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   nested_mri = step_mem->expforcing || step_mem->impforcing;
   if (ark_mem->fn == NULL || nested_mri)
   {
-    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->yn,
+    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->ycur,
                               ARK_FULLRHS_START);
     if (retval)
     {
@@ -2445,7 +2444,7 @@ int mriStep_TakeStepMRISR(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   }
   if (ark_mem->fn != NULL && !ark_mem->fn_is_current)
   {
-    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->yn, ark_mem->fn,
+    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->ycur, ark_mem->fn,
                              ARK_FULLRHS_START);
     if (retval)
     {
@@ -2881,7 +2880,7 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   /* for adaptive computations, reset the inner integrator to the beginning of this step */
   if (!ark_mem->fixedstep)
   {
-    retval = mriStepInnerStepper_Reset(step_mem->stepper, t0, ark_mem->yn);
+    retval = mriStepInnerStepper_Reset(step_mem->stepper, t0, ark_mem->ycur);
     if (retval != ARK_SUCCESS)
     {
       arkProcessError(ark_mem, ARK_INNERSTEP_FAIL, __LINE__, __func__, __FILE__,
@@ -2903,7 +2902,7 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   nested_mri = step_mem->expforcing || step_mem->impforcing;
   if (ark_mem->fn == NULL || nested_mri)
   {
-    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->yn,
+    retval = mriStep_UpdateF0(ark_mem, step_mem, ark_mem->tn, ark_mem->ycur,
                               ARK_FULLRHS_START);
     if (retval)
     {
@@ -2914,7 +2913,7 @@ int mriStep_TakeStepMERK(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
   }
   else if (ark_mem->fn != NULL && !ark_mem->fn_is_current)
   {
-    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->yn, ark_mem->fn,
+    retval = mriStep_FullRHS(ark_mem, ark_mem->tn, ark_mem->ycur, ark_mem->fn,
                              ARK_FULLRHS_START);
     if (retval)
     {