Skip to content

Feature/stageprocessing #813

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
drreynolds wants to merge 24 commits into
base: feature/fes-staging
Choose a base branch
from
Open

Feature/stageprocessing #813

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
328eea1
Initial implementation of step/stage preprocessing, as well as failed…
drreynolds Jan 9, 2026
9195cda
Regenerated SWIG
drreynolds Jan 13, 2026
876b016
Ran formatter
drreynolds Jan 13, 2026
0fa2494
updated pattern for stage preprocessing in ERKStep
drreynolds Jan 14, 2026
c8bf4cb
updated pattern for stage preprocessing in ARKStep
drreynolds Jan 14, 2026
4b24a20
updated pattern for stage preprocessing in SPRKStep
drreynolds Jan 14, 2026
e16f6f2
updated pattern for stage preprocessing in SPRKStep
drreynolds Jan 14, 2026
0e0b8d2
updated pattern for stage preprocessing in MRIStep; fixed ordering of…
drreynolds Jan 14, 2026
84f706b
updated pattern for stage preprocessing in LSRKStep
drreynolds Jan 14, 2026
98864f9
Renamed PreProcessStage to PreProcessRHS to clarify its role
drreynolds Jan 14, 2026
f8736ea
Regenerated SWIG
drreynolds Jan 14, 2026
d6d209e
Ran formatter
drreynolds Jan 14, 2026
801108c
Updated logging output to reflect reordering of explicit and implicit…
drreynolds Jan 14, 2026
a76be01
One more .out file
drreynolds Jan 14, 2026
1036279
Merge branch 'develop' into feature/stageprocessing
drreynolds Jan 26, 2026
2f6a251
Updated ARKODE time stepping attempt loop to initialize (tcur,ycur) w…
drreynolds Jan 27, 2026
019702b
Updated logging output files since tcur was previously inaccurate in …
drreynolds Jan 28, 2026
ed6d1f4
Updated answers submodule commit
drreynolds Jan 28, 2026
66985fa
Merge branch 'feature/fes-staging' into feature/stageprocessing
drreynolds Jan 28, 2026
24792e5
Merge branch 'feature/fes-staging' into feature/stageprocessing
drreynolds Feb 3, 2026
4fb8e68
Added preprocess-rhs before calls to the implicit RHS in the nonlinea…
drreynolds Feb 3, 2026
c9b0b06
Added some conditionals to only call PreProcessRHS if we will imminen…
drreynolds Feb 4, 2026
2920ca4
Formatting
drreynolds Feb 4, 2026
e9af994
Fixed parenthesis in conditional
drreynolds Feb 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions doc/arkode/guide/source/Constants.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -470,9 +470,9 @@ contains the ARKODE output constants.
Commented-out table rows:

+-------------------------------------+------+------------------------------------------------------------+
| :index:`ARK_POSTPROCESS_STEP_FAIL` | -37 | An error occurred when calling the user-provided |
| :index:`ARK_POSTPROCESS_STEP_FAIL` | -37 | An error occurred when calling a user-provided. |
| | | step-based :c:func:`ARKPostProcessFn` routine. |
+-------------------------------------+------+------------------------------------------------------------+
| :index:`ARK_POSTPROCESS_STAGE_FAIL` | -38 | An error occurred when calling the user-provided |
| :index:`ARK_POSTPROCESS_STAGE_FAIL` | -38 | An error occurred when calling a user-provided |
| | | stage-based :c:func:`ARKPostProcessFn` routine. |
+-------------------------------------+------+------------------------------------------------------------+
9 changes: 8 additions & 1 deletion include/arkode/arkode.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,8 @@ extern "C" {
#define ARK_INNERTOOUTER_FAIL -36

/* ARK_POSTPROCESS_FAIL equals ARK_POSTPROCESS_STEP_FAIL
for backwards compatibility */
for backwards compatibility. Note that we use these
same constants for step and stage preprocessing errors */
Copy link
Member

@gardner48 gardner48 Feb 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
same constants for step and stage preprocessing errors */
same constants for step and RHS preprocessing errors */

#define ARK_POSTPROCESS_FAIL -37
#define ARK_POSTPROCESS_STEP_FAIL -37
#define ARK_POSTPROCESS_STAGE_FAIL -38
Expand Down Expand Up @@ -276,8 +277,14 @@ SUNDIALS_EXPORT int ARKodeClearStopTime(void* arkode_mem);
SUNDIALS_EXPORT int ARKodeSetFixedStep(void* arkode_mem, sunrealtype hfixed);
SUNDIALS_EXPORT int ARKodeSetStepDirection(void* arkode_mem, sunrealtype stepdir);
SUNDIALS_EXPORT int ARKodeSetUserData(void* arkode_mem, void* user_data);
SUNDIALS_EXPORT int ARKodeSetPreprocessStepFn(void* arkode_mem,
ARKPostProcessFn ProcessStep);
SUNDIALS_EXPORT int ARKodeSetPostprocessStepFn(void* arkode_mem,
ARKPostProcessFn ProcessStep);
SUNDIALS_EXPORT int ARKodeSetPostprocessStepFailFn(void* arkode_mem,
ARKPostProcessFn ProcessStep);
SUNDIALS_EXPORT int ARKodeSetPreprocessRHSFn(void* arkode_mem,
ARKPostProcessFn ProcessRHS);
SUNDIALS_EXPORT int ARKodeSetPostprocessStageFn(void* arkode_mem,
ARKPostProcessFn ProcessStage);

Expand Down
38 changes: 31 additions & 7 deletions src/arkode/arkode.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -904,6 +904,18 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout,
"step = %li, tn = " SUN_FORMAT_G ", h = " SUN_FORMAT_G,
ark_mem->nst + 1, ark_mem->tn, ark_mem->h);

/* fill (tcur,ycur) with the last accepted step solution */
ark_mem->tcur = ark_mem->tn;
Copy link
Member

@gardner48 gardner48 Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note in the change log this fixes an error in the logging output where the time of the first stage would be incorrect when reattempting a step.

N_VScale(ONE, ark_mem->yn, ark_mem->ycur);
Copy link
Member

@gardner48 gardner48 Feb 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In some cases this copy is redundant (i.e., when taking internal steps in normal mode) because ycur was copied into yn in arkCompleteStep. A few ways we could handle this:

  1. Preprocess yn instead of ycur -- no copy needed, FASL methods will need to re-evaluate the RHS function if the vector is modified (this is true for any of the options, see comments below)
  2. Add a flag like ycur_is_copy_of_yn to signal when the copy is needed -- more bookkeeping, avoids any unnecessary copies but still has extra copies compared to passing yn
  3. Move the copy inside if (ark_mem->PreProcessStep != NULL) -- redundant copies only when a step preprocessing function is attached

Some related issues when preprocessing a step:

  1. Should we call the preprocessing function before any initial step computations (e.g., estimating the initial step size)? This might require more bookkeeping to avoid duplicate preprocessing calls on the same state.
  2. This currently does a copy for each step attempt. If yn was preprocessed rather than ycur, we might be able to call PreProcessStep only on the first attempt and skip it on subsequent step attempts.
  3. As noted above, FSAL methods will need to re-evaluate the RHS if the state is modified so the vales are consistent. I think the best solution here might be to add an output flag to the processing functions to indicate if the state has been modified. For now, this would require adding a new function type and we could deprecate the existing processing function type.


/* call the user-supplied step preprocessing function (if it exists) */
if (ark_mem->PreProcessStep != NULL)
{
retval = ark_mem->PreProcessStep(ark_mem->tcur, ark_mem->ycur,
ark_mem->ps_data);
if (retval != 0) { return (ARK_POSTPROCESS_STEP_FAIL); }
}

