Contents

 

Installation

I got the tarball. What do I do next?

Let's assume you downloaded sundials-x.y.z.tar.gz. First, you'll have to uncompress and unpack it:

% tar -xzf sundials-x.y.z.tar

This will create a directory sundials-x.y.z, which we will call srcdir. (Note that if you only downloaded an individual solver tarball, call it solver-x.y.z.tar.gz, uncompressing and unpacking it will result in the directory solver-x.y.z.)

Next, you need to configure and build the SUNDIALS libraries. You can get some information on how to do this here, but more details are available in the user guides (e.g. sundials/cvode/doc/cv_guide.pdf).

Starting with version 2.6.0 of sundials, CMake is the only supported method of installation.

The sundials build process requires CMake version 2.8.1 or higher and a working compiler. On Unix-like operating systems, it also requires Make (and curses, including its development libraries, for the GUI front end to CMake, ccmake), while on Windows it requires Visual Studio. While many Linux distributions offer CMake, the version included is probably out of date. Many new CMake features have been added recently, and you should download the latest version from http://www.cmake.org. Build instructions for CMake (only necessary for Unix-like systems) can be found on the CMake website. Once CMake is installed, Linux/Unix users will be able to use ccmake, while Windows users
will be able to use CMakeSetup.

The default CMake configuration will build all included solvers and associated examples and will build static and shared libraries. The instdir defaults to /usr/local and can be changed by setting the CMAKE INSTALL PREFIX variable. Support for FORTRAN and all other options are disabled. CMake can be used from the command line with the cmake command, or from a curses-based GUI by using the ccmake command. For the examples shown it is assumed that there is a top level sundials directory with appropriate source, build and install directories:

  % mkdir (...)sundials/instdir

  % mkdir (...)sundials/builddir

  % cd (...)sundials/builddir

CMake can be used with a GUI or from the command line. Using CMake from the command line is simply a matter of specifying CMake variable settings with the cmake command. The following will build the default configuration:

 % cmake -DCMAKE_INSTALL_PREFIX=/home/myname/sundials/instdir \ 

  -DEXAMPLES_INSTALL_PATH=/home/myname/sundials/instdir/examples \

  ../srcdir

  % make

  % make install

A complete list of all available options, plus details on using the CMake GUI, can be found in the installation chapter in each of the SUNDIALS user guides.

Can I use a non-default compiler to build SUNDIALS?

Yes, specific compilers can be specified on the CMake command line. The following example specifies gcc and g++ for the CXX compiler,

   
    % CC=gcc CXX=g++ cmake \

    -DCXX_ENABLE=TRUE \

    ../srcdir


This also works for the CMake Gui:
 

% CC=gcc CXX=g++ ccmake ../srcdir


NOTE: If you wish to change the compiler after already configuring, then you must start with a fresh build directory and specify the compiler on the command line as shown.

How do I install SUNDIALS under Windows?

One way of obtaining Windows libraries for the SUNDIALS solvers is to use cygwin (www.cygwin.com), in which case the installation procedure is the same as for any other *NIX system.

To build SUNDIALS for use with Visual Studio the following steps should be performed:

  1. Unzip the downloaded tar file(s) into a directory. This will be the srcdir
  2. Create a separate builddir
  3. Open a Visual Studio Command Prompt and cd to builddir
  4. Run cmake-gui ../srcdir
  5. Hit Configure
  6. Check/Uncheck solvers to be built
  7. Change CMAKE_INSTALL_PREFIX to instdir
  8. Set other options as desired
  9. Hit Generate

Back in the VS Command Window:

  1. Run msbuild ALL_BUILD.vcxproj
  2. Run msbuild INSTALL.vcxproj

The resulting libraries will be in the instdir. The sundials project can also now be opened in Visual Studio. Double click on the ALL_BUILD.vcxproj file to open the project. Build the whole solution to create the sundials libraries. To use the sundials libraries in your own projects, you must set the include directories for your project, add the sundials libraries to your project solution, and set the sundials libraries as dependencies for your project.

Everything installed fine! How do I link the SUNDIALS libraries to my own application?

