Skip to content

Commit

Permalink
implement Advance, OneStep, and TryStep
Browse files Browse the repository at this point in the history
  • Loading branch information
@balos1
balos1 committed May 6, 2024
1 parent 7eaf355 commit 381a90c
Show file tree
Hide file tree
Showing 9 changed files with 349 additions and 131 deletions.
5 changes: 5 additions & 0 deletions include/arkode/arkode_arkstep.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include <arkode/arkode_ls.h>
#include <sunadaptcontroller/sunadaptcontroller_imexgus.h>
#include <sunadaptcontroller/sunadaptcontroller_soderlind.h>
#include <sundials/sundials_stepper.h>

#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
Expand Down Expand Up @@ -339,6 +340,10 @@ SUNDIALS_EXPORT void ARKStepPrintMem(void* arkode_mem, FILE* outfile);
SUNDIALS_EXPORT int ARKStepCreateMRIStepInnerStepper(void* arkode_mem,
MRIStepInnerStepper* stepper);

/* SUNStepper interface functions */
SUNDIALS_EXPORT int ARKStepCreateSUNStepper(void* arkode_mem,
SUNStepper* stepper);

/* Relaxation functions */
SUNDIALS_EXPORT int ARKStepSetRelaxFn(void* arkode_mem, ARKRelaxFn rfn,
ARKRelaxJacFn rjac);
Expand Down
4 changes: 4 additions & 0 deletions include/arkode/arkode_mristep.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#include <arkode/arkode_butcher_dirk.h>
#include <arkode/arkode_butcher_erk.h>
#include <arkode/arkode_ls.h>
#include <sundials/sundials_stepper.h>

#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
Expand Down Expand Up @@ -304,6 +305,9 @@ SUNDIALS_EXPORT void MRIStepPrintMem(void* arkode_mem, FILE* outfile);
SUNDIALS_EXPORT int MRIStepInnerStepper_Create(SUNContext sunctx,
MRIStepInnerStepper* stepper);

SUNDIALS_EXPORT int MRIStepInnerStepper_CreateFromSUNStepper(
SUNContext sunctx, SUNStepper sunstepper, MRIStepInnerStepper* stepper);

SUNDIALS_EXPORT int MRIStepInnerStepper_Free(MRIStepInnerStepper* stepper);

SUNDIALS_EXPORT int MRIStepInnerStepper_SetContent(MRIStepInnerStepper stepper,
Expand Down
8 changes: 1 addition & 7 deletions include/sundials/sundials_stepper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,17 +21,14 @@ extern "C" {

typedef _SUNDIALS_STRUCT_ SUNStepper_s* SUNStepper;

typedef int (*SUNStepperEvolveFn)(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);

typedef int (*SUNStepperAdvanceFn)(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);

typedef int (*SUNStepperOneStepFn)(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);

typedef int (*SUNStepperTryStepFn)(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);
sunrealtype tout, N_Vector y, int* ark_flag);

typedef int (*SUNStepperFullRhsFn)(SUNStepper stepper, sunrealtype t,
N_Vector y, N_Vector f, int mode);
Expand All @@ -50,9 +47,6 @@ SUNErrCode SUNStepper_SetContent(SUNStepper stepper, void* content);
SUNDIALS_EXPORT
SUNErrCode SUNStepper_GetContent(SUNStepper stepper, void** content);

SUNDIALS_EXPORT
SUNErrCode SUNStepper_SetEvolveFn(SUNStepper stepper, SUNStepperEvolveFn fn);

SUNDIALS_EXPORT
SUNErrCode SUNStepper_SetAdvanceFn(SUNStepper stepper, SUNStepperAdvanceFn fn);

Expand Down
261 changes: 261 additions & 0 deletions src/arkode/arkode_arkstep.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,19 @@
#include <sundials/sundials_math.h>
#include <sunnonlinsol/sunnonlinsol_newton.h>

#include "arkode/arkode.h"
#include "arkode/arkode_butcher.h"
#include "arkode_arkstep_impl.h"
#include "arkode_impl.h"
#include "arkode_interp_impl.h"

#define FIXED_LIN_TOL

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

/*===============================================================
ARKStep Exported functions -- Required
===============================================================*/
Expand Down Expand Up @@ -3219,6 +3225,203 @@ int arkStep_ComputeSolutions_MassFixed(ARKodeMem ark_mem, sunrealtype* dsmPtr)
return (ARK_SUCCESS);
}

/*---------------------------------------------------------------
Utility routines for interfacing with SUNStepper
---------------------------------------------------------------*/

