A brief introduction to CVODE¶
CVODE is part of a software family called sundials for SUite of Nonlinear and DIfferential / ALgebraic equation Solver [LLNL2005].
The version used with PelePhysics for the tests presented in this document is v4.0.2. This does not mean this is the recommended version to link with PelePhysics. Always refer to Install CVODE and SuiteSparse to know which versions of CVODE/SuiteSparse to use !!
In the following, details pertaining to the methods implemented in PelePhysics are provided. The interested user is referred to the very exhaustive CVODE Userguide for more information.
Most of this section is adapted from the v4.1.0 Cvode Documentation.
Numerical methods overview¶
The methods implemented in CVODE are variableorder, variablestep multistep methods, based on formulas of the form
Here the \(y^n\) are computed approximations to \(y(t_n)\), and \(h_n = t_nt_{n1}\) is the step size. For stiff problems, CVODE includes the Backward Differentiation Formulas (BDF) in socalled fixedleading coefficient (FLC) form, given by \(K_1=q\) and \(K_2= 0\), with order \(q\) varying between 1 and 5. The coefficients are uniquely determined by the method type, its order, the recent history of the step sizes, and the normalization \(\alpha_{n,0}=1\) [BYRNE1975], [JAC1980].
A nonlinear system must be solved (approximately) at each integration step. This nonlinear system can be formulated as a rootfinding problem
where \(a_n = \sum_{i>0} (\alpha_{n,i} y^{ni} + h_n\beta_{n,i} \dot{y}^{ni})\). CVODE provides several nonlinear solver choices. By default, CVODE solves this problem with a Newton iteration, which requires the solution of linear systems
in which
Linear Algebra¶
To find the solution of the linear system (1); CVODE provides several linear solver choices. The linear solver modules distributed with Sundials are organized in two families, a direct family comprising direct linear solvers for dense, banded, or sparse matrices, and a spils family comprising scaled preconditioned iterative (Krylov) linear solvers. The solvers offered through these modules that are of interest to us are:
 a dense direct solver
 a sparse direct solver interface using the KLU sparse solver library
 SPGMR, a scaled possibly preconditioned GMRES (Generalized Minimal Residual method) solver [BROWN1990]
When using a dense direct solver, the user has the option to specify an Analytical Jacobian. If none is provided, a difference quotients is performed. When a sparse direct solver is employed however, the user must specify an analytical Jacobian. All of these options have been enabled in PelePhysics.
For large stiff systems, where direct methods are often not feasible, the combination of a BDF integrator and a preconditioned Krylov method yields a powerful tool. In this case, the linear solve is matrixfree, and the default Newton iteration is an Inexact Newton iteration, in which \(M\) is applied with matrixvector products \(Jv\) obtained by either difference quotients or a usersupplied routine. In PelePhysics, it is possible to use either a nonpreconditioned or a preconditioned GMRES solver. In the latter case, the preconditioner can be either specified in a dense or sparse format (if the KLU library is linked to CVODE), and it is provided in the form of a Jacobian approximation, based on the work of [McNenly2015].
Error control, stepsizing, order determination¶
In the process of controlling errors at various levels, CVODE uses a weighted rootmeansquare norm, denoted \( \bullet _{WRMS}\), for all errorlike quantities. The multiplicative weights used are based on the current solution and on the relative and absolute tolerances input by the user, namely
Because \(1/W_i\) represents a tolerance in the component \(y_i\), a vector whose norm is 1 is regarded as small. In PelePhysics, both these tolerances are fixed to a value of \(1.0e10\).
A critical part of CVODE  making it an ODE solver rather than just an ODE method, is its control of the local error. At every step, the local error is estimated and required to satisfy tolerance conditions, and the step is redone with reduced step size whenever that error test fails. Note that in PelePhysics, the first time step is always forced to \(1.0e9\).
In addition to adjusting the step size to meet the local error test, CVODE periodically adjusts the order, with the goal of maximizing the step size. The integration starts out at order 1 and varies the order dynamically after that. The basic idea is to pick the order \(q\) for which a polynomial of order \(q\) best fits the discrete data involved in the multistep method. In PelePhysics, the maximum order is limited to 2 for stability reasons.
The various algorithmic features of CVODE are inherited from VODE and VODPK, and are documented in [VODE1989] and [BROWN1990]. They are also summarized in the CVODE Userguide as well as in [LLNL2005].
[LLNL2005]  (1, 2)

[BYRNE1975] 

[JAC1980] 

[BROWN1990]  (1, 2)

[McNenly2015] 

[VODE1989] 