Each of the sundials solvers is distributed with a set of examples demonstrating basic usage. To build and install the examples, set both EXAMPLES ENABLE and EXAMPLES INSTALL to ON and specify the example installation directory EXAMPLES INSTALL PATH. CMake will generate a CMakeLists.txt configuration file (and Makefile files if on Linux/Unix) that reference the installed sundials headers and libraries. Either the CMakeLists.txt file or the traditional Makefile may be used to build the examples as well as serve as a template for creating user developed solutions. To use the supplied Makefile simply run make to compile and generate the executables. To use CMake, from within the installed example directory, run cmake (or ccmake to use the GUI) followed by make to compile the example code. Note that if CMake is used, it will overwrite the traditional Makefile with a new CMake generated Makefile.

CVODE / CVODES / ARKode

While written in terms of CVODE naming conventions, this section applies almost entirely also to ARKode.

How do I choose tolerances?

For many users, the appropriate choices for tolerance values are a concern. The following pieces of advice are relevant.

  1. The scalar relative tolerance reltol is to be set to control relative errors. So reltol = 1.0E-4 means that errors are controlled to .01%. We do not recommend using reltol larger than 1.0E-3. On the other hand, reltol should not be so small that it is comparable to the unit roundoff of the machine arithmetic (generally around 1.0E-15).
  2. The absolute tolerances abstol (whether scalar or vector) need to be set to control absolute errors when any components of the solution vector y may be so small that pure relative error control is meaningless. For example, if y[i] starts at some nonzero value, but in time decays to zero, then pure relative error control on y[i] makes no sense (and is overly costly) after y[i] is below some noise level. Then abstol (if scalar) or abstol[i] (if a vector) needs to be set to that noise level. If the different components have different noise levels, then abstol should be a vector. See the example cvdenx in the CVODE package, and the discussion of it in the CVODE Examples document. In that problem, the three components vary between 0 and 1, and have different noise levels; hence the abstol vector. It is impossible to give any general advice on abstol values, because the appropriate noise levels are completely problem-dependent. The user or modeler hopefully has some idea as to what those noise levels are.
  3. Finally, it is important to pick all the tolerance values conservatively, because they control the error committed on each individual time step. The final (global) errors are some sort of accumulation of those per-step errors. A good rule of thumb is to reduce the tolerances by a factor of .01 from the actual desired limits on errors. So if you want .01% accuracy (globally), a good choice is reltol = 1.0E-6. But in any case, it is a good idea to do a few experiments with the tolerances to see how the computed solution values vary as tolerances are reduced.

How do I choose what linear solver to use for the stiff case?

  1. If the problem is size is fairly small (say N < 100), then using the dense solver is probably best; it is the simplest to use, and reasonably inexpensive for small N. For larger N, it is important to take advantage of sparsity (zero-nonzero) structure within the problem.
  2. If there is local (nearest-neighbor) coupling, or if the coupling is local after a suitable reordering of y, then use the banded linear solver. Local coupling means that the i-th component of the RHS or residual function depends only on components y_j for which |i-j| is small relative to N. (Note that the dense and band solvers are only applicable for the serial version of the solver.)
  3. For even larger problems, consider one of the Krylov iterative methods. These are hardest to use, because for best results they usually require preconditioning. However they offer the best opportunity to exploit the sparsity structure in the problem. The preconditioner is a matrix which, at least crudely, approximates the actual matrix in the linear system to be solved, and is typically built from an approximation of the relevant Jacobian matrix. Typically, that approximation uses only part of the true Jacobian, but as a result is much less expensive to solve. If the Jacobian can be approximated by a matrix that is banded (serial case) or block-diagonal with banded blocks (parallel case), SUNDIALS includes preconditioner modules for such cases. In each of the user guides, the section 'Linear solver specification functions' and the section on preconditioner modules contain more detailed comments on preconditioning.
  4. On the construction of preconditioners for problems arising from the spatial discretization of time-dependent partial differential equation systems, there is considerable discussion in the paper P. N. Brown and A. C. Hindmarsh, "Reduced Storage Matrix Methods in Stiff ODE Systems," J. Appl. Math. & Comp., 31 (1989), pp.40-91.

How do I handle a data-defined function within the RHS function?

Often the RHS or residual function depends on some function A(t) that is data-defined, i.e. defined only at a set of discrete set of times t. The solver must be able to obtain values of the user-supplied functions at arbitrary times t in the integration interval. So the user must fit the data with a reasonably smooth function A(t) that is defined continuously for all relevant t, and incorporate an evaluation of that fit function in the user function involved. This may be as simple as a piecewise linear fit, but a smoother fit (e.g. spline) would make the integration more efficient. If there is noise in the data, the fit should be a least-squares fit instead of a straight interpolation. The same advice applies if the user function has a data-defined function A(y) that involves one or more components of the dependent variable vector y. Of course, if more that one component is involved, the fit is more complicated.

