diff --git a/CHANGELOG.md b/CHANGELOG.md index f60505d993..7d59c75830 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index e884275a61..8689fdf9f2 100644 --- a/doc/shared/RecentChanges.rst +++ b/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** diff --git a/src/arkode/arkode.c b/src/arkode/arkode.c index 22409fc94a..e63507a8ec 100644 --- a/src/arkode/arkode.c +++ b/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; diff --git a/src/arkode/arkode_arkstep.c b/src/arkode/arkode_arkstep.c index 1b64cefd73..6bb3590d57 100644 --- a/src/arkode/arkode_arkstep.c +++ b/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++; diff --git a/src/arkode/arkode_erkstep.c b/src/arkode/arkode_erkstep.c index 91b50a995e..f30f72a636 100644 --- a/src/arkode/arkode_erkstep.c +++ b/src/arkode/arkode_erkstep.c @@ -778,7 +778,7 @@ 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. */ @@ -786,7 +786,7 @@ int erkStep_TakeStep(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) 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) { diff --git a/src/arkode/arkode_mristep.c b/src/arkode/arkode_mristep.c index 43b799c621..a70a24475d 100644 --- a/src/arkode/arkode_mristep.c +++ b/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) {