/* Call time stepper module to attempt a step:
0 => step completed successfully
>0 => step encountered recoverable failure; reduce step if possible
Expand Down Expand Up @@ -1001,6 +1013,14 @@ int ARKodeEvolve(void* arkode_mem, sunrealtype tout, N_Vector yout,
ark_mem->h *= ark_mem->eta;
ark_mem->next_h = ark_mem->hprime = ark_mem->h;

/* since the previous step attempt failed, call the user-supplied step failure postprocessing function (if it exists) */
if (ark_mem->PostProcessStepFail != NULL)
{
retval = ark_mem->PostProcessStepFail(ark_mem->tcur, ark_mem->ycur,
Copy link
Member

@gardner48 gardner48 Feb 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On a failed step ycur could be an invalid state or an intermediate stage which it probably not what the user expects. Maybe another way to handle any logic needed for re-attempting a step would be to add the step attempt counter as an input to PreProcessStep.

This would also address the idea noted above about possibly not needing to call PreProcessStep each attempt if yn is passed as input by having user check the attempt number in their function.

ark_mem->ps_data);
if (retval != 0) { return (ARK_POSTPROCESS_STEP_FAIL); }
}

} /* end looping for step attempts */

/* If step attempt loop succeeded, complete step (update current time, solution,
Expand Down Expand Up @@ -1587,12 +1607,15 @@ ARKodeMem arkCreate(SUNContext sunctx)
ark_mem->VRabstolMallocDone = SUNFALSE;
ark_mem->MallocDone = SUNFALSE;

/* No user-supplied step postprocessing function yet */
ark_mem->ProcessStep = NULL;
ark_mem->ps_data = NULL;
/* No user-supplied step pre- or post-processing functions yet */
ark_mem->PreProcessStep = NULL;
ark_mem->PostProcessStep = NULL;
ark_mem->PostProcessStepFail = NULL;
ark_mem->ps_data = NULL;

/* No user-supplied stage postprocessing function yet */
ark_mem->ProcessStage = NULL;
/* No user-supplied stage pre- or post-processing functions yet */
ark_mem->PreProcessRHS = NULL;
ark_mem->PostProcessStage = NULL;

/* No user_data pointer yet */
ark_mem->user_data = NULL;
Expand Down Expand Up @@ -2727,9 +2750,10 @@ int arkCompleteStep(ARKodeMem ark_mem, sunrealtype dsm)
}

/* apply user-supplied step postprocessing function (if supplied) */
if (ark_mem->ProcessStep != NULL)
if (ark_mem->PostProcessStep != NULL)
{
retval = ark_mem->ProcessStep(ark_mem->tcur, ark_mem->ycur, ark_mem->ps_data);
retval = ark_mem->PostProcessStep(ark_mem->tcur, ark_mem->ycur,
ark_mem->ps_data);
if (retval != 0) { return (ARK_POSTPROCESS_STEP_FAIL); }
}
Copy link
Member

@gardner48 gardner48 Feb 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thinking of the state/RHS consistency for FSAL methods above, this probably needs to move inside the steppers so they can call PostProcessStage or PostProcessesStep as appropriate and remove a potentially redundant processing call


Expand Down
112 changes: 79 additions & 33 deletions src/arkode/arkode_arkstep.c
View file
Open in desktop
Copy link
Member

@gardner48 gardner48 Feb 7, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The RHS call in TakeStep_ERK_Adjoint is missing a PreProcessRHS call, does it need one? Similarly, is a PreProcessRHS call needed before a adj_fe call?

Original file line number Diff line number Diff line change
Expand Up @@ -1330,22 +1330,29 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
/* compute the full RHS */
if (!(ark_mem->fn_is_current))
{
/* compute the explicit component */
if (step_mem->explicit)
/* apply user-supplied stage preprocessing function (if supplied) */
if (ark_mem->PreProcessRHS != NULL)
{
retval = step_mem->fe(t, y, step_mem->Fe[0], ark_mem->user_data);
step_mem->nfe++;
retval = ark_mem->PreProcessRHS(t, y, ark_mem->user_data);
if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
}

/* compute the implicit component */
if (step_mem->implicit)
Copy link
Member

@gardner48 gardner48 Feb 3, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a particular reason to call the implicit RHS first or is this just to make the ordering consistent with what's done in TakeStep?

{
retval = step_mem->fi(t, y, step_mem->Fi[0], ark_mem->user_data);
step_mem->nfi++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__,
__FILE__, MSG_ARK_RHSFUNC_FAILED, t);
return (ARK_RHSFUNC_FAIL);
}

