Skip to content

Feature/stageprocessing#813

Open
drreynolds wants to merge 24 commits intofeature/fes-stagingfrom
feature/stageprocessing
Open

Feature/stageprocessing#813
drreynolds wants to merge 24 commits intofeature/fes-stagingfrom
feature/stageprocessing

Conversation

@drreynolds
Copy link
Collaborator

This PR adds step preprocessing, failed-step postprocessing, and stage RHS preprocessing to ARKODE.

While existing support for stage and step postprocessing has been useful for some advanced collaborator codes, the lack of symmetry increases challenges when applications must compute auxiliary data (e.g., solving Poisson equations for a potential field) but want do not wish to do this more than absolutely necessary. Additionally, when a step fails our existing step postprocessing is not called, eliminating the possibility for applications to rewind auxiliary data to the previous "accepted" step values.

As with previous postprocessing support, these should be considered advanced/experimental features. As such, they are not advertised in the SUNDIALS documentation

Step preprocessing vs postprocessing

The fundamental difference between these is when the processing routine is called inside the ARKODE time-stepping loop. It is unlikely that a user would use both, but this allows them to more precisely place their processing function.

Stage RHS preprocessing vs stage postprocessing

Stage RHS preprocessing occurs just prior to any set of RHS calls that would use an identical state (e.g., if both explicit and implicit RHS functions would be called on the same input), allowing auxiliary data (e.g., MPI exchange buffers) to be reused in both calls.

On the other hand, stage postprocessing is called just after updating the current stage solution.

Again, it is unlikely that a user would leverage both in the same application, but this allows increased precision over the timings when they are performed.

@drreynolds
Copy link
Collaborator Author

This is the first in a sequence of PRs for ARKODE to support our FES collaborators. Please review this one first. When the stack is complete, we'll do a final PR into develop

@gardner48 gardner48 added this to the SUNDIALS Next milestone Jan 20, 2026
@drreynolds
Copy link
Collaborator Author

This PR should be ready to review -- all CI issues resulting from the new release have been resolved.


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

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.

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

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

}

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

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?


/* NULL argument sets default, otherwise set inputs */
ark_mem->ProcessStage = ProcessStage;
ark_mem->PostProcessStage = ProcessStage;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
ark_mem->PostProcessStage = ProcessStage;
ark_mem->PostProcessStage = ProcessStage;
ark_mem->ps_data = ark_mem->user_data;

ark_mem = (ARKodeMem)arkode_mem;

/* NULL argument sets default, otherwise set inputs */
ark_mem->PreProcessRHS = PreprocessRHS;
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
ark_mem->PreProcessRHS = PreprocessRHS;
ark_mem->PreProcessRHS = PreprocessRHS;
ark_mem->ps_data = ark_mem->user_data;

if (retval != 0)
{
SUNLogInfo(ARK_LOGGER, "end-stages-list",
"status = failed preprocess stage, retval = %i", retval);
Copy link
Member

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);

if (retval != 0)
{
SUNLogInfo(ARK_LOGGER, "end-stages-list",
"status = failed preprocess stage, retval = %i", retval);
Copy link
Member

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);

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

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.

Copy link
Member

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_Adjoint is missing a PreProcessRHS call, does it need one? Similarly, should there be a PreProcessRHS call before adj_f calls?

Copy link
Member

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?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants