SUNStepper basics based on MRIStepInnerStepper #463

balos1 · 2024-04-29T22:37:04Z

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

src/arkode/arkode_arkstep.c

src/sundials/sundials_stepper.c

src/sundials/sundials_stepper_impl.h

balos1 · 2024-05-06T18:40:38Z

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

drreynolds

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

drreynolds · 2024-05-06T22:08:08Z

src/arkode/arkode_arkstep.c

+/* TryStep step result flags */
+#define TRYSTEP_FAILED  -1
+#define TRYSTEP_SUCCESS +0
+#define TRYSTEP_ADAPT   +1
+


Should these be in include/sundials/sundials_stepper.h instead?

drreynolds · 2024-05-06T22:11:23Z

src/arkode/xbraid/arkode_xbraid.c

-  /* Step was successful and passed the error test */
-  *ark_flag = STEP_SUCCESS;
-  return ARK_SUCCESS;
+  return arkStep_TryStep(arkode_mem, tstart, tstop, y, ark_flag);


This is a nice simplification!

drreynolds · 2024-05-06T22:15:36Z

src/sundials/sundials_stepper.c

+  if (!stepper->fused_vectors)
+  {
+    stepper->fused_vectors = N_VNewVectorArray(stepper->nforcing, SUNCTX_);
+    SUNCheckLastErr();
+  }
+
+  if (!stepper->fused_scalars)
+  {
+    stepper->fused_scalars =
+      (sunrealtype*)calloc(count + 1, sizeof(*stepper->fused_scalars));
+    SUNAssert(stepper->fused_scalars, SUN_ERR_MEM_FAIL);
+  }


Shouldn't both of these also check to see whether the length of the already-allocated fused_scalars and fused_vectors arrays have length at least count?

gardner48 · 2024-05-12T00:46:49Z

include/sundials/sundials_stepper.h

+#endif
+
+typedef _SUNDIALS_STRUCT_ SUNStepper_s* SUNStepper;
+


@balos1 what is stop_reason needed for? Wtih OneStep and Advance is it to distinguish between a successful return and a stop time return?

Yes, or a root return.

Do you need to differentiate between between a tstop return and a root found return or is a positive return sufficient to indicate a successful return that's different from a normal successful (zero) return?

I could imagine (but dont have a concrete example at the moment) that a root return might need special handling when checkpointing and/or during backwards integration.

It would be nice it this could just encoded in the return value but having a SUNErrCode also indicate a status would be mixing purposes.

Indeed, I think its best to have the error code and status separated.

gardner48 · 2024-05-12T00:50:13Z

include/sundials/sundials_stepper.h

+
+typedef int (*SUNStepperTryStepFn)(SUNStepper stepper, sunrealtype t0,
+                                   sunrealtype tout, N_Vector y,
+                                   sunrealtype* tret, int* stop_reason);


What does everyone think about removing t0 as an input and haveing y only be for output? That is, as stepper should always advance from it's last reached solution after a successful return unless something like Reset or ReInit is called.

Similarly, could tout be removed from TryStep and OneStep?

How would someone set the t0 to begin with? Are you assuming this would be done with the package create/init routines before it is wrapped as a stepper?

I would expect this to be done with the stepper Create or Init function or alternatively Reset could be called in the first step if we don't require a stepper to be initialized with the starting state.

I would expect this to be done with the stepper Create or Init

Do you mean the package "steppers"? The SUNStepper create does not currently have t0 (but it could) and it has no Init function.

alternatively Reset could be called in the first step if we don't require a stepper to be initialized with the starting state

So we would assume tR = t0 = 0 then?

Do you mean the package "steppers"? The SUNStepper create does not currently have t0 (but it could) and it has no Init function.

Yes, I was thinking the constructor for an implementation of a SUNStepper (i.e., MyStepper(t0, y0, ctx)) could take the initial condition. The package steppers are a little different since they wrap an integrator that was already created with the initial condition set. Alternatively having SUNStepper_Init(t0, y0) of SUNStepper_Init(t0, y0, yp0) might not be a bad idea.

So we would assume tR = t0 = 0 then?

Just tR = t0 since t0 is not necessarily 0

Just tR = t0 since t0 is not necessarily 0

Well, yes, but my point was that if there was no Init that took t0 (or Create) then we would not know t0.

The package steppers are a little different since they wrap an integrator that was already created with the initial condition set.

CVODE/IDA do not take the initial condition and t0 until their Init routines are called. So you could have the case where a stepper is created with a integrator that does not yet have the IC/t0 set. With that in mind, I think having SUNStepper_Init routine is cleaner.

CVODE/IDA do not take the initial condition and t0 until their Init routines are called. So you could have the case where a stepper is created with a integrator that does not yet have the IC/t0 set. With that in mind, I think having SUNStepper_Init routine is cleaner.

I think have an Init function is reasonable. As far as wrapping CVODE or IDA as a stepper, creating a stepper from a CVODE / IDA instance that was not initialized would be a user-error as the RHS / RES function would also not be set.

I think that David's alternate option of each consuming package calling Reset is perhaps the most appropriate. This already occurs in MRIStep (or it would take a minimal change to do so), and I anticipate that the operator-splitting stepper would also want to control the initial condition for any sub-stepper.

gardner48 · 2024-05-12T00:52:29Z

include/sundials/sundials_stepper.h

+#endif
+
+typedef _SUNDIALS_STRUCT_ SUNStepper_s* SUNStepper;
+


The functions below return an int but the corresponding SUNStepper_* functions return a SUNErrorCode. Should these have the same return type?

gardner48 · 2024-05-12T02:42:22Z

include/sundials/sundials_stepper.h

+#endif
+
+typedef _SUNDIALS_STRUCT_ SUNStepper_s* SUNStepper;
+


Should the SUNStepper support DAEs i.e., should there be y_ret and yp_ret for Advance, OneStep and TryStep? I'm thinking there are three ways to handle this:

Have SUNODEStepper and SUNDAEStepper with different versions of these functions

Have ODE and DAE versions of the functions i.e., ODETryStep and DAETryStep with different signatures

Add yp_ret to each of the functions as an optional output i.e., allow passing NULL for yp_ret

I'm thinking 3 might be the best approach.

I also like 3.

I like option 3

gardner48 · 2024-05-12T02:49:31Z

src/arkode/arkode_arkstep.c

+  retval = SUNStepper_SetContent(*stepper, inner_arkode_mem);
+  if (retval != ARK_SUCCESS) { return (retval); }
+
+  retval = SUNStepper_SetAdvanceFn(*stepper, arkStep_SUNStepperAdvance);


Attach OneStepFn and TryStepFn also?

gardner48 · 2024-05-12T03:10:29Z

src/arkode/arkode_arkstep.c

+ * Attempts one internal time step using ARKStepEvolve
+ * ---------------------------------------------------------------------------*/
+
+int arkStep_TryStep(void* arkode_mem, sunrealtype tstart, sunrealtype tstop,


The body of this function was taken from ARKBraid_TakeStep and some of the behavior is specific to XBraid's use (having to reset before every step and use a specific step size) and I'm not sure if it's how we want to SUNStepper to behave. This is related to my question above about having a stepper advance from its last reached state.

If TryStep is called in a loop (i.e., SUNStepper could provide an implementation of Advance that calls TryStep as needed) the call to Reset could lead to unnecessary RHS evaluations.

Do we want the step size to be determined by tstart and tstop or internally by the stepper?

gardner48 · 2024-05-12T03:18:00Z

src/sundials/sundials_stepper.c

+  *stepper = (SUNStepper)malloc(sizeof(**stepper));
+  SUNAssert(stepper, SUN_ERR_MALLOC_FAIL);
+
+  memset(*stepper, 0, sizeof(**stepper));


Could make the malloc above a calloc and remove the memset

gardner48 · 2024-05-12T03:30:42Z

src/sundials/sundials_stepper.c

+  (*stepper)->ops = (SUNStepper_Ops)malloc(sizeof(*((*stepper)->ops)));
+  SUNAssert((*stepper)->ops, SUN_ERR_MALLOC_FAIL);
+
+  memset((*stepper)->ops, 0, sizeof(*((*stepper)->ops)));


Could make the malloc above a calloc and remove the memset

gardner48 · 2024-05-12T03:39:34Z

src/sundials/sundials_stepper.c

+  return SUN_SUCCESS;
+}
+
+SUNErrCode SUNStepper_Advance(SUNStepper stepper, sunrealtype t0,


If TryStep advances from the last time reached and if we add a SetStopTime function, there could be a default implementation of Advance that calls TryStep

gardner48 · 2024-05-12T03:40:39Z

src/sundials/sundials_stepper.c

+  return SUN_ERR_NOT_IMPLEMENTED;
+}
+
+SUNErrCode SUNStepper_Step(SUNStepper stepper, sunrealtype t0, sunrealtype tout,


Like Advance we could have a default implementation that calls TryStep

src/arkode/arkode_arkstep.c

Steven-Roberts · 2024-05-13T16:11:57Z

include/sundials/sundials_stepper.h

+#endif
+
+typedef _SUNDIALS_STRUCT_ SUNStepper_s* SUNStepper;
+


I also like 3.

Steven-Roberts · 2024-05-13T16:18:27Z

src/arkode/arkode_arkstep.c

+  if (retval != ARK_SUCCESS) { return (retval); }
+
+  /* set the stop time */
+  retval = ARKStepSetStopTime(arkode_mem, tout);


Is this needed? Maybe it could be an option. For operator splitting, it would be preferable not to have this.

In our team discussions we defined the Advance function to operate in this way. That is, we said Advance would step to precisely tout as opposed to Evolve which would potentially step past tout and interpolate back. MRIStep needs this. It sounds like you need an Evolve function for the operator splitting stepper.

Ah I didn't realize we would have both Advance and Evolve. I that case, this is not an issue.

In the last discussion we dropped Evolve because it seemed like we did not need it. @Steven-Roberts for the splitting methods, would Advance still be sufficient or do you need to step past the time the method is supposed to integrate to?

Advance alone would be sufficient, but then I don't see a way to handle symmetric operator splitting methods most efficiently. As an example of the functionality needed, consider integration a stepper from t=0.5 to t=1.5 with output needed at t=1. Two calls to Evolve would be ideal for this, while Advance might add more function evaluations by forcing the first call to end exactly at t=1. Originally, I thought adding an interpolation function to the stepper was the way to fix this, but on second thought I don't think it's reasonable to have interpolation data across the entire tspan passed to Advance.

Right, it's sounding like we need Evolve + SetStopTime rather than Advance or Advance with multiple output times.

Originally we talked about having four functions:

Evolve: step to tout, possibly passing it and returning interpolated solution at t=tout

Advance: step to tout, stopping exactly at tout and return

Step: take one successful step (or fail trying to) and return

TryStep: take exactly one internal step and return

We said we did not have a use for Evolve yet, so I removed it from the SUNStepper API for now. However, it seems the operator splitting stepper does necessitate an Evolve function PLUS a SUNStepper_SetTStop function.

If we add a SUNStepper_SetTStop function, then it should only be used with Evolve, NOT Advance otherwise there is no reason to have Advance as a function since it would work just like Evolve. Furthermore, I believe MRIStep relies on Advance working as defined above (if it does not, then I further question why we need Advance at all).

An alternative to adding SUNStepper_SetTStop+Evolve is to let Evolve and/or Advance return outputs at multiple time (as David suggested). This is different than any of our packages (as Steven points out), but it is also something that pretty much every other time integration package lets you do.

I'm leaning toward Evolve + SetStopTime since Advance (or even Evolve) with multiple outputs could just be a wrapper looping over the output times and calling the one output Evolve with the right output vector.

However, it seems the operator splitting stepper does necessitate an Evolve function PLUS a SUNStepper_SetTStop function.

It does need the Evolve function. SetTStop would only be needed if Advance was removed (as it could be achieved by combining Evolve and SetTStop).

So I think the question is which set of functions to use

Evolve, Advance, Step, TryStep or

Evolve, SetTStop, Step, TryStep

I'm a bit late to the party, but it sounds to me like Evolve + SetStopTime is the way to go. I'll note that like Steven, in the MERK multirate time-stepping methods I evolve over a longer duration, periodically storing the solution for some intermediate stages. It would be perfectly handled by this approach, whereas right now I anticipate that the MRIStepInnerStepper may take a small number of extra steps as it stops to return intermediate solutions. I'll defer to the rest of the group on the need for Step and TryStep usage modes.

Steven-Roberts · 2024-05-13T16:29:33Z

src/sundials/sundials_stepper.c

+{
+  SUNFunctionBegin(sunctx);
+
+  *stepper = NULL;


Unnecessary

Steven-Roberts · 2024-05-13T16:35:55Z

src/sundials/sundials_stepper.c

+  sunrealtype tau, taui;
+  int i;


These seem unnecessary

Yeah they are leftover from copy/pasting from MRIStepInnerStepper. I will remove them.

Steven-Roberts · 2024-05-13T16:43:37Z

src/sundials/sundials_stepper.c

+                              sunrealtype tout, N_Vector y, sunrealtype* tret,
+                              int* stop_reason)
+{
+  SUNFunctionBegin(stepper->sunctx);


Do all the functions in this file need to check if stepper is NULL?

No, we do not need to check that. If the stepper is NULL we cannot handle the error properly anyways.

balos1 added 2 commits April 29, 2024 14:12

move MRIStepInnerStepper to SUNStepper

8b8b7ef

add advance, onestep, trystep functions

c31bf03

balos1 assigned Steven-Roberts Apr 29, 2024

drreynolds reviewed Apr 30, 2024

View reviewed changes

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 added 5 commits May 6, 2024 10:27

revert arkode changes

7eaf355

implement Advance, OneStep, and TryStep

381a90c

cleanup

536c14b

wrap SUNStepper as MRIStepInnerStepper

9cc9567

Merge remote-tracking branch 'origin/develop' into feature/sunstepper

0d3e9b0

balos1 changed the title ~~move MRIStepInnerStepper to SUNStepper~~ SUNStepper basics based on MRIStepInnerStepper May 6, 2024

drreynolds reviewed May 6, 2024

View reviewed changes

update sunstepper evolve functions to pass tret and return reason

f97e71b

gardner48 reviewed May 12, 2024

View reviewed changes

Steven-Roberts reviewed May 13, 2024

View reviewed changes

SUNStepper basics based on MRIStepInnerStepper #463

Are you sure you want to change the base?

SUNStepper basics based on MRIStepInnerStepper #463

Conversation

balos1 commented Apr 29, 2024

balos1 commented May 6, 2024

drreynolds left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gardner48 May 13, 2024 • edited

Choose a reason for hiding this comment

balos1 May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Steven-Roberts May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 May 13, 2024 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

balos1 May 13, 2024 •

edited

balos1 May 13, 2024 •

edited

balos1 May 13, 2024 •

edited

gardner48 May 13, 2024 •

edited

balos1 May 13, 2024 •

edited

Steven-Roberts May 13, 2024 •

edited

balos1 May 13, 2024 •

edited