int ARKStepCreateSUNStepper(void* inner_arkode_mem, SUNStepper* stepper)
{
int retval;
ARKodeMem ark_mem;
ARKodeARKStepMem step_mem;

retval = arkStep_AccessStepMem(inner_arkode_mem, "ARKStepCreateSUNStepper",
&ark_mem, &step_mem);
if (retval)
{
arkProcessError(NULL, ARK_ILL_INPUT, __LINE__, __func__, __FILE__,
"The ARKStep memory pointer is NULL");
return ARK_ILL_INPUT;
}

retval = SUNStepper_Create(ark_mem->sunctx, stepper);
if (retval != ARK_SUCCESS) { return (retval); }

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

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

retval = SUNStepper_SetFullRhsFn(*stepper, arkStep_SUNStepperFullRhs);
if (retval != ARK_SUCCESS) { return (retval); }

retval = SUNStepper_SetResetFn(*stepper, arkStep_SUNStepperReset);
if (retval != ARK_SUCCESS) { return (retval); }

return (ARK_SUCCESS);
}

/*------------------------------------------------------------------------------
arkStep_SUNStepperAdvance
Implementation of SUNStepperAdvanceFn to advance the inner (fast)
ODE IVP.
----------------------------------------------------------------------------*/

int arkStep_SUNStepperAdvance(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y)
{
void* arkode_mem; /* arkode memory */
sunrealtype tret; /* return time */
sunrealtype tshift, tscale; /* time normalization values */
N_Vector* forcing; /* forcing vectors */
int nforcing; /* number of forcing vectors */
int retval; /* return value */

/* extract the ARKODE memory struct */
retval = SUNStepper_GetContent(stepper, &arkode_mem);
if (retval != ARK_SUCCESS) { return (retval); }

/* get the forcing data */
retval = SUNStepper_GetForcingData(stepper, &tshift, &tscale, &forcing,
&nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

/* set the inner forcing data */
retval = arkStep_SetInnerForcing(arkode_mem, tshift, tscale, forcing, nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

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

/* evolve inner ODE */
retval = ARKStepEvolve(arkode_mem, tout, y, &tret, ARK_NORMAL);
if (retval < 0) { return (retval); }

/* disable inner forcing */
retval = arkStep_SetInnerForcing(arkode_mem, ZERO, ONE, NULL, 0);
if (retval != ARK_SUCCESS) { return (retval); }

return (ARK_SUCCESS);
}

int arkStep_SUNStepperOneStep(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y)
{
void* arkode_mem; /* arkode memory */
sunrealtype tret; /* return time */
sunrealtype tshift, tscale; /* time normalization values */
N_Vector* forcing; /* forcing vectors */
int nforcing; /* number of forcing vectors */
int retval; /* return value */

/* extract the ARKODE memory struct */
retval = SUNStepper_GetContent(stepper, &arkode_mem);
if (retval != ARK_SUCCESS) { return (retval); }

/* get the forcing data */
retval = SUNStepper_GetForcingData(stepper, &tshift, &tscale, &forcing,
&nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

/* set the inner forcing data */
retval = arkStep_SetInnerForcing(arkode_mem, tshift, tscale, forcing, nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

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

/* evolve inner ODE */
retval = ARKStepEvolve(arkode_mem, tout, y, &tret, ARK_ONE_STEP);
if (retval < 0) { return (retval); }

/* disable inner forcing */
retval = arkStep_SetInnerForcing(arkode_mem, ZERO, ONE, NULL, 0);
if (retval != ARK_SUCCESS) { return (retval); }

return (ARK_SUCCESS);
}

int arkStep_SUNStepperTryStep(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y, int* ark_flag)
{
void* arkode_mem; /* arkode memory */
sunrealtype tret; /* return time */
sunrealtype tshift, tscale; /* time normalization values */
N_Vector* forcing; /* forcing vectors */
int nforcing; /* number of forcing vectors */
int retval; /* return value */

/* extract the ARKODE memory struct */
retval = SUNStepper_GetContent(stepper, &arkode_mem);
if (retval != ARK_SUCCESS) { return (retval); }

/* get the forcing data */
retval = SUNStepper_GetForcingData(stepper, &tshift, &tscale, &forcing,
&nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

/* set the inner forcing data */
retval = arkStep_SetInnerForcing(arkode_mem, tshift, tscale, forcing, nforcing);
if (retval != ARK_SUCCESS) { return (retval); }

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

/* try to evolve inner ODE */
retval = arkStep_TryStep(arkode_mem, t0, tout, y, ark_flag);
if (retval != ARK_SUCCESS) { return (retval); }

/* disable inner forcing */
retval = arkStep_SetInnerForcing(arkode_mem, ZERO, ONE, NULL, 0);
if (retval != ARK_SUCCESS) { return (retval); }

return ARK_SUCCESS;
}

/*------------------------------------------------------------------------------
arkStep_SUNStepperFullRhs
Implementation of SUNStepperFullRhsFn to compute the full inner
(fast) ODE IVP RHS.
----------------------------------------------------------------------------*/

int arkStep_SUNStepperFullRhs(SUNStepper stepper, sunrealtype t, N_Vector y,
N_Vector f, int mode)
{
void* arkode_mem;
int retval;

/* extract the ARKODE memory struct */
retval = SUNStepper_GetContent(stepper, &arkode_mem);
if (retval != ARK_SUCCESS) { return (retval); }

return (arkStep_FullRHS(arkode_mem, t, y, f, mode));
}

/*------------------------------------------------------------------------------
arkStep_SUNStepperReset
Implementation of SUNStepperResetFn to reset the inner (fast) stepper
state.
----------------------------------------------------------------------------*/

int arkStep_SUNStepperReset(SUNStepper stepper, sunrealtype tR, N_Vector yR)
{
void* arkode_mem;
int retval;

/* extract the ARKODE memory struct */
retval = SUNStepper_GetContent(stepper, &arkode_mem);
if (retval != ARK_SUCCESS) { return (retval); }

return (ARKStepReset(arkode_mem, tR, yR));
}

/*---------------------------------------------------------------
Utility routines for interfacing with MRIStep
---------------------------------------------------------------*/
Expand Down Expand Up @@ -3654,6 +3857,64 @@ int arkStep_GetOrder(ARKodeMem ark_mem)
return step_mem->q;
}

/* -----------------------------------------------------------------------------
* arkStep_TryStep
*
* Attempts one internal time step using ARKStepEvolve
* ---------------------------------------------------------------------------*/

int arkStep_TryStep(void* arkode_mem, sunrealtype tstart, sunrealtype tstop,
N_Vector y, int* ark_flag)
{
int flag; /* generic return flag */
int tmp_flag; /* evolve return flag */
sunrealtype tret; /* return time */

/* Check inputs */
if (arkode_mem == NULL) { return ARK_MEM_NULL; }
if (y == NULL) { return ARK_ILL_INPUT; }

/* Reset ARKStep state */
flag = ARKStepReset(arkode_mem, tstart, y);
if (flag != ARK_SUCCESS) { return flag; }

/* Set the time step size */
flag = ARKStepSetInitStep(arkode_mem, tstop - tstart);
if (flag != ARK_SUCCESS) { return flag; }

/* Ignore temporal error test result and force step to pass */
flag = arkSetForcePass(arkode_mem, SUNTRUE);
if (flag != ARK_SUCCESS) { return flag; }

/* Take step, check flag below */
tmp_flag = ARKStepEvolve(arkode_mem, tstop, y, &tret, ARK_ONE_STEP);

/* Re-enable temporal error test check */
flag = arkSetForcePass(arkode_mem, SUNFALSE);
if (flag != ARK_SUCCESS) { return flag; }

/* Check if evolve call failed */
if (tmp_flag < 0)
{
*ark_flag = TRYSTEP_FAILED;
return ARK_SUCCESS;
}

/* Check if temporal error test failed */
flag = arkGetLastKFlag(arkode_mem, &tmp_flag);
if (flag != ARK_SUCCESS) { return flag; }

if (tmp_flag > 0)
{
*ark_flag = TRYSTEP_ADAPT;
return ARK_SUCCESS;
}

/* Step was successful and passed the error test */
*ark_flag = TRYSTEP_SUCCESS;
return ARK_SUCCESS;
}

/*===============================================================
EOF
===============================================================*/
15 changes: 15 additions & 0 deletions src/arkode/arkode_arkstep_impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,6 +215,17 @@ int arkStep_NlsLSolve(N_Vector delta, void* arkode_mem);
int arkStep_NlsConvTest(SUNNonlinearSolver NLS, N_Vector y, N_Vector del,
sunrealtype tol, N_Vector ewt, void* arkode_mem);

/* private functions for interfacing with SUNStepper */
int arkStep_SUNStepperAdvance(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);
int arkStep_SUNStepperOneStep(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y);
int arkStep_SUNStepperTryStep(SUNStepper stepper, sunrealtype t0,
sunrealtype tout, N_Vector y, int* ark_flag);
int arkStep_SUNStepperFullRhs(SUNStepper stepper, sunrealtype t, N_Vector y,
N_Vector f, int mode);
int arkStep_SUNStepperReset(SUNStepper stepper, sunrealtype tR, N_Vector yR);

/* private functions for interfacing with MRIStep */
int arkStep_SetInnerForcing(void* arkode_mem, sunrealtype tshift,
sunrealtype tscale, N_Vector* f, int nvecs);
Expand All @@ -230,6 +241,10 @@ int arkStep_RelaxDeltaE(ARKodeMem ark_mem, ARKRelaxJacFn relax_jac_fn,
long int* relax_jac_fn_evals, sunrealtype* delta_e_out);
int arkStep_GetOrder(ARKodeMem ark_mem);

/* private utility functions */
int arkStep_TryStep(void* arkode_mem, sunrealtype tstart, sunrealtype tstop,
N_Vector y, int* ark_flag);

/*===============================================================
Reusable ARKStep Error Messages
===============================================================*/
Expand Down
Loading

0 comments on commit 381a90c

Please sign in to comment.