How do I control unphysical negative values?

In many applications, some components in the true solution are always positive or non-negative, though at times very small. In the numerical solution, however, small negative (hence unphysical) values can then occur. In most cases, these values are harmless, and simply need to be controlled, not eliminated. The following pieces of advice are relevant. (See also the Note on this subject.)

  1. The way to control the size of unwanted negative computed values is with tighter absolute tolerances. Again this requires some knowledge of the noise level of these components, which may or may not be different for different components. Some experimentation may be needed.
  2. If output plots or tables are being generated, and it is important to avoid having negative numbers appear there (for the sake of avoiding a long explanation of them, if nothing else), then eliminate them, but only in the context of the output medium. Then the internal values carried by the solver are unaffected. Remember that a small negative value in y returned by CVODE, with magnitude comparable to abstol or less, is equivalent to zero as far as the computation is concerned.
  3. The user's right-hand side routine f should never change a negative value in the solution vector y to a non-negative value, as a "solution" to this problem. This can cause instability. If the f routine cannot tolerate a zero or negative value (e.g. because there is a square root or log of it), then the offending value should be changed to zero or a tiny positive number in a temporary variable (not in the input y vector) for the purposes of computing f(t,y).

How do I treat discontinuities in the RHS function?

If the jumps at the discontinuities are relatively small, simply keep them in the RHS function, and let the integrator respond to them (possibly taking smaller steps through each point of discontinuity). If the jumps are large, it is more efficient to stop at the point of discontinuity and restart the integrator with a readjusted ODE model. To stop when the location of the discontinuity is known, simply make that location a value of tout. To stop when the location of the discontinuity is determined by the solution, use the rootfinding feature. In either case, it is critical that the RHS function not incorporate the discontinuity, but rather have a smooth extension over the discontinuity, so that the step across it (and subsequent rootfinding, if used) can be done efficiently. Then use a switch within the RHS function that can be flipped between the stopping of the integration and the restart, so that the restarted problem uses the new values (which have jumped).

When is it advantageous to supply my own EwtFn function?

The main situation where this is a good idea is where the problem needs something "in between" the cases covered by CV_SS and CV_SV. Namely, suppose there are a few groups of variables (relative to the total number of variables) such that all the variables in each group require the same value of abstol, but these values are very different from one group to another. Then a user EwtFn function can keep an array of those values and construct the ewt vector without any additional storage. Also, in rare cases, one may want to use this option to apply different values of reltol to different variables (or groups of variables).

How do switch on/off forward sensitivity computations in CVODES?

If you want to turn on and off forward sensitivity calculations during several successive integrations (such as if you were using CVODES within a dynamically-constrained optimization loop, when sometimes you want to only integrate the states and sometimes you also need sensitivities computed), it is most efficient to use CVodeSensToggleOff.

Sensitivity calculations are enabled by the following functions: CVodeSensMalloc and CVodeSensReInit and are disabled by CVodeSensFree (after calling this one, they can be re-enabled only by calling CVodeSensMalloc) and CVodeSensToggleOff.

What is the role of plist in CVODES?

The argument plist to CVodeSetSensParams is used to specify the problem parameters with respect to which solution sensitivities are to be computed.

plist is used only if the sensitivity right-hand sides are evaluated using the internal difference-quotient approximation function. In that case, plist should be declared as an array of Ns integers and should contain the indices in the array of problem parameters p with respect to which sensitivities are desired. For example, if you want to compute sensitivities with respect to the first and third parameters in the p array, p[0] and p[2], you need to set

plist[0] = 0
plist[1] = 2

If plist is not provided, CVODES will compute sensitivities with respect to the first Ns parameters in the array p (i.e. it will use plist[i]=i, i=0,1,...Ns).
If the user provides a function to evaluate the right-hand sides of the sensitivity equations or if the default values are desired, a NULL pointer can be passed to CVodeSetSensParams.

What is the role of pbar in CVODES?

The argument pbar to CVodeSetSensParams is used to specify scaling factors for the problem parameters.