/* compute and store M(t)^{-1} fe */
/* compute and store M(t)^{-1} fi */
if (step_mem->mass_type == MASS_TIMEDEP)
{
retval = step_mem->msolve((void*)ark_mem, step_mem->Fe[0],
retval = step_mem->msolve((void*)ark_mem, step_mem->Fi[0],
step_mem->nlscoef / ark_mem->h);
if (retval)
{
Expand All @@ -1356,22 +1363,22 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
}
}

/* compute the implicit component */
if (step_mem->implicit)
/* compute the explicit component */
if (step_mem->explicit)
{
retval = step_mem->fi(t, y, step_mem->Fi[0], ark_mem->user_data);
step_mem->nfi++;
retval = step_mem->fe(t, y, step_mem->Fe[0], ark_mem->user_data);
step_mem->nfe++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__,
__FILE__, MSG_ARK_RHSFUNC_FAILED, t);
return (ARK_RHSFUNC_FAIL);
}

/* compute and store M(t)^{-1} fi */
/* compute and store M(t)^{-1} fe */
if (step_mem->mass_type == MASS_TIMEDEP)
{
retval = step_mem->msolve((void*)ark_mem, step_mem->Fi[0],
retval = step_mem->msolve((void*)ark_mem, step_mem->Fe[0],
step_mem->nlscoef / ark_mem->h);
if (retval)
{
Expand Down Expand Up @@ -1455,11 +1462,18 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
/* recompute RHS functions */
if (recomputeRHS)
{
/* compute the explicit component */
if (step_mem->explicit)
/* apply user-supplied stage preprocessing function (if supplied) */
if (ark_mem->PreProcessRHS != NULL)
{
retval = step_mem->fe(t, y, step_mem->Fe[0], ark_mem->user_data);
step_mem->nfe++;
retval = ark_mem->PreProcessRHS(t, y, ark_mem->user_data);
if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
}

/* compute the implicit component */
if (step_mem->implicit)
{
retval = step_mem->fi(t, y, step_mem->Fi[0], ark_mem->user_data);
step_mem->nfi++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__,
Expand All @@ -1470,7 +1484,7 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
/* compute and store M(t)^{-1} fi */
if (step_mem->mass_type == MASS_TIMEDEP)
{
retval = step_mem->msolve((void*)ark_mem, step_mem->Fe[0],
retval = step_mem->msolve((void*)ark_mem, step_mem->Fi[0],
step_mem->nlscoef / ark_mem->h);
if (retval)
{
Expand All @@ -1481,11 +1495,11 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
}
}

/* compute the implicit component */
if (step_mem->implicit)
/* compute the explicit component */
if (step_mem->explicit)
{
retval = step_mem->fi(t, y, step_mem->Fi[0], ark_mem->user_data);
step_mem->nfi++;
retval = step_mem->fe(t, y, step_mem->Fe[0], ark_mem->user_data);
step_mem->nfe++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__,
Expand All @@ -1496,7 +1510,7 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
/* compute and store M(t)^{-1} fi */
if (step_mem->mass_type == MASS_TIMEDEP)
{
retval = step_mem->msolve((void*)ark_mem, step_mem->Fi[0],
retval = step_mem->msolve((void*)ark_mem, step_mem->Fe[0],
step_mem->nlscoef / ark_mem->h);
if (retval)
{
Expand Down Expand Up @@ -1564,11 +1578,18 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,

case ARK_FULLRHS_OTHER:

/* compute the explicit component and store in ark_tempv2 */
if (step_mem->explicit)
/* apply user-supplied stage preprocessing function (if supplied) */
if (ark_mem->PreProcessRHS != NULL)
{
retval = step_mem->fe(t, y, ark_mem->tempv2, ark_mem->user_data);
step_mem->nfe++;
retval = ark_mem->PreProcessRHS(t, y, ark_mem->user_data);
if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
}

/* compute the implicit component and store in sdata */
if (step_mem->implicit)
{
retval = step_mem->fi(t, y, step_mem->sdata, ark_mem->user_data);
step_mem->nfi++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__, __FILE__,
Expand All @@ -1577,11 +1598,11 @@ int arkStep_FullRHS(ARKodeMem ark_mem, sunrealtype t, N_Vector y, N_Vector f,
}
}

/* compute the implicit component and store in sdata */
if (step_mem->implicit)
/* compute the explicit component and store in ark_tempv2 */
if (step_mem->explicit)
{
retval = step_mem->fi(t, y, step_mem->sdata, ark_mem->user_data);
step_mem->nfi++;
retval = step_mem->fe(t, y, ark_mem->tempv2, ark_mem->user_data);
step_mem->nfe++;
if (retval != 0)
{
arkProcessError(ark_mem, ARK_RHSFUNC_FAIL, __LINE__, __func__, __FILE__,
Expand Down Expand Up @@ -1845,6 +1866,13 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
}
else
{
/* apply user-supplied stage preprocessing function (if supplied) */
if (ark_mem->PreProcessRHS != NULL)
{
retval = ark_mem->PreProcessRHS(ark_mem->tn, ark_mem->yn,
ark_mem->user_data);
if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
Copy link
Member

@gardner48 gardner48 Feb 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
if (retval != 0) { return (ARK_POSTPROCESS_STAGE_FAIL); }
if (retval != 0)
{
SUNLogInfo(ARK_LOGGER, "end-stages-list",
"status = failed preprocess rhs, retval = %i", retval);
return (ARK_POSTPROCESS_STAGE_FAIL);
}

}
retval = step_mem->fi(ark_mem->tn, ark_mem->yn, step_mem->Fi[0],
ark_mem->user_data);
step_mem->nfi++;
Expand Down Expand Up @@ -2015,10 +2043,10 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
/* apply user-supplied stage postprocessing function (if supplied) */
/* NOTE: with internally inconsistent IMEX methods (c_i^E != c_i^I) the value
of tcur corresponds to the stage time from the implicit table (c_i^I). */
if (ark_mem->ProcessStage != NULL)
if (ark_mem->PostProcessStage != NULL)
{
retval = ark_mem->ProcessStage(ark_mem->tcur, ark_mem->ycur,
ark_mem->user_data);
retval = ark_mem->PostProcessStage(ark_mem->tcur, ark_mem->ycur,
ark_mem->user_data);
if (retval != 0)
{
SUNLogInfo(ARK_LOGGER, "end-stages-list",
Expand Down Expand Up @@ -2068,6 +2096,24 @@ int arkStep_TakeStep_Z(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr)
}
}

/* apply user-supplied stage preprocessing function (if supplied) */
/* NOTE: with internally inconsistent IMEX methods (c_i^E != c_i^I) the value
of tcur corresponds to the stage time from the implicit table (c_i^I). */
if (ark_mem->PreProcessRHS != NULL)
{
if ((step_mem->implicit && !deduce_stage) || (step_mem->explicit))
{
retval = ark_mem->PreProcessRHS(ark_mem->tcur, ark_mem->ycur,
ark_mem->user_data);
if (retval != 0)
{
SUNLogInfo(ARK_LOGGER, "end-stages-list",
"status = failed preprocess stage, retval = %i", retval);
Copy link
Member

@gardner48 gardner48 Feb 6, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
"status = failed preprocess stage, retval = %i", retval);
"status = failed preprocess rhs, retval = %i", retval);

return (ARK_POSTPROCESS_STAGE_FAIL);
}
}
}

/* store implicit RHS (value in Fi[is] is from preceding nonlinear iteration) */
if (step_mem->implicit)
{
Expand Down
Loading
Loading