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 ARKStep module. Most of the steps are independent of the NVECTOR, SUNMATRIX, SUNLINSOL and SUNNONLINSOL implementations used. For the steps that are not, refer to the sections Vector Data Structures, Matrix Data Structures, Description of the SUNLinearSolver module, and Nonlinear Solver Data Structures 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. Set problem dimensions, etc.

    This generally includes the problem size, N, and may include the local vector length Nlocal.

    Note

    The variables N and Nlocal should be of type sunindextype.

  3. Set vector of initial values

    To set the vector y0 of initial values, use the appropriate functions defined by the particular NVECTOR implementation.

    For native SUNDIALS vector implementations (except the CUDA and RAJA based ones), use a call of the form

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

    if the realtype array ydata containing the initial values of \(y\) already exists. Otherwise, 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);
    

    See the sections The NVECTOR_SERIAL Module through The NVECTOR_PTHREADS Module for details.

    For the HYPRE and PETSc vector wrappers, first create and initialize the underlying vector, and then create the NVECTOR wrapper with a call of the form

    y0 = N_VMake_***(yvec);
    

    where yvec is a HYPRE or PETSc vector. Note that calls like N_VNew_***(...) and N_VGetArrayPointer_***(...) are not available for these vector wrappers. See the sections The NVECTOR_PARHYP Module and The NVECTOR_PETSC Module for details.

    If using either the CUDA- or RAJA-based vector implementations use a call of the form

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

    where c is a pointer to a suncudavec or sunrajavec vector class if this class already exists. Otherwise, create a new vector by making a call of the form

    N_VGetDeviceArrayPointer_***
    

    or

    N_VGetHostArrayPointer_***
    

    Note that the vector class will allocate memory on both the host and device when instantiated. See the sections The NVECTOR_CUDA Module and The NVECTOR_RAJA Module for details.

  4. Create ARKStep object

    Call arkode_mem = ARKStepCreate(...) to create the ARKStep memory block. ARKStepCreate() returns a void* pointer to this memory structure. See the section ARKStep initialization and deallocation functions for details.

  5. Specify integration tolerances

    Call ARKStepSStolerances() or ARKStepSVtolerances() to specify either a scalar relative tolerance and scalar absolute tolerance, or a scalar relative tolerance and a vector of absolute tolerances, respectively. Alternatively, call ARKStepWFtolerances() to specify a function which sets directly the weights used in evaluating WRMS vector norms. See the section ARKStep tolerance specification functions for details.

    If a problem with non-identity mass matrix is used, and the solution units differ considerably from the equation units, absolute tolerances for the equation residuals (nonlinear and linear) may be specified separately through calls to ARKStepResStolerance(), ARKStepResVtolerance(), or ARKStepResFtolerance().

  6. Create matrix object

    If a nonlinear solver requiring a linear solver will be used (e.g., a Newton iteration) and the linear solver will be a matrix-based linear solver, then a template Jacobian matrix must be created by using the appropriate functions defined by the particular SUNMATRIX implementation.

    For the SUNDIALS-supplied SUNMATRIX implementations, the matrix object may be created using a call of the form

    SUNMatrix A = SUNBandMatrix(...);
    

    or

    SUNMatrix A = SUNDenseMatrix(...);
    

    or

    SUNMatrix A = SUNSparseMatrix(...);
    

    Similarly, if the problem involves a non-identity mass matrix, and the mass-matrix linear systems will be solved using a direct linear solver, then a template mass matrix must be created by using the appropriate functions defined by the particular SUNMATRIX implementation.

    NOTE: The dense, banded, and sparse matrix objects are usable only in a serial or threaded environment.

  7. Create linear solver object

    If a nonlinear solver requiring a linear solver will be used (e.g., a Newton iteration), or if the problem involves a non-identity mass matrix, then the desired linear solver object(s) must be created by using the appropriate functions defined by the particular SUNLINSOL implementation.

    For any of the SUNDIALS-supplied SUNLINSOL implementations, the linear solver object may be created using a call of the form

    SUNLinearSolver LS = SUNLinSol_*(...);
    

    where * can be replaced with “Dense”, “SPGMR”, or other options, as discussed in the sections Linear solver interface functions and Description of the SUNLinearSolver module.

  8. Set linear solver optional inputs

    Call *Set* functions from the selected linear solver module to change optional inputs specific to that linear solver. See the documentation for each SUNLINSOL module in the section Description of the SUNLinearSolver module for details.

  9. Attach linear solver module

    If a linear solver was created above for implicit stage solves, initialize the ARKLS linear solver interface by attaching the linear solver object (and Jacobian matrix object, if applicable) with the call (for details see the section Linear solver interface functions):

    ier = ARKStepSetLinearSolver(...);
    

    Similarly, if the problem involves a non-identity mass matrix, initialize the ARKLS mass matrix linear solver interface by attaching the mass linear solver object (and mass matrix object, if applicable) with the call (for details see the section Linear solver interface functions):

    ier = ARKStepSetMassLinearSolver(...);
    
  10. Set optional inputs

    Call ARKStepSet* functions to change any optional inputs that control the behavior of ARKStep from their default values. See the section Optional input functions for details.

  11. Create nonlinear solver object

    If the problem involves an implicit component, and if a non-default nonlinear solver object will be used for implicit stage solves (see the section Nonlinear solver interface functions), then the desired nonlinear solver object must be created by using the appropriate functions defined by the particular SUNNONLINSOL implementation (e.g., NLS = SUNNonlinSol_***(...); where *** is the name of the nonlinear solver (see the section Nonlinear Solver Data Structures for details).

    For the SUNDIALS-supplied SUNNONLINSOL implementations, the nonlinear solver object may be created using a call of the form

    SUNNonlinearSolver NLS = SUNNonlinSol_Newton(...);
    

    or

    SUNNonlinearSolver NLS = SUNNonlinSol_FixedPoint(...);
    
  12. Attach nonlinear solver module

    If a nonlinear solver object was created above, then it must be attached to ARKStep using the call (for details see the section Nonlinear solver interface functions):

    ier = ARKStepSetNonlinearSolver(...);
    
  13. Set nonlinear solver optional inputs

    Call the appropriate set functions for the selected nonlinear solver module to change optional inputs specific to that nonlinear solver. These must be called after attaching the nonlinear solver to ARKStep, otherwise the optional inputs will be overridden by ARKStep defaults. See the section Nonlinear Solver Data Structures for more information on optional inputs.

  14. Specify rootfinding problem

    Optionally, call ARKStepRootInit() to initialize a rootfinding problem to be solved during the integration of the ODE system. See the section Rootfinding initialization function for general details, and the section Optional input functions for relevant optional input calls.

  15. Advance solution in time

    For each point at which output is desired, call

    ier = ARKStepEvolve(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 the section ARKStep solver function for details.

  16. Get optional outputs

    Call ARKStepGet* functions to obtain optional output. See the section Optional output functions for details.

  17. Deallocate memory for solution vector

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

    N_VDestroy(y);
    
  18. Free solver memory

    Call ARKStepFree(&arkode_mem) to free the memory allocated for the ARKStep module (and any nonlinear solver module).

  19. Free linear solver and matrix memory

    Call SUNLinSolFree() and (possibly) SUNMatDestroy() to free any memory allocated for the linear solver and matrix objects created above.

  20. Finalize MPI, if used

    Call MPI_Finalize to terminate MPI.

SUNDIALS provides some linear solvers only as a means for users to get problems running and not as highly efficient solvers. For example, if solving a dense system, we suggest using the LAPACK solvers if the size of the linear system is \(> 50,000\) (thanks to A. Nicolai for his testing and recommendation). The table below shows the linear solver interfaces available as SUNLinearSolver modules and the vector implementations required for use. As an example, one cannot use the dense direct solver interfaces with the MPI-based vector implementation. However, as discussed in section Description of the SUNLinearSolver module the SUNDIALS packages operate on generic SUNLinearSolver objects, allowing a user to develop their own solvers should they so desire.

SUNDIALS linear solver interfaces and vector implementations that can be used for each

Linear Solver Interface Serial Parallel (MPI) OpenMP pThreads hypre Vec. PETSc Vec. CUDA RAJA User Suppl.
Dense X   X X         X
Band X   X X         X
LapackDense X   X X         X
LapackBand X   X X         X
KLU X   X X         X
SuperLU_MT X   X X         X
SPGMR X X X X X X X X X
SPFGMR X X X X X X X X X
SPBCGS X X X X X X X X X
SPTFQMR X X X X X X X X X
PCG X X X X X X X X X
User supplied X X X X X X X X X