pbar is used only if
♦ the internal difference-quotient functions are used for the evaluation of the sensitivity right-hand sides, in which case pbar is used in computing an appropriate perturbation for the finite-difference approximation
or
♦ the tolerances for the sensitivity variables are estimated automatically by CVODES from those specified for the state variables.

If provided, pbar should be declared as an array of Ns realtypes and should contain non-zero scaling factors for the Ns parameters with respect to which sensitivities are to be computed. For non-zero problem parameters, a good choice is

pbar[i] = p[plist[i]]

If pbar is not provided, CVODES will use pbar[i]=1.0, i=0,1,...Ns-1.

If the user provides a function to evaluate the right-hand sides of the sensitivity equations and also specifies tolerances for the sensitivity variables (through CVodeSetSensTolerances) or if the default values are desired, a NULL pointer can be passed to CVodeSetSensParams.

What is pure quadrature integration?

Suppose your ODE is y'=F(t,y) and you integrate it from 0 to T and that you are also interested in computing an integral of the form

G = int_0^T g(t,y(t)) dt

for some function g. The most efficient way of computing z is by appending one additional differential equation to your ODE system:

z' = g(t,y)

with initial condition z(0)=0, in which case

G = z(T)

This additional equation is "a pure quadrature equation" and its main characteristic is that the new differential variable z does not appear in the right hand side of the extended ODE system. If CVODES is notified of such "pure quadrature equations", it can take advantage of this property and do less work than if it didn't know about them (these variables need not be considered in the nonlinear system solution).

The main reason for the special treatment of "pure quadrature equations" in CVODES is that such integrals (very often a large number of them) need to be computed for adjoint sensitivity.

IDA

How do I choose tolerances in IDA?

See How do I choose tolerances? in the CVODE/CVODES section.

How do I choose what linear solver to use?

See How do I choose what linear solver to use for the stiff case? in the CVODE/CVODES section.

How do I handle a data-defined function within the residual function?

See How do I handle a data-defined function within the RHS function? in the CVODE/CVODES section.

How do I control unphysical negative values in IDA?

See How do I control unphysical negative values? in the CVODE/CVODES section. (See also the Note on this subject.)
In addition, IDA provides the option of enforcing positivity or non-negativity on components. But these constraint options should only be exercised if the use of absolute tolerances to control the computed values is unsuccessful, because they involve some extra overhead cost.

How do I treat discontinuities in the residual function in IDA?

See How do I treat discontinuities in the RHS function? in the CVODE/CVODES section.

When is it advantageous to supply my own EwtFn function in IDA?

See How do I treat discontinuities in the RHS function? in the CVODE/CVODES section.

KINSOL

How do I reinitialize KINSOL within a C/C++ program?

Although KINSOL does not provide a reinitialization function, it is possible to reinitialize the solver (meaning reuse a KINSOL object), but only if the problem size remains unchanged. To reinitialize KINSOL, begin by making any necessary changes to the problem definition by calling the appropriate KINSet* functions (e.g., KINSetSysFunc). Next, if you would like to use a different linear solver, call the appropriate function (e.g., KINDense) followed by any calls to the corresponding KIN*Set* functions. Then you can call the KINSol function to solve the updated nonlinear algebraic system.

Why is the system function being evaluated at points that violate the constraints?

If you have not supplied a function to compute either J(u) (of type KINDenseJacFn or KINBandJacFn) or J(u)v (of type KINSpilsJacTimesVecFn), then the internal function may be the culprit. The default function used to compute a difference quotient approximation to the Jacobian (direct methods) or Jacobian matrix-vector product (Kylov methods) evaluates the user-supplied system function at a slightly perturbed point, but does not check if that point violates the constraints.

Miscellaneous

How do I determine which version of SUNDIALS I have?

If you still have access to the distribution files, then the SUNDIALS release number is indicated in the header of sundials/README and the corresponding solver versions can be determined by reading the appropriate row of the "Release History" table.

The specific version number of each solver is contained in corresponding README files: sundials/src/<solver>/README.

How do I apply a patch under UNIX or Linux?

To apply patch foobar.patch and rebuild SUNDIALS complete the following steps:

% cp foobar.patch <sundials_source_dir>
% cd <sundials_source_dir>
% patch -p1 < foobar.patch
% cd <sundials_source_dir>
% make
% make install