Skip to content
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.

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

SUNStepper basics based on MRIStepInnerStepper #463

Open
wants to merge 69 commits into
base: develop
Choose a base branch
from

Conversation

balos1
Copy link
Member

@balos1 balos1 commented Apr 29, 2024

I am opening this PR as a draft to gather feedback. This renames/moves MRIStepInnerStepper to SUNStepper so that @Steven-Roberts and I can extend it for other purposes.

src/arkode/arkode.c Outdated Show resolved Hide resolved
src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper_impl.h Outdated Show resolved Hide resolved
@balos1
Copy link
Member Author

balos1 commented May 6, 2024

@drreynolds @gardner48 I have made the changes we discussed on Tuesday. Let me know what you think.

@balos1 balos1 changed the title move MRIStepInnerStepper to SUNStepper SUNStepper basics based on MRIStepInnerStepper May 6, 2024
Copy link
Collaborator

@drreynolds drreynolds left a comment

Choose a reason for hiding this comment

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

I like how clean this is now. I noticed a few items for comment (see below).

src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
src/arkode/xbraid/arkode_xbraid.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
include/sundials/sundials_stepper.h Show resolved Hide resolved
include/sundials/sundials_stepper.h Outdated Show resolved Hide resolved
include/sundials/sundials_stepper.h Show resolved Hide resolved
include/sundials/sundials_stepper.h Show resolved Hide resolved
src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
include/sundials/sundials_stepper.h Show resolved Hide resolved
src/arkode/arkode_arkstep.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Show resolved Hide resolved
Copy link
Collaborator

@Steven-Roberts Steven-Roberts left a comment

Choose a reason for hiding this comment

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

Initial pass at latest changes

src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
src/sundials/sundials_stepper_impl.h Show resolved Hide resolved
include/sundials/sundials_stepper.h Outdated Show resolved Hide resolved
src/sundials/sundials_stepper.c Outdated Show resolved Hide resolved
include/sundials/sundials_stepper.h Outdated Show resolved Hide resolved
@Steven-Roberts
Copy link
Collaborator

Steven-Roberts commented Aug 21, 2024

After some more thought and discussion with David on IDA and passing around yp, the primary issue I see is keeping it consistent with y. With operator splitting, for example, y jumps around between calls to evolve, and I don't see how to keep yp consistent. With MRI, the change in the forcing between stages can cause a similar issue. This is not something the outer methods can do and needs to be done in the inner SunStepper integrator.

I'd like to propose the following to make SunStepper self starting (all you need is y to reset then evolve):

  • Remove yp as arguments
  • To support an IDA SunStepper, have a constructor like
int IDACreateSUNStepper(void* inner_ida_mem, ConsistentYPrimeFn f, SUNStepper* stepper)

where there's a function which can be specified by users (we could have a default that uses IDACalcIC) to get yp given y.

ARKodeMem ark_mem = (ARKodeMem)arkode_mem;

SUNErrCode err = SUNStepper_Create(ark_mem->sunctx, stepper);
if (err != SUN_SUCCESS) { return ARK_SUNSTEPPER_ERR; }
Copy link
Collaborator

Choose a reason for hiding this comment

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

Should this invoke arkProcessError?

Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't feel strongly on this question either way. While it is an ARKODE-provided routine (and so arkProcessError makes sense), it is a utility routine provided for SUNStepper users, so I could equivalently see it returning a SUNErrCode instead of an int, and "feeling" more similar to the SUNStepper routines.

swig/sundials/fsundials_stepper.i Outdated Show resolved Hide resolved
SUNStepperFullRhsFn fullrhs;
SUNStepperResetFn reset;
SUNStepperSetStopTimeFn setstoptime;
SUNStepperSetForcingFn setforcing;
Copy link
Collaborator

Choose a reason for hiding this comment

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

Note for reviewers: unlike MRIStepInnerStepper, SUNStepper does not include properties for forcing. It defers that to the implementation. This streamlines this struct, removes the need for a forcing "getter," and can reduce storage (saves a N_Vector for ForcingStep over previous version).

.. _SUNStepper.Description.BaseMethods.SteppingFunctions:

Stepping Functions
^^^^^^^^^^^^^^^^^^
Copy link
Collaborator

Choose a reason for hiding this comment

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

Note for reviewers: Unlike MRIStepInnerStepper, the functions in this section are publicly documented. They can be hidden if folks think that's preferable.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I think that these functions should be publicly documented, in case a user wishes to construct a custom SUNStepper. I'm particularly thinking of operator-splitting methods where analytical solutions exist for some partitions, in which case each of these would be easy to implement but may have subtle differences.

Attaching and Accessing the Content Pointer
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. c:function:: SUNErrCode SUNStepper_SetContent(SUNStepper stepper, void *content)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Note for reviewers: these functions are new to SUNStepper compared to MRIStepInnerStepper. Since the SUNStepper struct is private, I think this is the proper way to access the last_flag.

SUNErrCode SUNStepper_SetOneStepFn(SUNStepper stepper, SUNStepperOneStepFn fn);

SUNDIALS_EXPORT
SUNErrCode SUNStepper_SetFullRhsFn(SUNStepper stepper, SUNStepperFullRhsFn fn);
Copy link
Collaborator

Choose a reason for hiding this comment

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

Note to reviewers: there is no analogous SUNStepper_FullRhs. This is following what MRI does where it was presumably excluded because the full RHS functions are private in ARKODE. I can make a public SUNStepper_FullRhs function if folks want.

This section describes the virtual methods defined by the :c:type:`SUNStepper`
abstract base class.

An :c:type:`SUNStepper` *must* provide implementations of the following
Copy link
Collaborator

Choose a reason for hiding this comment

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

I've updated the docs to indicate they are optional and determined by the "consumer" of the SUNStepper. I'll need to specify the required functionality in the operator splitting PR.

Copy link
Collaborator

@drreynolds drreynolds left a comment

Choose a reason for hiding this comment

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

A few minor items, but once those are resolved I think this is good to merge.

.. _SUNStepper.Description.BaseMethods.SteppingFunctions:

Stepping Functions
^^^^^^^^^^^^^^^^^^
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think that these functions should be publicly documented, in case a user wishes to construct a custom SUNStepper. I'm particularly thinking of operator-splitting methods where analytical solutions exist for some partitions, in which case each of these would be easy to implement but may have subtle differences.

Comment on lines 4699 to 4710
.. c:function:: int ARKodeCreateSUNStepper(void *inner_arkode_mem, SUNStepper *stepper)

Wraps an ARKODE integrator as a :c:type:`SUNStepper`.

:param arkode_mem: pointer to the ARKODE memory block.
:param stepper: the :c:type:`SUNStepper` object.

:retval ARK_SUCCESS: the function exited successfully.
:retval ARK_MEM_FAIL: a memory allocation failed.
:retval ARK_SUNSTEPPER_ERR: the :c:type:`SUNStepper` initialization failed.

.. versionadded:: x.y.z
Copy link
Collaborator

Choose a reason for hiding this comment

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

Are there any assumed requirements about the ARKODE solver module? Specifically, since some steppers "support" all ARKODE functionality, whereas other steppers only support a tiny subset, then it would be good to clarify which steppers this function would work for.

Copy link
Collaborator

Choose a reason for hiding this comment

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

The only functionality not supported by all of ARKODE is SUNStepper_SetForcing. I'll make note of that here and in ARKStep, ERKStep, etc.

doc/shared/sunstepper/SUNStepper_Description.rst Outdated Show resolved Hide resolved
doc/shared/sunstepper/SUNStepper_Description.rst Outdated Show resolved Hide resolved
Comment on lines +18 to +19
Implementing a SUNStepper
=========================
Copy link
Collaborator

Choose a reason for hiding this comment

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

It would be great if we had an example that implemented a custom SUNStepper, in which case this file could also mention an example where users can see these instructions "in action." I'm imagining a simple operator-splitting example where one partition has an analytical solution.

Copy link
Collaborator

Choose a reason for hiding this comment

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

In the operator splitting PR, currently only a unit test creates a custom SUNStepper. I can see about modifying the advection-diffusion-reaction example so we have some something to point to.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Or this could be added in the near future (not required for approval of this PR).

ARKodeMem ark_mem = (ARKodeMem)arkode_mem;

SUNErrCode err = SUNStepper_Create(ark_mem->sunctx, stepper);
if (err != SUN_SUCCESS) { return ARK_SUNSTEPPER_ERR; }
Copy link
Collaborator

Choose a reason for hiding this comment

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

I don't feel strongly on this question either way. While it is an ARKODE-provided routine (and so arkProcessError makes sense), it is a utility routine provided for SUNStepper users, so I could equivalently see it returning a SUNErrCode instead of an int, and "feeling" more similar to the SUNStepper routines.

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.

4 participants