3.4.4.1. A skeleton of the user’s main program

The following is a skeleton of the user’s main program (or calling program) for the integration of an ODE IVP using the SPRKStep module. Most of the steps are independent of the NVECTOR implementation used. For the steps that are not, refer to §9 for the specific name of the function to be called or macro to be referenced.

  1. Initialize parallel or multi-threaded environment, if appropriate.

    For example, call MPI_Init to initialize MPI if used, or set num_threads, the number of threads to use within the threaded vector functions, if used.

  2. Create the SUNDIALS simulation context object.

    Call SUNContext_Create() to allocate the SUNContext object.

  3. Set problem dimensions, etc.

    This generally includes the problem size, N, and may include the local vector length Nlocal. The problem size N is the size including both the q and p variables.

    Note

    The variables N and Nlocal should be of type sunindextype.

  4. Set vector of initial values

    To set the vector y0 of initial values, use the appropriate functions defined by the particular NVECTOR implementation. The vector should include both the q and p variables.

    For most native SUNDIALS vector implementations, use a call of the form

    y0 = N_VMake_***(..., ydata);
    

    if the sunrealtype array ydata containing the initial values of \(y\) already exists. For some GPU-enabled vectors, a similar constructor can be used to provide host and device data pointers. If the data array does not already exist, create a new vector by making a call of the form

    y0 = N_VNew_***(...);
    

    and then set its elements by accessing the underlying data where it is located with a call of the form

    ydata = N_VGetArrayPointer_***(y0);
    

    For details on each of SUNDIALS’ provided vector implementations, see the corresponding sections in §9 for details.

  5. Create SPRKStep object

    Call arkode_mem = SPRKStepCreate(...) to create the SPRKStep memory block. SPRKStepCreate() returns a void* pointer to this memory structure. See §3.4.4.2.1 for details.

  6. Specify time step size

    Call SPRKStepSetFixedStep() to set the fixed time step size. .. or SPRKStepAdaptivityFn() to specify either a fixed time-step .. size or a callback function that adapts the time-step size. SPRKStep .. does not support error-based adaptivity like other ARKODE time-stepper .. modules due to the incompatibility with symplectic methods.

  7. Set optional inputs

    Call SPRKStepSet* functions to change any optional inputs that control the behavior of SPRKStep from their default values. See §3.4.4.2.4 for details.

  8. Specify rootfinding problem

    Optionally, call SPRKStepRootInit() to initialize a rootfinding problem to be solved during the integration of the ODE system. See §3.4.4.2.2 for general details, and §3.4.4.2.4 for relevant optional input calls.

  9. Advance solution in time

    For each point at which output is desired, call

    ier = SPRKStepEvolve(arkode_mem, tout, yout, &tret, itask);
    

    Here, itask specifies the return mode. The vector yout (which can be the same as the vector y0 above) will contain \(y(t_\text{out})\). See §3.4.4.2.3 for details.

  10. Get optional outputs

    Call SPRKStepGet* functions to obtain optional output. See §3.4.4.2.6 for details.

  11. Deallocate memory for solution vector

    Upon completion of the integration, deallocate memory for the vector y (or yout) by calling the NVECTOR destructor function:

    N_VDestroy(y);
    
  12. Free solver memory

    Call SPRKStepFree() to free the memory allocated for the SPRKStep module.

  13. Finalize MPI, if used

    Call MPI_Finalize to terminate MPI.