6.4. Using IDA for IVP Solution

This chapter is concerned with the use of IDA for the integration of DAEs.

The following sections treat the header files and the layout of the user’s main program, and provide descriptions of the IDA user-callable functions and user-supplied functions. The sample programs described in the companion document [74] may also be helpful. Those codes may be used as templates (with the removal of some lines used in testing) and are included in the IDA package.

IDA uses various constants for both input and output. These are defined as needed in this chapter, but for convenience are also listed separately in §6.5.

The user should be aware that not all SUNLinearSolver and SUNMatrix objects are compatible with all N_Vector implementations. Details on compatibility are given in the documentation for each SUNMatrix (Chapter §10) and SUNLinearSolver (Chapter §11) implementation. For example, NVECTOR_PARALLEL is not compatible with the dense, banded, or sparse SUNMatrix types, or with the corresponding dense, banded, or sparse SUNLinearSolver objects. Please check Chapters §10 and §11 to verify compatibility between these objects. In addition to that documentation, we note that the IDABBDPRE preconditioner can only be used with NVECTOR_PARALLEL. It is not recommended to use a threaded vector object with SuperLU_MT unless it is the NVECTOR_OPENMP module, and SuperLU_MT is also compiled with OpenMP.

6.4.1. Access to library and header files

At this point, it is assumed that the installation of IDA, following the procedure described in §2.1, has been completed successfully. In the proceeding text, the directories libdir and incdir are the installation library and include directories, respectively. For a default installation, these are instdir/lib and instdir/include, respectively, where instdir is the directory where SUNDIALS was installed.

Regardless of where the user’s application program resides, its associated compilation and load commands must make reference to the appropriate locations for the library and header files required by IDA. IDA symbols are found in libdir/libsundials_ida.lib. Thus, in addition to linking to libdir/libsundials_core.lib, IDA users need to link to the IDA library. Symbols for additional SUNDIALS modules, vectors and algebraic solvers, are found in

<libdir>/libsundials_nvec*.lib
<libdir>/libsundials_sunmat*.lib
<libdir>/libsundials_sunlinsol*.lib
<libdir>/libsundials_sunnonlinsol*.lib
<libdir>/libsundials_sunmem*.lib

The file extension .lib is typically .so for shared libraries and .a for static libraries.

The relevant header files for IDA are located in the subdirectories incdir/include/ida. To use IDA the application needs to include the header file for IDA in addition to the SUNDIALS core header file:

#include <sundials/sundials_core.h> // Provides core SUNDIALS types
#include <ida/ida.h>                // IDA provides methods for DAEs

The calling program must also include an N_Vector implementation header file, of the form nvector/nvector_*.h. See §9 for the appropriate name.

If using a non-default nonlinear solver object, or when interacting with a SUNNonlinearSolver object directly, the calling program must also include a SUNNonlinearSolver implementation header file, of the form sunnonlinsol/sunnonlinsol_*.h where * is the name of the nonlinear solver (see Chapter §12 for more information).

If using a nonlinear solver that requires the solution of a linear system of the form (6.3) (e.g., the default Newton iteration), the calling program must also include a SUNLinearSolver implementation header file, of the from sunlinsol/sunlinsol_*.h where * is the name of the linear solver (see Chapter §11 for more information).

If the linear solver is matrix-based, the linear solver header will also include a header file of the from sunmatrix/sunmatrix_*.h where * is the name of the matrix implementation compatible with the linear solver (see Chapter §10 for more information).

Other headers may be needed, according to the choice of preconditioner, etc. For example, in the example idaFoodWeb_kry_p (see [74]), preconditioning is done with a block-diagonal matrix. For this, even though the SUNLINSOL_SPGMR linear solver is used, the header sundials/sundials_dense.h is included for access to the underlying generic dense matrix arithmetic routines.

6.4.2. 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 a DAE IVP. Most of the steps are independent of the N_Vector, SUNMatrix, SUNLinearSolver, and SUNNonlinearSolver implementations used. For the steps that are not, refer to Chapters §9, §10, §11, and §12 for the specific name of the function to be called or macro to be referenced.

Initialize parallel or multi-threaded environment (if appropriate)

For example, call MPI_Init to initialize MPI if used.
Create the SUNDIALS context object

Call SUNContext_Create() to allocate the SUNContext object.
Create the vector of initial values

Construct an N_Vector of initial values using the appropriate functions defined by the particular N_Vector implementation (see §9 for details).

For native SUNDIALS vector implementations, use a call of the form y0 = N_VMake_***(..., ydata) if the array containing the initial values of \(y\) already exists. Otherwise, create a new vector by making a call of the form N_VNew_***(...), and then set its elements by accessing the underlying data with a call of the form ydata = N_VGetArrayPointer(y0). Here, *** is the name of the vector implementation.

For hypre, PETSc, and Trilinos vector wrappers, first create and initialize the underlying vector, and then create an N_Vector wrapper with a call of the form y0 = N_VMake_***(yvec), where yvec is a hypre, PETSc, or Trilinos vector. Note that calls like N_VNew_***(...) and N_VGetArrayPointer(...) are not available for these vector wrappers.

Set the vector yp0 of initial conditions for \(\dot{y}\) similarly.
Create matrix object (if appropriate)

If a linear solver is required (e.g., when using the default Newton solver) and the linear solver will be a matrix-based linear solver, then a template Jacobian matrix must be created by calling the appropriate constructor defined by the particular SUNMatrix implementation.

For the native SUNDIALS SUNMatrix implementations, the matrix object may be created using a call of the form SUN***Matrix(...) where *** is the name of the matrix (see §10 for details).
Create linear solver object (if appropriate)

If a linear solver is required (e.g., when using the default Newton solver), then the desired linear solver object must be created by calling the appropriate constructor defined by the particular SUNLinearSolver implementation.

For any of the native SUNDIALS SUNLinearSolver implementations, the linear solver object may be created using a call of the form SUNLinearSolver LS = SUNLinSol_***(...); where *** is the name of the linear solver (see §11 for details).
Create nonlinear solver object (if appropriate)

If using a non-default nonlinear solver, then the desired nonlinear solver object must be created by calling the appropriate constructor defined by the particular SUNNonlinearSolver implementation.

For any of the native SUNDIALS SUNNonLinearSolver implementations, the nonlinear solver object may be created using a call of the form SUNNonlinearSolver NLS = SUNNonlinSol_***(...); where *** is the name of the nonlinear solver (see §12 for details).
Create IDA object

Call IDACreate() to create the IDA solver object.
Initialize IDA solver

Call IDAInit() to provide the initial condition vectors created above, set the DAE residual function, and initialize IDA.
Specify integration tolerances

Call one of the following functions to set the integration tolerances:
- IDASStolerances() to specify scalar relative and absolute tolerances.
- IDASVtolerances() to specify a scalar relative tolerance and a vector of absolute tolerances.
- IDAWFtolerances() to specify a function which sets directly the weights used in evaluating WRMS vector norms.
See §6.4.3.3 for general advice on selecting tolerances and §6.4.3.4 for advice on controlling unphysical values.
Attach the linear solver (if appropriate)

If a linear solver was created above, initialize the IDALS linear solver interface by attaching the linear solver object (and matrix object, if applicable) with IDASetLinearSolver().
Set linear solver optional inputs (if appropriate)

See Table 6.2 for IDALS optional inputs and Chapter §11 for linear solver specific optional inputs.
Attach nonlinear solver module (if appropriate)

If a nonlinear solver was created above, initialize the IDANLS nonlinear solver interface by attaching the nonlinear solver object with IDASetNonlinearSolver().
Set nonlinear solver optional inputs (if appropriate)

See Table 6.3 for IDANLS optional inputs and Chapter §12 for nonlinear solver specific optional inputs. Note, solver specific optional inputs must be called after IDASetNonlinearSolver(), otherwise the optional inputs will be overridden by IDA defaults.
Specify rootfinding problem (optional)

Call IDARootInit() to initialize a rootfinding problem to be solved during the integration of the ODE system. See Table 6.6 for relevant optional input calls.
Set optional inputs

Call IDASet*** functions to change any optional inputs that control the behavior of IDA from their default values. See §6.4.3.10 for details.
Correct initial values (optional)

Call IDACalcIC() to correct the initial values y0 and yp0 passed to IDAInit(). See Table 6.4 for relevant optional input calls.
Advance solution in time

For each point at which output is desired, call ier = IDASolve(ida_mem, tout, &tret, yret, ypret, itask). Here itask specifies the return mode. The vector yret (which can be the same as the vector y0 above) will contain \(y(t)\), while the vector ypret (which can be the same as the vector yp0 above) will contain \(\dot{y}(t)\).

See IDASolve() for details.
Get optional outputs

Call IDAGet*** functions to obtain optional output. See §6.4.3.12 for details.
Destroy objects

Upon completion of the integration call the following functions, as necessary, to destroy any objects created above:
- Call N_VDestroy() to free vector objects.
- Call SUNMatDestroy() to free matrix objects.
- Call SUNLinSolFree() to free linear solvers objects.
- Call SUNNonlinSolFree() to free nonlinear solvers objects.
- Call IDAFree() to free the memory allocated by IDA.
- Call SUNContext_Free() to free the SUNDIALS context.
Finalize MPI, if used

Call MPI_Finalize to terminate MPI.

6.4.3. User-callable functions

This section describes the IDA functions that are called by the user to setup and then solve an IVP. Some of these are required. However, starting with §6.4.3.10, the functions listed involve optional inputs/outputs or restarting, and those paragraphs may be skipped for a casual use of IDA. In any case, refer to §6.4.2 for the correct order of these calls.

On an error, each user-callable function returns a negative value and sends an error message to the error handler routine, which prints the message on stderr by default. However, the user can set a file as error output or can provide his own error handler function (see §6.4.3.10.1).

6.4.3.1. IDA initialization and deallocation functions

void *IDACreate(SUNContext sunctx)

The function IDACreate instantiates an IDA solver object.

Arguments:

sunctx – the SUNContext object (see §2.4)

Return value:

void* pointer the IDA solver object.

int IDAInit(void *ida_mem, IDAResFn res, sunrealtype t0, N_Vector y0, N_Vector yp0)

The function IDAInit provides required problem and solution specifications, allocates internal memory, and initializes IDA.

Arguments:

ida_mem – pointer to the IDA solver object.
res – is the function which computes the residual function \(F(t, y, \dot{y})\) for the DAE. For full details see IDAResFn.
t0 – is the initial value of \(t\).
y0 – is the initial value of \(y\).
yp0 – is the initial value of \(\dot{y}\).

Return value:

IDA_SUCCESS – The call was successful.
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_MEM_FAIL – A memory allocation request has failed.
IDA_ILL_INPUT – An input argument to IDAInit() has an illegal value.

Notes:

If an error occurred, IDAInit() also sends an error message to the error handler function.

void IDAFree(void **ida_mem)

The function IDAFree frees the pointer allocated by a previous call to IDACreate().

Arguments:

ida_mem – pointer to the IDA solver object.

Return value:

void

6.4.3.2. IDA tolerance specification functions

One of the following three functions must be called to specify the integration tolerances (or directly specify the weights used in evaluating WRMS vector norms). Note that this call must be made after the call to IDAInit().

int IDASStolerances(void *ida_mem, sunrealtype reltol, sunrealtype abstol)

The function IDASStolerances specifies scalar relative and absolute tolerances.

Arguments:

ida_mem – pointer to the IDA solver object.
reltol – is the scalar relative error tolerance.
abstol – is the scalar absolute error tolerance.

Return value:

IDA_SUCCESS – The call was successful
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_NO_MALLOC – The allocation function IDAInit() has not been called.
IDA_ILL_INPUT – One of the input tolerances was negative.

int IDASVtolerances(void *ida_mem, sunrealtype reltol, N_Vector abstol)

The function IDASVtolerances specifies scalar relative tolerance and vector absolute tolerances.

Arguments:

ida_mem – pointer to the IDA solver object.
reltol – is the scalar relative error tolerance.
abstol – is the vector of absolute error tolerances.

Return value:

IDA_SUCCESS – The call was successful
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_NO_MALLOC – The allocation function IDAInit() has not been called.
IDA_ILL_INPUT – The relative error tolerance was negative or the absolute tolerance vector had a negative component.

Notes:

This choice of tolerances is important when the absolute error tolerance needs to be different for each component of the state vector \(y\).

int IDAWFtolerances(void *ida_mem, IDAEwtFn efun)

The function IDAWFtolerances specifies a user-supplied function efun that sets the multiplicative error weights \(W_i\) for use in the weighted RMS norm, which are normally defined by (6.4).

Arguments:

ida_mem – pointer to the IDA solver object. IDACreate()
efun – is the function which defines the ewt vector. For full details see IDAEwtFn.

Return value:

IDA_SUCCESS – The call was successful
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_NO_MALLOC – The allocation function IDAInit() has not been called.

6.4.3.3. General advice on choice of tolerances

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

The scalar relative tolerance reltol is to be set to control relative errors. So reltol of \(10^{-4}\) means that errors are controlled to .01%. We do not recommend using reltol larger than \(10^{-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 \(10^{-15}\)).
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 a 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 idaRoberts_dns in the IDA package, and the discussion of it in the IDA Examples document [74]. In that problem, the three components vary betwen 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.
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 to is a reltol of \(10^{-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.

6.4.3.4. Advice on controlling 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.

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.
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 yret returned by IDA, with magnitude comparable to abstol or less, is equivalent to zero as far as the computation is concerned.
The user’s residual function res should never change a negative value in the solution vector yy to a non-negative value, as a “solution” to this problem. This can cause instability. If the res 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 yy vector) for the purposes of computing \(F(t,y,\dot{y})\).
IDA provides the option of enforcing positivity or non-negativity on components. Also, such constraints can be enforced by use of the recoverable error return feature in the user-supplied residual function. However, because these options involve some extra overhead cost, they should only be exercised if the use of absolute tolerances to control the computed values is unsuccessful.

6.4.3.5. Linear solver interface functions

As previously explained, if the nonlinear solver requires the solution of linear systems of the form (6.5), e.g., the default Newton solver, then the solution of these linear systems is handled with the IDALS linear solver interface. This interface supports all valid SUNLinearSolver objects. Here, a matrix-based SUNLinearSolver utilizes SUNMatrix objects to store the Jacobian matrix \(J = \dfrac{\partial{F}}{\partial{y}} + \alpha \dfrac{\partial{F}}{\partial{\dot{y}}}\) and factorizations used throughout the solution process. Conversely, matrix-free SUNLinearSolver object instead use iterative methods to solve the linear systems of equations, and only require the action of the Jacobian on a vector, \(Jv\).

With most iterative linear solvers, preconditioning can be done on the left only, on the right only, on both the left and the right, or not at all. The exceptions to this rule are SPFGMR that supports right preconditioning only and PCG that performs symmetric preconditioning. However, in IDA only left preconditioning is supported. For the specification of a preconditioner, see the iterative linear solver sections in §6.4.3.10 and §6.4.4. A preconditioner matrix \(P\) must approximate the Jacobian \(J\), at least crudely.

To attach a generic linear solver to IDA, after the call to IDACreate() but before any calls to IDASolve(), the user’s program must create the appropriate SUNLinearSolver object and call the function IDASetLinearSolver(). To create the SUNLinearSolver object, the user may call one of the SUNDIALS-packaged SUNLinearSolver constructors via a call of the form

SUNLinearSolver LS = SUNLinSol_*(...);

Alternately, a user-supplied SUNLinearSolver object may be created and used instead. The use of each of the generic linear solvers involves certain constants, functions and possibly some macros, that are likely to be needed in the user code. These are available in the corresponding header file associated with the specific SUNMatrix or SUNLinearSolver object in question, as described in Chapters §10 and §11.

Once this solver object has been constructed, the user should attach it to IDA via a call to IDASetLinearSolver(). The first argument passed to this function is the IDA memory pointer returned by IDACreate(); the second argument is the desired SUNLinearSolver object to use for solving systems. The third argument is an optional SUNMatrix object to accompany matrix-based SUNLinearSolver inputs (for matrix-free linear solvers, the third argument should be NULL). A call to this function initializes the IDALS linear solver interface, linking it to the main IDA integrator, and allows the user to specify additional parameters and routines pertinent to their choice of linear solver.

int IDASetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix J)

The function IDASetLinearSolver attaches a SUNLinearSolver object LS and corresponding template Jacobian SUNMatrix object J (if applicable) to IDA, initializing the IDALS linear solver interface.

Arguments:

ida_mem – pointer to the IDA solver object.
LS – SUNLinearSolver object to use for solving linear systems of the form (6.5).
J – SUNMatrix object for used as a template for the Jacobian or NULL if not applicable.

Return value:

IDALS_SUCCESS – The IDALS initialization was successful.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_ILL_INPUT – The IDALS interface is not compatible with the LS or J input objects or is incompatible with the N_Vector object passed to IDAInit().
IDALS_SUNLS_FAIL – A call to the LS object failed.
IDALS_MEM_FAIL – A memory allocation request failed.

Notes:

If LS is a matrix-based linear solver, then the template Jacobian matrix J will be used in the solve process, so if additional storage is required within the SUNMatrix object (e.g., for factorization of a banded matrix), ensure that the input object is allocated with sufficient size (see the documentation of the particular SUNMatrix in Chapter §10 for further information).

Warning

The previous routines IDADlsSetLinearSolver and IDASpilsSetLinearSolver are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

6.4.3.6. Nonlinear solver interface function

By default IDA uses the SUNNonlinearSolver implementation of Newton’s method (see §12.7). To attach a different nonlinear solver in IDA, the user’s program must create a SUNNonlinearSolver object by calling the appropriate constructor routine. The user must then attach the SUNNonlinearSolver object to IDA by calling IDASetNonlinearSolver().

When changing the nonlinear solver in IDA, IDASetNonlinearSolver() must be called after IDAInit(). If any calls to IDASolve() have been made, then IDA will need to be reinitialized by calling IDAReInit() to ensure that the nonlinear solver is initialized correctly before any subsequent calls to IDASolve().

The first argument passed to IDASetNonlinearSolver() is the IDA memory pointer returned by IDACreate() and the second argument is the SUNNonlinearSolver object to use for solving the nonlinear system (6.3). A call to this function attaches the nonlinear solver to the main IDA integrator. We note that at present, the SUNNonlinearSolver object must be of type SUNNONLINEARSOLVER_ROOTFIND.

int IDASetNonlinearSolver(void *ida_mem, SUNNonlinearSolver NLS)

The function IDASetNonLinearSolver attaches a SUNNonlinearSolver object (NLS) to IDA.

Arguments:

ida_mem – pointer to the IDA solver object.
NLS – SUNNonlinearSolver object to use for solving nonlinear systems.

Return value:

IDA_SUCCESS – The nonlinear solver was successfully attached.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The SUNNonlinearSolver object is NULL , does not implement the required nonlinear solver operations, is not of the correct type, or the residual function, convergence test function, or maximum number of nonlinear iterations could not be set.

6.4.3.7. Initial condition calculation function

IDACalcIC() calculates corrected initial conditions for the DAE system for certain index-one problems including a class of systems of semi-implicit form (see §6.2.2 and [24]). It uses a Newton iteration combined with a linesearch algorithm. Calling IDACalcIC() is optional. It is only necessary when the initial conditions do not satisfy the given system. Thus if y0 and yp0 are known to satisfy \(F(t_0, y_0, \dot{y}_0) = 0\), then a call to IDACalcIC() is generally not necessary.

A call to the function IDACalcIC() must be preceded by successful calls to IDACreate() and IDAInit() (or IDAReInit()), and by a successful call to the linear system solver specification function. The call to IDACalcIC() should precede the call(s) to IDASolve() for the given problem.

int IDACalcIC(void *ida_mem, int icopt, sunrealtype tout1)

The function IDACalcIC corrects the initial values y0 and yp0 at time t0.

Arguments:

ida_mem – pointer to the IDA solver object.
icopt – is one of the following two options for the initial condition calculation.
- IDA_YA_YDP_INIT directs IDACalcIC() to compute the algebraic components of \(y\) and differential components of \(\dot{y}\), given the differential components of \(y\). This option requires that the N_Vector id was set through IDASetId(), specifying the differential and algebraic components.
- IDA_Y_INIT directs IDACalcIC() to compute all components of \(y\), given \(\dot{y}\). In this case, id is not required.
tout1 – is the first value of \(t\) at which a solution will be requested (from IDASolve()). This value is needed here only to determine the direction of integration and rough scale in the independent variable \(t\).

Return value:

IDA_SUCCESS – IDACalcIC() succeeded.
IDA_MEM_NULL – The argument ida_mem was NULL.
IDA_NO_MALLOC – The allocation function IDAInit() has not been called.
IDA_ILL_INPUT – One of the input arguments was illegal.
IDA_LSETUP_FAIL – The linear solver’s setup function failed in an unrecoverable manner.
IDA_LINIT_FAIL – The linear solver’s initialization function failed.
IDA_LSOLVE_FAIL – The linear solver’s solve function failed in an unrecoverable manner.
IDA_BAD_EWT – Some component of the error weight vector is zero (illegal), either for the input value of y0 or a corrected value.
IDA_FIRST_RES_FAIL – The user’s residual function returned a recoverable error flag on the first call, but IDACalcIC() was unable to recover.
IDA_RES_FAIL – The user’s residual function returned a nonrecoverable error flag.
IDA_NO_RECOVERY – The user’s residual function, or the linear solver’s setup or solve function had a recoverable error, but IDACalcIC() was unable to recover.
IDA_CONSTR_FAIL – IDACalcIC() was unable to find a solution satisfying the inequality constraints.
IDA_LINESEARCH_FAIL – The linesearch algorithm failed to find a solution with a step larger than steptol in weighted RMS norm, and within the allowed number of backtracks.
IDA_CONV_FAIL – IDACalcIC() failed to get convergence of the Newton iterations.

Notes:

IDACalcIC() will correct the values of \(y(t_0)\) and \(\dot{y}(t_0)\) which were specified in the previous call to IDAInit() or IDAReInit(). To obtain the corrected values, call IDAGetConsistentIC().

6.4.3.8. Rootfinding initialization function

While solving the IVP, IDA has the capability to find the roots of a set of user-defined functions. To activate the root finding algorithm, call the following function. This is normally called only once, prior to the first call to IDASolve(), but if the rootfinding problem is to be changed during the solution, IDARootInit() can also be called prior to a continuation call to IDASolve().

int IDARootInit(void *ida_mem, int nrtfn, IDARootFn g)

The function IDARootInit specifies that the roots of a set of functions \(g_i(t,y)\) are to be found while the IVP is being solved.

Arguments:

ida_mem – pointer to the IDA solver object
nrtfn – is the number of root functions
g – is the function which defines the nrtfn functions \(g_i(t,y,\dot{y})\) whose roots are sought. See IDARootFn for more details.

Return value:

IDA_SUCCESS – The call was successful
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_MEM_FAIL – A memory allocation failed.
IDA_ILL_INPUT – The function g is NULL, but nrtfn > 0.

Notes:

If a new IVP is to be solved with a call to IDAReInit(), where the new IVP has no rootfinding problem but the prior one did, then call IDARootInit() with nrtfn = 0.

6.4.3.9. IDA solver function

This is the central step in the solution process, the call to perform the integration of the DAE. The input arguments (itask) specifies one of two modes as to where IDA is to return a solution. These modes are modified if the user has set a stop time (with IDASetStopTime()) or requested rootfinding (with IDARootInit()).

int IDASolve(void *ida_mem, sunrealtype tout, sunrealtype *tret, N_Vector yret, N_Vector ypret, int itask)

The function IDASolve integrates the DAE over an interval in t.

Arguments:

ida_mem – pointer to the IDA solver object.
tout – the next time at which a computed solution is desired.
tret – the time reached by the solver output.
yret – the computed solution vector y.
ypret – the computed solution vector \(\dot{y}\).
itask – a flag indicating the job of the solver for the next user step
- IDA_NORMAL – the solver will take internal steps until it has reached or just passed the user specified tout parameter. The solver then interpolates in order to return approximate values of \(y(t_{out})\) and \(\dot{y}(t_{out})\).
- IDA_ONE_STEP – the solver will just take one internal step and return the solution at the point reached by that step.

Return value:

IDA_SUCCESS – The call was successful.
IDA_TSTOP_RETURN – IDASolve() succeeded by reaching the stop point specified through the optional input function IDASetStopTime().
IDA_ROOT_RETURN – IDASolve() succeeded and found one or more roots. In this case, tret is the location of the root. If nrtfn >1, call IDAGetRootInfo() to see which \(g_i\) were found to have a root.
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_ILL_INPUT – One of the inputs to IDASolve() was illegal, or some other input to the solver was either illegal or missing. The latter category includes the following situations:
- The tolerances have not been set.
- A component of the error weight vector became zero during internal time-stepping.
- The linear solver initialization function called by the user after calling IDACreate() failed to set the linear solver-specific lsolve field in ida_mem.
- A root of one of the root functions was found both at a point \(t\) and also very near \(t\).
In any case, the user should see the printed error message for details.
IDA_TOO_MUCH_WORK – The solver took mxstep internal steps but could not reach tout. The default value for mxstep is MXSTEP_DEFAULT = 500.
IDA_TOO_MUCH_ACC – The solver could not satisfy the accuracy demanded by the user for some internal step.
IDA_ERR_FAIL – Error test failures occurred too many times (MXNEF = 10) during one internal time step or occurred with \(|h| = h_{\text{min}}\).
IDA_CONV_FAIL – Convergence test failures occurred too many times (MXNCF = 10) during one internal time step or occurred with \(|h| = h_{\text{min}}\).
IDA_LINIT_FAIL – The linear solver’s initialization function failed.
IDA_LSETUP_FAIL – The linear solver’s setup function failed in an unrecoverable manner.
IDA_LSOLVE_FAIL – The linear solver’s solve function failed in an unrecoverable manner.
IDA_CONSTR_FAIL – The inequality constraints were violated and the solver was unable to recover.
IDA_REP_RES_ERR – The user’s residual function repeatedly returned a recoverable error flag, but the solver was unable to recover.
IDA_RES_FAIL – The user’s residual function returned a nonrecoverable error flag.
IDA_RTFUNC_FAIL – The rootfinding function failed.

Notes:

The vectors yret and ypret can occupy the same space as the initial condition vectors y0 and yp0, respectively, that were passed to IDAInit().

In the IDA_ONE_STEP mode, tout is used on the first call only, and only to get the direction and rough scale of the independent variable.

If a stop time is enabled (through a call to IDASetStopTime()), then IDASolve() returns the solution at tstop. Once the integrator returns at a stop time, any future testing for tstop is disabled (and can be reenabled only though a new call to IDASetStopTime()).

All failure return values are negative and therefore a test flag < 0 will trap all IDASolve() failures.

On any error return in which one or more internal steps were taken by IDASolve(), the returned values of tret, yret, and ypret correspond to the farthest point reached in the integration. On all other error returns, these values are left unchanged from the previous IDASolve() return.

6.4.3.10. Optional input functions

There are numerous optional input parameters that control the behavior of the IDA solver. IDA provides functions that can be used to change these optional input parameters from their default values. The main inputs are divided in the following categories:

Table 6.1 list the main IDA optional input functions,
Table 6.2 lists the IDALS linear solver interface optional input functions,
Table 6.3 lists the IDANLS nonlinear solver interface optional input functions,
Table 6.4 lists the initial condition calculation optional input functions,
Table 6.5 lists the IDA step size adaptivity optional input functions, and
Table 6.6 lists the rootfinding optional input functions.

These optional inputs are described in detail in the remainder of this section. For the most casual use of IDA, the reader can skip to §6.4.4.

We note that, on an error return, all of the optional input functions also send an error message to the error handler function. All error return values are negative, so the test flag < 0 will catch all errors.

The optional input calls can, unless otherwise noted, be executed in any order. Finally, a call to an IDASet*** function can, unless otherwise noted, be made at any time from the user’s calling program and, if successful, takes effect immediately.

6.4.3.10.1. Main solver optional input functions

Table 6.1 Optional inputs for IDA
Optional input	Function name	Default
User data	`IDASetUserData()`	NULL
Maximum order for BDF method	`IDASetMaxOrd()`	5
Maximum no. of internal steps before \(t_{{\scriptsize out}}\)	`IDASetMaxNumSteps()`	500
Initial step size	`IDASetInitStep()`	estimated
Minimum absolute step size \(h_{\text{min}}\)	`IDASetMinStep()`	0
Maximum absolute step size \(h_{\text{max}}\)	`IDASetMaxStep()`	\(\infty\)
Value of \(t_{stop}\)	`IDASetStopTime()`	undefined
Disable the stop time	`IDAClearStopTime()`	N/A
Maximum no. of error test failures	`IDASetMaxErrTestFails()`	10
Suppress alg. vars. from error test	`IDASetSuppressAlg()`	`SUNFALSE`
Variable types (differential/algebraic)	`IDASetId()`	NULL
Inequality constraints on solution	`IDASetConstraints()`	NULL

int IDASetUserData(void *ida_mem, void *user_data)

The function IDASetUserData attaches a user-defined data pointer to the main IDA solver object.

Arguments:

ida_mem – pointer to the IDA solver object.
user_data – pointer to the user data.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

If specified, the pointer to user_data is passed to all user-supplied functions that have it as an argument. Otherwise, a NULL pointer is passed.

Warning

If user_data is needed in user linear solver or preconditioner functions, the call to IDASetUserData() must be made before the call to specify the linear solver.

int IDASetMaxOrd(void *ida_mem, int maxord)

The function IDASetMaxOrd specifies the maximum order of the linear multistep method.

Arguments:

ida_mem – pointer to the IDA solver object.
maxord – value of the maximum method order. This must be positive.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The input value maxord is \(\leq\) 0 , or larger than the max order value when IDAInit() was called.

Notes:

The default value is 5. If the input value exceeds 5, the value 5 will be used. If called before IDAInit(), maxord limits the memory requirements for the internal IDA memory block and its value cannot be increased past the value set when IDAInit() was called.

int IDASetMaxNumSteps(void *ida_mem, long int mxsteps)

The function IDASetMaxNumSteps specifies the maximum number of steps to be taken by the solver in its attempt to reach the next output time.

Arguments:

ida_mem – pointer to the IDA solver object.
mxsteps – maximum allowed number of steps.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

Passing mxsteps = 0 results in IDA using the default value (500). Passing mxsteps < 0 disables the test (not recommended).

int IDASetInitStep(void *ida_mem, sunrealtype hin)

The function IDASetInitStep specifies the initial step size.

Arguments:

ida_mem – pointer to the IDA solver object.
hin – value of the initial step size to be attempted. Pass 0.0 to have IDA use the default value.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

By default, IDA estimates the initial step as the solution of \(\|h \dot{y} \|_{{\scriptsize WRMS}} = 1/2\), with an added restriction that \(|h| \leq .001|t_{\text{out}} - t_0|\).

int IDASetMinStep(void *ida_mem, sunrealtype hmin)

The function IDASetMinStep specifies the minimum absolute value of the step size.

Pass hmin = 0 to obtain the default value of 0.

Arguments:

ida_mem – pointer to the IDA solver object.
hmin – minimum absolute value of the step size.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – hmin is negative.

New in version 6.2.0.

int IDASetMaxStep(void *ida_mem, sunrealtype hmax)

The function IDASetMaxStep specifies the maximum absolute value of the step size.

Arguments:

ida_mem – pointer to the IDA solver object.
hmax – maximum absolute value of the step size.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – Either hmax is not positive or it is smaller than the minimum allowable step.

Notes:

Pass hmax = 0 to obtain the default value \(\infty\).

int IDASetStopTime(void *ida_mem, sunrealtype tstop)

The function IDASetStopTime specifies the value of the independent variable \(t\) past which the solution is not to proceed.

Arguments:

ida_mem – pointer to the IDA solver object.
tstop – value of the independent variable past which the solution should not proceed.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The value of tstop is not beyond the current \(t\) value, \(t_n\).

Notes:

The default, if this routine is not called, is that no stop time is imposed. Once the integrator returns at a stop time, any future testing for tstop is disabled (and can be reenabled only though a new call to IDASetStopTime()).

A stop time not reached before a call to IDAReInit() will remain active but can be disabled by calling IDAClearStopTime().

int IDAClearStopTime(void *ida_mem)

Disables the stop time set with IDASetStopTime().

Arguments:

ida_mem – pointer to the IDA memory block.

Return value:

IDA_SUCCESS if successful
IDA_MEM_NULL if the IDA memory is NULL

Notes:

The stop time can be reenabled though a new call to IDASetStopTime().

New in version 6.5.1.

int IDASetMaxErrTestFails(void *ida_mem, int maxnef)

The function IDASetMaxErrTestFails specifies the maximum number of error test failures in attempting one step.

Arguments:

ida_mem – pointer to the IDA solver object.
maxnef – maximum number of error test failures allowed on one step (>0).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The default value is 10.

int IDASetSuppressAlg(void *ida_mem, sunbooleantype suppressalg)

The function IDASetSuppressAlg indicates whether or not to suppress algebraic variables in the local error test.

Arguments:

ida_mem – pointer to the IDA solver object.
suppressalg – indicates whether to suppress (SUNTRUE) or include (SUNFALSE) the algebraic variables in the local error test.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The default value is SUNFALSE. If suppressalg = SUNTRUE is selected, then the id vector must be set (through IDASetId()) to specify the algebraic components. In general, the use of this option (with suppressalg = SUNTRUE) is discouraged when solving DAE systems of index 1, whereas it is generally encouraged for systems of index 2 or more. See pp. 146-147 of [19] for more on this issue.

int IDASetId(void *ida_mem, N_Vector id)

The function IDASetId specifies algebraic/differential components in the \(y\) vector.

Arguments:

ida_mem – pointer to the IDA solver object.
id – a vector of values identifying the components of \(y\) as differential or algebraic variables. A value of 1.0 indicates a differential variable, while 0.0 indicates an algebraic variable.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The vector id is required if the algebraic variables are to be suppressed from the local error test (see IDASetSuppressAlg()) or if IDACalcIC() is to be called with icopt = IDA_YA_YDP_INIT.

int IDASetConstraints(void *ida_mem, N_Vector constraints)

The function IDASetConstraints specifies a vector defining inequality constraints for each component of the solution vector \(y\).

Arguments:

ida_mem – pointer to the IDA solver object.
constraints – vector of constraint flags.
- If constraints[i] = 0, no constraint is imposed on \(y_i\).
- If constraints[i] = 1, \(y_i\) will be constrained to be \(y_i \ge 0.0\).
- If constraints[i] = -1, \(y_i\) will be constrained to be \(y_i \le 0.0\).
- If constraints[i] = 2, \(y_i\) will be constrained to be \(y_i > 0.0\).
- If constraints[i] = -2, \(y_i\) will be constrained to be \(y_i < 0.0\).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The constraints vector contains illegal values.

Notes:

The presence of a non-NULL constraints vector that is not 0.0 in all components will cause constraint checking to be performed. However, a call with 0.0 in all components of constraints vector will result in an illegal input return. A NULL input will disable constraint checking.

6.4.3.10.2. Linear solver interface optional input functions

Table 6.2 Optional inputs for the IDALS linear solver interface
Optional input	Function name	Default
Jacobian function	`IDASetJacFn()`	DQ
Set parameter determining if a \(c_j\) change requires a linear solver setup call	`IDASetDeltaCjLSetup()`	0.25
Enable or disable linear solution scaling	`IDASetLinearSolutionScaling()`	on
Jacobian-times-vector function	`IDASetJacTimes()`	NULL, DQ
Preconditioner functions	`IDASetPreconditioner()`	NULL, NULL
Ratio between linear and nonlinear tolerances	`IDASetEpsLin()`	0.05
Increment factor used in DQ \(Jv\) approx.	`IDASetIncrementFactor()`	1.0
Jacobian-times-vector DQ Res function	`IDASetJacTimesResFn()`	NULL
Newton linear solve tolerance conversion factor	`IDASetLSNormFactor()`	vector length

The mathematical explanation of the linear solver methods available to IDA is provided in §6.2.2. We group the user-callable routines into four categories: general routines concerning the overall IDALS linear solver interface, optional inputs for matrix-based linear solvers, optional inputs for matrix-free linear solvers, and optional inputs for iterative linear solvers. We note that the matrix-based and matrix-free groups are mutually exclusive, whereas the “iterative” tag can apply to either case.

When using matrix-based linear solver modules, the IDALS solver interface needs a function to compute an approximation to the Jacobian matrix \(J(t,y,\dot{y})\). This function must be of type IDALsJacFn. The user can supply a Jacobian function or, if using a SUNMATRIX_DENSE or SUNMATRIX_BAND matrix \(J\), can use the default internal difference quotient approximation that comes with the IDALS interface. To specify a user-supplied Jacobian function jac, IDALS provides the function IDASetJacFn(). The IDALS interface passes the pointer user_data to the Jacobian function. This allows the user to create an arbitrary structure with relevant problem data and access it during the execution of the user-supplied Jacobian function, without using global data in the program. The pointer user_data may be specified through IDASetUserData().

int IDASetJacFn(void *ida_mem, IDALsJacFn jac)

The function IDASetJacFn specifies the Jacobian approximation function to be used for a matrix-based solver within the IDALS interface.

Arguments:

ida_mem – pointer to the IDA solver object.
jac – user-defined Jacobian approximation function. See IDALsJacFn for more details.

Return value:

IDALS_SUCCESS – The optional value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver interface has not been initialized.

Notes:

This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver(). By default, IDALS uses an internal difference quotient function for the SUNMATRIX_DENSE and SUNMATRIX_BAND modules. If NULL is passed to jac, this default function is used. An error will occur if no jac is supplied when using other matrix types.

Warning

The previous routine IDADlsSetJacFn is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

When using a matrix-based linear solver the matrix information will be updated infrequently to reduce matrix construction and, with direct solvers, factorization costs. As a result the value of \(\alpha\) may not be current and a scaling factor is applied to the solution of the linear system to account for the lagged value of \(\alpha\). See §11.5.1 for more details. The function IDASetLinearSolutionScaling() can be used to disable this scaling when necessary, e.g., when providing a custom linear solver that updates the matrix using the current \(\alpha\) as part of the solve.

int IDASetDeltaCjLSetup(void *ida_mem, sunrealtype dcj)

The function IDASetDeltaCjLSetup specifies the parameter that determines the bounds on a change in \(c_j\) that require a linear solver setup call. If cj_current / cj_previous < (1 - dcj) / (1 + dcj) or cj_current / cj_previous > (1 + dcj) / (1 - dcj), the linear solver setup function is called.

If dcj is \(< 0\) or \(\geq 1\) then the default value (0.25) is used.

Arguments:

ida_mem – pointer to the IDA memory block.
dcj – the \(c_j\) change threshold.

Return value:

IDA_SUCCESS – The flag value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetLinearSolutionScaling(void *ida_mem, sunbooleantype onoff)

The function IDASetLinearSolutionScaling enables or disables scaling the linear system solution to account for a change in \(\alpha\) in the linear system. For more details see §11.5.1.

Arguments:

ida_mem – pointer to the IDA solver object.
onoff – flag to enable (SUNTRUE) or disable (SUNFALSE) scaling

Return value:

IDALS_SUCCESS – The flag value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver interface has not been initialized.
IDALS_ILL_INPUT – The attached linear solver is not matrix-based.

Notes:

This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver(). By default scaling is enabled with matrix-based linear solvers.

When using matrix-free linear solver modules, the IDALS solver interface requires a function to compute an approximation to the product between the Jacobian matrix \(J(t,y,\dot{y})\) and a vector \(v\). The user can supply a Jacobian-times-vector approximation function, or use the default internal difference quotient function that comes with the IDALS solver interface.

A user-defined Jacobian-vector product function must be of type IDALsJacTimesVecFn and can be specified through a call to IDASetJacTimes(). The evaluation and processing of any Jacobian-related data needed by the user’s Jacobian-vector product function may be done in the optional user-supplied function jtsetup (see §6.4.4.6 for specification details). The pointer user_data received through IDASetUserData() (or a pointer to NULL if user_data was not specified) is passed to the Jacobian-vector product setup and product functions, jtsetup and jtimes, each time they are called. This allows the user to create an arbitrary structure with relevant problem data and access it during the execution of the user-supplied functions without using global data in the program.

int IDASetJacTimes(void *ida_mem, IDALsJacTimesSetupFn jsetup, IDALsJacTimesVecFn jtimes)

The function IDASetJacTimes specifies the Jacobian-vector product setup and product functions.

Arguments:

ida_mem – pointer to the IDA solver object.
jtsetup – user-defined function to set up the Jacobian-vector product. See IDALsJacTimesSetupFn for more details. Pass NULL if no setup is necessary.
jtimes – user-defined Jacobian-vector product function. See IDALsJacTimesVecFn for more details.

Return value:

IDALS_SUCCESS – The optional value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.
IDALS_SUNLS_FAIL – An error occurred when setting up the system matrix-times-vector routines in the SUNLinearSolver object used by the IDALS interface.

Notes:

The default is to use an internal finite difference quotient for jtimes and to omit jtsetup. If NULL is passed to jtimes, these defaults are used. A user may specify non-NULL jtimes and NULL jtsetup inputs. This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver().

Warning

The previous routine IDASpilsSetJacTimes is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

When using the default difference-quotient approximation to the Jacobian-vector product, the user may specify the factor to use in setting increments for the finite-difference approximation, via a call to IDASetIncrementFactor().

int IDASetIncrementFactor(void *ida_mem, sunrealtype dqincfac)

The function IDASetIncrementFactor specifies the increment factor to be used in the difference-quotient approximation to the product \(Jv\). Specifically, \(Jv\) is approximated via the formula

\[Jv = \frac{1}\sigma\left[F(t,\tilde{y},\tilde{\dot{y}}) - F(t,y,\dot{y})\right],\]

where \(\tilde{y} = y + \sigma v\), \(\tilde{\dot{y}} = \dot{y} + c_j \sigma v\), \(c_j\) is a BDF parameter proportional to the step size, \(\sigma = \mathtt{dqincfac} \sqrt{N}\), and \(N\) is the number of equations in the DAE system.

Arguments:

ida_mem – pointer to the IDA solver object.
dqincfac – user-specified increment factor positive.

Return value:

IDALS_SUCCESS – The optional value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.
IDALS_ILL_INPUT – The specified value of dqincfac is \(\le 0\).

Notes:

The default value is 1.0. This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver().

Warning

The previous routine IDASpilsSetIncrementFactor() is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

Additionally, when using the internal difference quotient, the user may also optionally supply an alternative residual function for use in the Jacobian-vector product approximation by calling IDASetJacTimesResFn(). The alternative residual function should compute a suitable (and differentiable) approximation to the residual function provided to IDAInit(). For example, as done in [46] for an ODE in explicit form, the alternative function may use lagged values when evaluating a nonlinearity to avoid differencing a potentially non-differentiable factor.

int IDASetJacTimesResFn(void *ida_mem, IDAResFn jtimesResFn)

The function IDASetJacTimesResFn specifies an alternative DAE residual function for use in the internal Jacobian-vector product difference quotient approximation.

Arguments:

ida_mem – pointer to the IDA solver object.
jtimesResFn – is the function which computes the alternative DAE residual function to use in Jacobian-vector product difference quotient approximations. See IDAResFn for more details.

Return value:

IDALS_SUCCESS – The optional value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.
IDALS_ILL_INPUT – The internal difference quotient approximation is disabled.

Notes:

The default is to use the residual function provided to IDAInit() in the internal difference quotient. If the input resudual function is NULL, the default is used. This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver().

When using an iterative linear solver, the user may supply a preconditioning operator to aid in solution of the system. This operator consists of two user-supplied functions, psetup and psolve, that are supplied to IDA using the function IDASetPreconditioner(). The psetup function supplied to this routine should handle evaluation and preprocessing of any Jacobian data needed by the user’s preconditioner solve function, psolve. Both of these functions are fully specified in §6.4.4.7 and §6.4.4.8). The user data pointer received through IDASetUserData() (or NULL if a user data pointer was not specified) is passed to the psetup and psolve functions. This allows the user to create an arbitrary structure with relevant problem data and access it during the execution of the user-supplied preconditioner functions without using global data in the program.

int IDASetPreconditioner(void *ida_mem, IDALsPrecSetupFn psetup, IDALsPrecSolveFn psolve)

The function IDASetPreconditioner specifies the preconditioner setup and solve functions.

Arguments:

ida_mem – pointer to the IDA solver object.
psetup – user-defined function to set up the preconditioner. See IDALsPrecSetupFn for more details. Pass NULL if no setup is necessary.
psolve – user-defined preconditioner solve function. See IDALsPrecSolveFn for more details.

Return value:

IDALS_SUCCESS – The optional values have been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.
IDALS_SUNLS_FAIL – An error occurred when setting up preconditioning in the SUNLinearSolver object used by the IDALS interface.

Notes:

The default is NULL for both arguments (i.e., no preconditioning). This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver().

Warning

The previous routine IDASpilsSetPreconditioner is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

Also, as described in §6.2.2, the IDALS interface requires that iterative linear solvers stop when the norm of the preconditioned residual satisfies

\[\|r\| \le \frac{\epsilon_L \epsilon}{10}\]

where \(\epsilon\) is the nonlinear solver tolerance, and the default \(\epsilon_L = 0.05\); this value may be modified by the user through the IDASetEpsLin() function.

int IDASetEpsLin(void *ida_mem, sunrealtype eplifac)

The function IDASetEpsLin specifies the factor by which the Krylov linear solver’s convergence test constant is reduced from the nonlinear iteration test constant.

Arguments:

ida_mem – pointer to the IDA solver object.
eplifac – linear convergence safety factor \(\geq 0.0\).

Return value:

IDALS_SUCCESS – The optional value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.
IDALS_ILL_INPUT – The factor eplifac is negative.

Notes:

The default value is \(0.05\). This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver(). If eplifac \(= 0.0\) is passed, the default value is used.

Warning

The previous routine IDASpilsSetEpsLin is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDASetLSNormFactor(void *ida_mem, sunrealtype nrmfac)

The function IDASetLSNormFactor specifies the factor to use when converting from the integrator tolerance (WRMS norm) to the linear solver tolerance (L2 norm) for Newton linear system solves e.g., tol_L2 = fac * tol_WRMS.

Arguments:

ida_mem – pointer to the IDA solver object.
nrmfac – the norm conversion factor.
- If nrmfac > 0, the provided value is used.
- If nrmfac = 0 then the conversion factor is computed using the vector length i.e., nrmfac = N_VGetLength(y) (default).
- If nrmfac < 0 then the conversion factor is computed using the vector dot product nrmfac = N_VDotProd(v,v) where all the entries of v are one.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

This function must be called after the IDALS linear solver interface has been initialized through a call to IDASetLinearSolver(). Prior to the introduction of N_VGetLength() in SUNDIALS v5.0.0 (IDA v5.0.0) the value of nrmfac was computed using N_VDotProd() i.e., the nrmfac < 0 case.

6.4.3.10.3. Nonlinear solver interface optional input functions

Table 6.3 Optional inputs for the IDANLS nonlinear solver interface
Optional input	Function name	Default
Maximum no. of nonlinear iterations	`IDASetMaxNonlinIters()`	4
Maximum no. of convergence failures	`IDASetMaxConvFails()`	10
Coeff. in the nonlinear convergence test	`IDASetNonlinConvCoef()`	0.33
Residual function for nonlinear system evaluations	`IDASetNlsResFn()`	`NULL`

The following functions can be called to set optional inputs controlling the nonlinear solver.

int IDASetMaxNonlinIters(void *ida_mem, int maxcor)

The function IDASetMaxNonlinIters specifies the maximum number of nonlinear solver iterations in one solve attempt.

Arguments:

ida_mem – pointer to the IDA solver object.
maxcor – maximum number of nonlinear solver iterations allowed in one solve attempt (>0).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_MEM_FAIL – The SUNNonlinearSolver object is NULL.

Notes:

The default value is 4.

int IDASetMaxConvFails(void *ida_mem, int maxncf)

The function IDASetMaxConvFails specifies the maximum number of nonlinear solver convergence failures in one step.

Arguments:

ida_mem – pointer to the IDA solver object.
maxncf – maximum number of allowable nonlinear solver convergence failures in one step (>0).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The default value is 10.

int IDASetNonlinConvCoef(void *ida_mem, sunrealtype nlscoef)

The function IDASetNonlinConvCoef specifies the safety factor in the nonlinear convergence test; see (6.7).

Arguments:

ida_mem – pointer to the IDA solver object.
nlscoef – coefficient in nonlinear convergence test (>0.0).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The value of nlscoef is \(\leq 0.0\).

Notes:

The default value is 0.33.

int IDASetNlsResFn(void *ida_mem, IDAResFn res)

The function IDASetNlsResFn specifies an alternative residual function for use in nonlinear system function evaluations.

Arguments:

ida_mem – pointer to the IDA solver object.
res – the alternative function which computes the DAE residual function \(F(t, y, \dot{y})\). See IDAResFn for more details.

Return value:

IDA_SUCCESS – The optional function has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The default is to use the residual function provided to IDAInit() in nonlinear system function evaluations. If the input residual function is NULL, the default is used.

When using a non-default nonlinear solver, this function must be called after IDASetNonlinearSolver().

6.4.3.10.4. Initial condition calculation optional input functions

Table 6.4 Optional inputs for IDA initial condition calculation
Optional input	Function name	Default
Coeff. in the nonlinear convergence test	`IDASetNonlinConvCoefIC()`	0.0033
Maximum no. of steps	`IDASetMaxNumStepsIC()`	5
Maximum no. of Jacobian/precond. evals.	`IDASetMaxNumJacsIC()`	4
Maximum no. of Newton iterations	`IDASetMaxNumItersIC()`	10
Max. linesearch backtracks per Newton iter.	`IDASetMaxBacksIC()`	100
Turn off linesearch	`IDASetLineSearchOffIC()`	`SUNFALSE`
Lower bound on Newton step	`IDASetStepToleranceIC()`	uround\(^{2/3}\)

The following functions can be called just prior to calling IDACalcIC() to set optional inputs controlling the initial condition calculation.

int IDASetNonlinConvCoefIC(void *ida_mem, sunrealtype epiccon)

The function IDASetNonlinConvCoefIC specifies the positive constant in the Newton iteration convergence test within the initial condition calculation.

Arguments:

ida_mem – pointer to the IDA solver object.
epiccon – coefficient in the Newton convergence test \((>0)\).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The epiccon factor is \(\leq 0.0\).

Notes:

The default value is \(0.01 \cdot 0.33\). This test uses a weighted RMS norm (with weights defined by the tolerances). For new initial value vectors \(y\) and \(\dot{y}\) to be accepted, the norm of \(J^{-1}F(t_0, y, \dot{y})\) must be \(\leq \mathtt{epiccon}\), where \(J\) is the system Jacobian.

int IDASetMaxNumStepsIC(void *ida_mem, int maxnh)

The function IDASetMaxNumStepsIC specifies the maximum number of steps allowed when icopt = IDA_YA_YDP_INIT in IDACalcIC(), where \(h\) appears in the system Jacobian, \(J = \dfrac{\partial F}{\partial y} + \left(\dfrac1h\right)\dfrac{\partial F}{\partial \dot{y}}\).

Arguments:

ida_mem – pointer to the IDA solver object.
maxnh – maximum allowed number of values for \(h\).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – maxnh is non-positive.

Notes:

The default value is \(5\).

int IDASetMaxNumJacsIC(void *ida_mem, int maxnj)

The function IDASetMaxNumJacsIC specifies the maximum number of the approximate Jacobian or preconditioner evaluations allowed when the Newton iteration appears to be slowly converging.

Arguments:

ida_mem – pointer to the IDA solver object.
maxnj – maximum allowed number of Jacobian or preconditioner evaluations.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – maxnj is non-positive.

Notes:

The default value is \(4\).

int IDASetMaxNumItersIC(void *ida_mem, int maxnit)

The function IDASetMaxNumItersIC specifies the maximum number of Newton iterations allowed in any one attempt to solve the initial conditions calculation problem.

Arguments:

ida_mem – pointer to the IDA solver object.
maxnit – maximum number of Newton iterations.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – maxnit is non-positive.

Notes:

The default value is \(10\).

int IDASetMaxBacksIC(void *ida_mem, int maxbacks)

The function IDASetMaxBacksIC specifies the maximum number of linesearch backtracks allowed in any Newton iteration, when solving the initial conditions calculation problem.

Arguments:

ida_mem – pointer to the IDA solver object.
maxbacks – maximum number of linesearch backtracks per Newton step.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – maxbacks is non-positive.

Notes:

The default value is \(100\).

int IDASetLineSearchOffIC(void *ida_mem, sunbooleantype lsoff)

The function IDASetLineSearchOffIC specifies whether to turn on or off the linesearch algorithm.

Arguments:

ida_mem – pointer to the IDA solver object.
lsoff – a flag to turn off (SUNTRUE) or keep (SUNFALSE) the linesearch algorithm.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The default value is SUNFALSE.

int IDASetStepToleranceIC(void *ida_mem, int steptol)

The function IDASetStepToleranceIC specifies a positive lower bound on the Newton step.

Arguments:

ida_mem – pointer to the IDA solver object.
steptol – Minimum allowed WRMS-norm of the Newton step \((> 0.0)\).

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – The steptol tolerance is \(\leq 0.0\).

Notes:

The default value is \((\text{unit roundoff})^{2/3}\).

6.4.3.10.5. Time step adaptivity optional input functions

Table 6.5 Optional inputs for IDA time step adaptivity
Optional input	Function name	Default
Fixed step size bounds \(\eta_{\mathrm{min\_fx}}\) and \(\eta_{\mathrm{max\_fx}}\)	`IDASetEtaFixedStepBounds()`	1.0 and 2.0
Maximum step size growth factor \(\eta_{\mathrm{max}}\)	`IDASetEtaMax()`	2.0
Minimum step size reduction factor \(\eta_{\mathrm{min}}\)	`IDASetEtaMin()`	0.5
Maximum step size reduction factor \(\eta_{\mathrm{low}}\)	`IDASetEtaLow()`	0.9
Minimum step size reduction factor after an error test failure \(\eta_{\mathrm{min\_ef}}\)	`IDASetEtaMinErrFail()`	0.25
Step size reduction factor after a nonlinear solver convergence failure \(\eta_{\mathrm{cf}}\)	`IDASetEtaConvFail()`	0.25

The following functions can be called to set optional inputs to control the step size adaptivity.

Note

The default values for the step size adaptivity tuning parameters have a long history of success and changing the values is generally discouraged. However, users that wish to experiment with alternative values should be careful to make changes gradually and with testing to determine their effectiveness.

int IDASetEtaFixedStepBounds(void *ida_mem, sunrealtype eta_min_fx, sunrealtype eta_max_fx)

The function IDASetEtaFixedStepBounds specifies the bounds \(\eta_{\mathrm{min\_fx}}\) and \(\eta_{\mathrm{max\_fx}}\). If step size change factor \(\eta\) satisfies \(\eta_{\mathrm{min\_fx}} < \eta < \eta_{\mathrm{max\_fx}}\) the current step size is retained.

The default values are \(\eta_{\mathrm{fxmin}} = 1\) and \(\eta_{\mathrm{fxmax}} = 2\).

eta_fxmin should satisfy \(0 < \eta_{\mathrm{fxmin}} \leq 1\), otherwise the default value is used. eta_fxmax should satisfy \(\eta_{\mathrm{fxmin}} \geq 1\), otherwise the default value is used.

Arguments:

ida_mem – pointer to the IDA solver object.
eta_min_fx – value of the fixed step size lower bound.
eta_max_fx – value of the fixed step size upper bound.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetEtaMax(void *ida_mem, sunrealtype eta_max)

The function IDASetEtaMax specifies the maximum step size growth factor \(\eta_{\mathrm{max}}\).

The default value is \(\eta_{\mathrm{max}} = 2\).

Arguments:

ida_mem – pointer to the IDA solver object.
eta_max – maximum step size growth factor.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetEtaMin(void *ida_mem, sunrealtype eta_min)

The function IDASetEtaMin specifies the minimum step size reduction factor \(\eta_{\mathrm{min}}\).

The default value is \(\eta_{\mathrm{min}} = 0.5\).

eta_min should satisfy \(0 < \eta_{\mathrm{min}} < 1\), otherwise the default value is used.

Arguments:

ida_mem – pointer to the IDA solver object.
eta_min – value of the minimum step size reduction factor.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetEtaLow(void *ida_mem, sunrealtype eta_low)

The function IDASetEtaLow specifies the maximum step size reduction factor \(\eta_{\mathrm{low}}\).

The default value is \(\eta_{\mathrm{low}} = 0.9\).

eta_low should satisfy \(0 < \eta_{\mathrm{low}} \leq 1\), otherwise the default value is used.

Arguments:

ida_mem – pointer to the IDA solver object.
eta_low – value of the maximum step size reduction factor.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetEtaMinErrFail(void *ida_mem, sunrealtype eta_min_ef)

The function IDASetEtaMinErrFail specifies the minimum step size reduction factor \(\eta_{\mathrm{min\_ef}}\) after an error test failure.

The default value is \(\eta_{\mathrm{min\_ef}} = 0.25\).

If eta_min_ef is \(\leq 0\) or \(\geq 1\), the default value is used.

Arguments:

ida_mem – pointer to the IDA solver object.
eta_low – value of the minimum step size reduction factor.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

int IDASetEtaConvFail(void *ida_mem, sunrealtype eta_cf)

The function IDASetEtaConvFail specifies the step size reduction factor \(\eta_{\mathrm{cf}}\) after a nonlinear solver convergence failure.

The default value is \(\eta_{\mathrm{cf}} = 0.25\).

If eta_cf is \(\leq 0\) or \(\geq 1\), the default value is used.

Arguments:

ida_mem – pointer to the IDA solver object.
eta_low – value of the step size reduction factor.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.2.0.

6.4.3.10.6. Rootfinding optional input functions

Table 6.6 Optional inputs for IDA rootfinding
Optional input	Function name	Default
Direction of zero-crossing	`IDASetRootDirection()`	both
Disable rootfinding warnings	`IDASetNoInactiveRootWarn()`	none

The following functions can be called to set optional inputs to control the rootfinding algorithm.

int IDASetRootDirection(void *ida_mem, int *rootdir)

The function IDASetRootDirection specifies the direction of zero-crossings to be located and returned to the user.

Arguments:

ida_mem – pointer to the IDA solver object.
rootdir – state array of length nrtfn , the number of root functions \(g_i\) , as specified in the call to the function IDARootInit().
- A value of \(0\) for rootdir[i] indicates that crossing in either direction should be reported for \(g_i\).
- A value of \(+1\) or \(-1\) for rootdir[i] indicates that the solver should report only zero-crossings where \(g_i\) is increasing or decreasing, respectively.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – rootfinding has not been activated through a call to IDARootInit().

Notes:

The default behavior is to locate both zero-crossing directions.

int IDASetNoInactiveRootWarn(void *ida_mem)

The function IDASetNoInactiveRootWarn disables issuing a warning if some root function appears to be identically zero at the beginning of the integration.

Arguments:

ida_mem – pointer to the IDA solver object.

Return value:

IDA_SUCCESS – The optional value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

IDA will not report the initial conditions as a possible zero-crossing (assuming that one or more components \(g_i\) are zero at the initial time). However, if it appears that some \(g_i\) is identically zero at the initial time (i.e., \(g_i\) is zero at the initial time and after the first step), IDA will issue a warning which can be disabled with this optional input function.

6.4.3.11. Interpolated output function

An optional function IDAGetDky() is available to obtain additional output values. This function must be called after a successful return from IDASolve() and provides interpolated values of \(y\) or its derivatives of order up to the last internal order used for any value of \(t\) in the last internal step taken by IDA.

int IDAGetDky(void *ida_mem, sunrealtype t, int k, N_Vector dky)

The function IDAGetDky computes the interpolated values of the \(k^{th}\) derivative of \(y\) for any value of \(t\) in the last internal step taken by IDA. The value of \(k\) must be non-negative and smaller than the last internal order used. A value of \(0\) for \(k\) means that the \(y\) is interpolated. The value of \(t\) must satisfy \(t_n - h_u \le t \le t_n\), where \(t_n\) denotes the current internal time reached, and \(h_u\) is the last internal step size used successfully.

Arguments:

ida_mem – pointer to the IDA solver object.
t – time at which to interpolate.
k – integer specifying the order of the derivative of \(y\) wanted.
dky – vector containing the interpolated \(k^{th}\) derivative of \(y(t)\).

Return value:

IDA_SUCCESS – IDAGetDky succeeded.
IDA_MEM_NULL – The ida_mem argument was NULL.
IDA_BAD_T – t is not in the interval \([t_n - h_u , t_n]\).
IDA_BAD_K – k is not one of \({0, 1, \ldots, k_{\text{last}}}\).
IDA_BAD_DKY – dky is NULL.

Notes:

It is only legal to call the function IDAGetDky() after a successful return from IDASolve(). Functions IDAGetCurrentTime(), IDAGetLastStep() and IDAGetLastOrder() can be used to access \(t_n\), \(h_u\), and \(k_{\text{last}}\).

6.4.3.12. Optional output functions

IDA provides an extensive list of functions that can be used to obtain solver performance information. Table 6.7 lists all optional output functions in IDA, which are then described in detail in the remainder of this section.

Some of the optional outputs, especially the various counters, can be very useful in determining how successful the IDA solver is in doing its job. For example, the counters nsteps and nrevals provide a rough measure of the overall cost of a given run, and can be compared among runs with differing input options to suggest which set of options is most efficient. The ratio nniters/nsteps measures the performance of the nonlinear solver in solving the nonlinear systems at each time step; typical values for this range from 1.1 to 1.8. The ratio njevals/nniters (in the case of a matrix-based linear solver), and the ratio npevals/nniters (in the case of an iterative linear solver) measure the overall degree of nonlinearity in these systems, and also the quality of the approximate Jacobian or preconditioner being used. Thus, for example, njevals/nniters can indicate if a user-supplied Jacobian is inaccurate, if this ratio is larger than for the case of the corresponding internal Jacobian. The ratio nliters/nniters measures the performance of the Krylov iterative linear solver, and thus (indirectly) the quality of the preconditioner.

Table 6.7 Optional outputs for IDA, IDALS, and IDANLS
Optional output	Function name
Size of IDA real and integer workspace	`IDAGetWorkSpace()`
Cumulative number of internal steps	`IDAGetNumSteps()`
No. of calls to residual function	`IDAGetNumResEvals()`
No. of calls to linear solver setup function	`IDAGetNumLinSolvSetups()`
No. of local error test failures that have occurred	`IDAGetNumErrTestFails()`
No. of failed steps due to a nonlinear solver failure	`IDAGetNumStepSolveFails()`
Order used during the last step	`IDAGetLastOrder()`
Order to be attempted on the next step	`IDAGetCurrentOrder()`
Actual initial step size used	`IDAGetActualInitStep()`
Step size used for the last step	`IDAGetLastStep()`
Step size to be attempted on the next step	`IDAGetCurrentStep()`
Current internal time reached by the solver	`IDAGetCurrentTime()`
Suggested factor for tolerance scaling	`IDAGetTolScaleFactor()`
Error weight vector for state variables	`IDAGetErrWeights()`
Estimated local errors	`IDAGetEstLocalErrors()`
All IDA integrator statistics	`IDAGetIntegratorStats()`
No. of nonlinear solver iterations	`IDAGetNumNonlinSolvIters()`
No. of nonlinear convergence failures	`IDAGetNumNonlinSolvConvFails()`
IDA nonlinear solver statistics	`IDAGetNonlinSolvStats()`
User data pointer	`IDAGetUserData()`
Array showing roots found	`IDAGetRootInfo()`
No. of calls to user root function	`IDAGetNumGEvals()`
Print all statistics	`IDAPrintAllStats()`
Name of constant associated with a return flag	`IDAGetReturnFlagName()`
Number of backtrack operations	`IDAGetNumBacktrackOps()`
Corrected initial conditions	`IDAGetConsistentIC()`
Stored Jacobian of the DAE residual function	`IDAGetJac()`
\(c_j\) value used in the Jacobian evaluation	`IDAGetJacCj()`
Time at which the Jacobian was evaluated	`IDAGetJacTime()`
Step number at which the Jacobian was evaluated	`IDAGetJacNumSteps()`
Size of real and integer workspace	`IDAGetLinWorkSpace()`
No. of Jacobian evaluations	`IDAGetNumJacEvals()`
No. of residual calls for finite diff. Jacobian-vector evals.	`IDAGetNumLinResEvals()`
No. of linear iterations	`IDAGetNumLinIters()`
No. of linear convergence failures	`IDAGetNumLinConvFails()`
No. of preconditioner evaluations	`IDAGetNumPrecEvals()`
No. of preconditioner solves	`IDAGetNumPrecSolves()`
No. of Jacobian-vector setup evaluations	`IDAGetNumJTSetupEvals()`
No. of Jacobian-vector product evaluations	`IDAGetNumJtimesEvals()`
Last return from a linear solver function	`IDAGetLastLinFlag()`
Name of constant associated with a return flag	`IDAGetLinReturnFlagName()`

6.4.3.12.1. Main solver optional output functions

IDA provides several user-callable functions that can be used to obtain different quantities that may be of interest to the user, such as solver workspace requirements, solver performance statistics, as well as additional data from the IDA solver object (a suggested tolerance scaling factor, the error weight vector, and the vector of estimated local errors). Also provided are functions to extract statistics related to the performance of the nonlinear solver being used. As a convenience, additional extraction functions provide the optional outputs in groups. These optional output functions are described next.

int IDAGetWorkSpace(void *ida_mem, long int lenrw, long int leniw)

The function IDAGetWorkSpace returns the IDA real and integer workspace sizes.

Arguments:

ida_mem – pointer to the IDA solver object.
lenrw – number of real values in the IDA workspace.
leniw – number of integer values in the IDA workspace.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

In terms of the problem size \(N\), the maximum method order maxord, and the number of root functions nrtfn (see §6.4.3.8), the actual size of the real workspace, in sunrealtype words, is given by the following:

base value: \(\mathtt{lenrw} = 55 + (m + 6) * N_r + 3 * \mathtt{nrtfn}\);
with IDASVtolerances(): \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);
with constraint checking (see IDASetConstraints()): \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);
with id specified (see IDASetId()): \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);

where \(m = \max(\mathtt{maxord}, 3)\), and \(N_r\) is the number of real words in one N_Vector \((\approx N)\).

The size of the integer workspace (without distinction between int and long int words) is given by:

base value: \(\mathtt{leniw} = 38 + (m + 6) * N_i + \mathtt{nrtfn}\);
with IDASVtolerances(): \(\mathtt{leniw} = \mathtt{leniw} + N_i\);
with constraint checking: \(\mathtt{lenrw} = \mathtt{lenrw} + N_i\);
with id specified (see IDASetId()): \(\mathtt{lenrw} = \mathtt{lenrw} + N_i\);

where \(N_i\) is the number of integer words in one N_Vector (= 1 for the serial N_Vector and 2 * npes for the parallel N_Vector on npes processors). For the default value of maxord, with no rootfinding, no id, no constraints, and with no call to IDASVtolerances(), these lengths are given roughly by \(\mathtt{lenrw} = 55 + 11 * N\) and \(\mathtt{leniw} = 49\).

int IDAGetNumSteps(void *ida_mem, long int *nsteps)

The function IDAGetNumSteps returns the cumulative number of internal steps taken by the solver (total so far).

Arguments:

ida_mem – pointer to the IDA solver object.
nsteps – number of steps taken by IDA.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetNumResEvals(void *ida_mem, long int *nrevals)

The function IDAGetNumResEvals returns the number of calls to the user’s residual evaluation function.

Arguments:

ida_mem – pointer to the IDA solver object.
nrevals – number of calls to the user’s res function.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

The nrevals value returned by IDAGetNumResEvals() does not account for calls made to res from a linear solver or preconditioner module.

int IDAGetNumLinSolvSetups(void *ida_mem, long int *nlinsetups)

The function IDAGetNumLinSolvSetups returns the cumulative number of calls made to the linear solver’s setup function (total so far).

Arguments:

ida_mem – pointer to the IDA solver object.
nlinsetups – number of calls made to the linear solver setup function.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetNumErrTestFails(void *ida_mem, long int *netfails)

The function IDAGetNumErrTestFails returns the cumulative number of local error test failures that have occurred (total so far).

Arguments:

ida_mem – pointer to the IDA solver object.
netfails – number of error test failures.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetNumStepSolveFails(void *ida_mem, long int *ncnf)

Returns the number of failed steps due to a nonlinear solver failure.

Arguments:

ida_mem – pointer to the IDA solver object.
ncnf – number of step failures.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetLastOrder(void *ida_mem, int *klast)

The function IDAGetLastOrder returns the integration method order used during the last internal step.

Arguments:

ida_mem – pointer to the IDA solver object.
klast – method order used on the last internal step.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetCurrentOrder(void *ida_mem, int *kcur)

The function IDAGetCurrentOrder returns the integration method order to be used on the next internal step.

Arguments:

ida_mem – pointer to the IDA solver object.
kcur – method order to be used on the next internal step.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetLastStep(void *ida_mem, sunrealtype *hlast)

The function IDAGetLastStep returns the integration step size taken on the last internal step (if from IDASolve()), or the last value of the artificial step size \(h\) (if from IDACalcIC()).

Arguments:

ida_mem – pointer to the IDA solver object.
hlast – step size taken on the last internal step by IDA, or last artificial step size used in IDACalcIC() , whichever was called last.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetCurrentStep(void *ida_mem, sunrealtype *hcur)

The function IDAGetCurrentStep returns the integration step size to be attempted on the next internal step.

Arguments:

ida_mem – pointer to the IDA solver object.
hcur – step size to be attempted on the next internal step.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetActualInitStep(void *ida_mem, sunrealtype *hinused)

The function IDAGetActualInitStep returns the value of the integration step size used on the first step.

Arguments:

ida_mem – pointer to the IDA solver object.
hinused – actual value of initial step size.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

Even if the value of the initial integration step size was specified by the user through a call to IDASetInitStep(), this value might have been changed by IDA to ensure that the step size is within the prescribed bounds \((h_{min} \le h_0 \le h_{max})\), or to meet the local error test.

int IDAGetCurrentTime(void *ida_mem, sunrealtype *tcur)

The function IDAGetCurrentTime returns the current internal time reached by the solver.

Arguments:

ida_mem – pointer to the IDA solver object.
tcur – current internal time reached.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetTolScaleFactor(void *ida_mem, sunrealtype *tolsfac)

The function IDAGetTolScaleFactor returns a suggested factor by which the user’s tolerances should be scaled when too much accuracy has been requested for some internal step.

Arguments:

ida_mem – pointer to the IDA solver object.
tolsfac – suggested scaling factor for user tolerances.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetErrWeights(void *ida_mem, N_Vector eweight)

The function IDAGetErrWeights returns the solution error weights at the current time. These are the \(W_i\) given by (6.4) (or by the user’s IDAEwtFn).

Arguments:

ida_mem – pointer to the IDA solver object.
eweight – solution error weights at the current time.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Warning

The user must allocate space for eweight.

int IDAGetEstLocalErrors(void *ida_mem, N_Vector ele)

The function IDAGetEstLocalErrors returns the estimated local errors.

Arguments:

ida_mem – pointer to the IDA solver object.
ele – estimated local errors at the current time.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Warning

The user must allocate space for ele. The values returned in ele are only valid if IDASolve() returned a non-negative value.

Note

The ele vector, togther with the eweight vector from IDAGetErrWeights(), can be used to determine how the various components of the system contributed to the estimated local error test. Specifically, that error test uses the RMS norm of a vector whose components are the products of the components of these two vectors. Thus, for example, if there were recent error test failures, the components causing the failures are those with largest values for the products, denoted loosely as eweight[i]*ele[i].

int IDAGetIntegratorStats(void *ida_mem, long int *nsteps, long int *nrevals, long int *nlinsetups, long int *netfails, int *qlast, int *qcur, sunrealtype *hinused, sunrealtype *hlast, sunrealtype *hcur, sunrealtype *tcur)

The function IDAGetIntegratorStats returns the IDA integrator stats in one function call.

Arguments:

ida_mem – pointer to the IDA solver object.
nsteps – cumulative number of steps taken by IDA.
nrevals – cumulative number of calls to the user’s res functions.
nlinsetups – cumulative number of calls made to the linear solver setup function.
netfails – cumulative number of error test failures.
klast – method order used on the last internal step.
kcur – method order to be used on the next internal step.
hinused – actual value of initial step size.
hlast – step sized taken on the last internal step.
hcur – step size to be attempted on the next internal step.
tcur – current internal time reached.

Return value:

IDA_SUCCESS – The optional output values have been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetNumNonlinSolvIters(void *ida_mem, long int *nniters)

The function IDAGetNumNonlinSolvIters returns the cumulative number of nonlinear iterations performed.

Arguments:

ida_mem – pointer to the IDA solver object.
nniters – number of nonlinear iterations performed.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_MEM_FAIL – The SUNNonlinearSolver object is NULL.

int IDAGetNumNonlinSolvConvFails(void *ida_mem, long int *nncfails)

The function IDAGetNumNonlinSolvConvFails returns the cumulative number of nonlinear convergence failures that have occurred.

Arguments:

ida_mem – pointer to the IDA solver object.
nncfails – number of nonlinear convergence failures.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetNonlinSolvStats(void *ida_mem, long int *nniters, long int *nncfails)

The function IDAGetNonlinSolvStats returns the IDA nonlinear solver statistics as a group.

Arguments:

ida_mem – pointer to the IDA solver object.
nniters – cumulative number of nonlinear iterations performed.
nncfails – cumulative number of nonlinear convergence failures.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_MEM_FAIL – The SUNNonlinearSolver object is NULL.

int IDAGetUserData(void *ida_mem, void **user_data)

The function IDAGetUserData returns the user data pointer provided to IDASetUserData().

Arguments:

ida_mem – pointer to the IDA memory block.
user_data – memory reference to a user data pointer.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

New in version 6.3.0.

int IDAPrintAllStats(void *ida_mem, FILE *outfile, SUNOutputFormat fmt)

The function IDAPrintAllStats outputs all of the integrator, nonlinear solver, linear solver, and other statistics.

Arguments:

ida_mem – pointer to the IDA memory block.
outfile – pointer to output file.
fmt – the output format:
- SUN_OUTPUTFORMAT_TABLE – prints a table of values
- SUN_OUTPUTFORMAT_CSV – prints a comma-separated list of key and value pairs e.g., key1,value1,key2,value2,...

Return value:

IDA_SUCCESS – The output was successfully.
IDA_MEM_NULL – The ida_mem pointer is NULL.
IDA_ILL_INPUT – An invalid formatting option was provided.

Note

The file scripts/sundials_csv.py provides python utility functions to read and output the data from a SUNDIALS CSV output file using the key and value pair format.

New in version 6.2.0.

char *IDAGetReturnFlagName(long int flag)

The function IDAGetReturnFlagName returns the name of the IDA constant corresponding to flag.

Arguments:

flag – the flag returned by a call to an IDA function

Return value:

char* – the flag name string

6.4.3.12.2. Initial condition calculation optional output functions

int IDAGetNumBacktrackOps(void *ida_mem, long int *nbacktr)

The function IDAGetNumBacktrackOps returns the number of backtrack operations done in the linesearch algorithm in IDACalcIC().

Arguments:

ida_mem – pointer to the IDA solver object.
nbacktr – the cumulative number of backtrack operations.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

int IDAGetConsistentIC(void *ida_mem, N_Vector yy0_mod, N_Vector yp0_mod)

The function IDAGetConsistentIC returns the corrected initial conditions calculated by IDACalcIC().

Arguments:

ida_mem – pointer to the IDA solver object.
yy0_mod – consistent solution vector.
yp0_mod – consistent derivative vector.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_ILL_INPUT – The function was not called before the first call to IDASolve().
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

If the consistent solution vector or consistent derivative vector is not desired, pass NULL for the corresponding argument.

Warning

The user must allocate space for yy0_mod and yp0_mod (if not NULL).

6.4.3.12.3. Rootfinding optional output functions

There are two optional output functions associated with rootfinding.

int IDAGetRootInfo(void *ida_mem, int *rootsfound)

The function IDAGetRootInfo returns an array showing which functions were found to have a root.

Arguments:

ida_mem – pointer to the IDA solver object.
rootsfound – array of length nrtfn with the indices of the user functions \(g_i\) found to have a root. For \(\mathtt{i} = 0, \ldots, \mathtt{nrtfn} -1\), \(\mathtt{rootsfound[i]} \ne 0\) if \(g_i\) has a root, and \(= 0\) if not.

Return value:

IDA_SUCCESS – The optional output values have been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

Notes:

Note that, for the components \(g_i\) for which a root was found, the sign of rootsfound[i] indicates the direction of zero-crossing. A value of \(+1\) indicates that \(g_i\) is increasing, while a value of \(-1\) indicates a decreasing \(g_i\).

Warning

The user must allocate memory for the vector rootsfound.

int IDAGetNumGEvals(void *ida_mem, long int *ngevals)

The function IDAGetNumGEvals returns the cumulative number of calls to the user root function \(g\).

Arguments:

ida_mem – pointer to the IDA solver object.
ngevals – number of calls to the user’s function \(g\) so far.

Return value:

IDA_SUCCESS – The optional output value has been successfully set.
IDA_MEM_NULL – The ida_mem pointer is NULL.

6.4.3.12.4. IDALS linear solver interface optional output functions

The following optional outputs are available from the IDALS modules:

int IDAGetJac(void *ida_mem, SUNMatrix *J)

Returns the internally stored copy of the Jacobian matrix of the DAE residual function.

Parameters:

ida_mem – the IDA memory structure
J – the Jacobian matrix

Return values:

IDALS_SUCCESS – the output value has been successfully set
IDALS_MEM_NULL – ida_mem was NULL
IDALS_LMEM_NULL – the linear solver interface has not been initialized

Warning

With linear solvers that overwrite the input Jacobian matrix as part of the linear solver setup (e.g., performing an in-place LU factorization) the matrix returned by IDAGetJac() may differ from the matrix returned by the last Jacobian evaluation.

Warning

This function is provided for debugging purposes and the values in the returned matrix should not be altered.

int IDAGetJacCj(void *ida_mem, sunrealtype *cj_J)

Returns the \(c_j\) value used to compute the internally stored copy of the Jacobian matrix of the DAE residual function.

Parameters:

ida_mem – the IDA memory structure
cj_J – the \(c_j\) value used in the Jacobian was evaluation

Return values:

IDALS_SUCCESS – the output value has been successfully set
IDALS_MEM_NULL – ida_mem was NULL
IDALS_LMEM_NULL – the linear solver interface has not been initialized

int IDAGetJacTime(void *ida_mem, sunrealtype *t_J)

Returns the time at which the internally stored copy of the Jacobian matrix of the DAE residual function was evaluated.

Parameters:

ida_mem – the IDA memory structure
t_J – the time at which the Jacobian was evaluated

Return values:

IDALS_SUCCESS – the output value has been successfully set
IDALS_MEM_NULL – ida_mem was NULL
IDALS_LMEM_NULL – the linear solver interface has not been initialized

int IDAGetJacNumSteps(void *ida_mem, long int *nst_J)

Returns the value of the internal step counter at which the internally stored copy of the Jacobian matrix of the DAE residual function was evaluated.

Parameters:

ida_mem – the IDA memory structure
nst_J – the value of the internal step counter at which the Jacobian was evaluated

Return values:

IDALS_SUCCESS – the output value has been successfully set
IDALS_MEM_NULL – ida_mem was NULL
IDALS_LMEM_NULL – the linear solver interface has not been initialized

int IDAGetLinWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS)

The function IDAGetLinWorkSpace returns the sizes of the real and integer workspaces used by the IDALS linear solver interface.

Arguments:

ida_mem – pointer to the IDA solver object.
lenrwLS – the number of real values in the IDALS workspace.
leniwLS – the number of integer values in the IDALS workspace.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Notes:

The workspace requirements reported by this routine correspond only to memory allocated within this interface and to memory allocated by the SUNLinearSolver object attached to it. The template Jacobian matrix allocated by the user outside of IDALS is not included in this report.

Warning

The previous routines IDADlsGetWorkspace and IDASpilsGetWorkspace are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJacEvals(void *ida_mem, long int *njevals)

The function IDAGetNumJacEvals returns the cumulative number of calls to the IDALS Jacobian approximation function.

Arguments:

ida_mem – pointer to the IDA solver object.
njevals – the cumulative number of calls to the Jacobian function total so far.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDADlsGetNumJacEvals is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinResEvals(void *ida_mem, long int *nrevalsLS)

The function IDAGetNumLinResEvals returns the cumulative number of calls to the user residual function due to the finite difference Jacobian approximation or finite difference Jacobian-vector product approximation.

Arguments:

ida_mem – pointer to the IDA solver object.
nrevalsLS – the cumulative number of calls to the user residual function.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Notes:

The value nrevalsLS is incremented only if one of the default internal difference quotient functions is used.

Warning

The previous routines IDADlsGetNumRhsEvals and IDASpilsGetNumRhsEvals are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinIters(void *ida_mem, long int *nliters)

The function IDAGetNumLinIters returns the cumulative number of linear iterations.

Arguments:

ida_mem – pointer to the IDA solver object.
nliters – the current number of linear iterations.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumLinIters is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinConvFails(void *ida_mem, long int *nlcfails)

The function IDAGetNumLinConvFails returns the cumulative number of linear convergence failures.

Arguments:

ida_mem – pointer to the IDA solver object.
nlcfails – the current number of linear convergence failures.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumConvFails is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumPrecEvals(void *ida_mem, long int *npevals)

The function IDAGetNumPrecEvals returns the cumulative number of preconditioner evaluations, i.e., the number of calls made to psetup.

Arguments:

ida_mem – pointer to the IDA solver object.
npevals – the cumulative number of calls to psetup.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumPrecEvals is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumPrecSolves(void *ida_mem, long int *npsolves)

The function IDAGetNumPrecSolves returns the cumulative number of calls made to the preconditioner solve function, psolve.

Arguments:

ida_mem – pointer to the IDA solver object.
npsolves – the cumulative number of calls to psolve.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumPrecSolves is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJTSetupEvals(void *ida_mem, long int *njtsetup)

The function IDAGetNumJTSetupEvals returns the cumulative number of calls made to the Jacobian-vector product setup function jtsetup.

Arguments:

ida_mem – pointer to the IDA solver object.
njtsetup – the current number of calls to jtsetup.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumJTSetupEvals is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJtimesEvals(void *ida_mem, long int *njvevals)

The function IDAGetNumJtimesEvals returns the cumulative number of calls made to the Jacobian-vector product function, jtimes.

Arguments:

ida_mem – pointer to the IDA solver object.
njvevals – the cumulative number of calls to jtimes.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Warning

The previous routine IDASpilsGetNumJtimesEvals is now a wrapper for this routine, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetLastLinFlag(void *ida_mem, long int *lsflag)

The function IDAGetLastLinFlag returns the last return value from an IDALS routine.

Arguments:

ida_mem – pointer to the IDA solver object.
lsflag – the value of the last return flag from an IDALS function.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer is NULL.
IDALS_LMEM_NULL – The IDALS linear solver has not been initialized.

Notes:

If the IDALS setup function failed (i.e., IDASolve() returned IDA_LSETUP_FAIL) when using the SUNLINSOL_DENSE or SUNLINSOL_BAND modules, then the value of lsflag is equal to the column index (numbered from one) at which a zero diagonal element was encountered during the LU factorization of the (dense or banded) Jacobian matrix. If the IDALS setup function failed when using another SUNLinearSolver object, then lsflag will be SUNLS_PSET_FAIL_UNREC, SUNLS_ASET_FAIL_UNREC, or SUN_ERR_EXT_FAIL. If the IDALS solve function failed (IDASolve() returned IDA_LSOLVE_FAIL), lsflag contains the error return flag from the SUNLinearSolver object, which will be one of: SUN_ERR_ARG_CORRUPTRRUPT, indicating that the SUNLinearSolver memory is NULL; SUNLS_ATIMES_FAIL_UNREC, indicating an unrecoverable failure in the \(J*v\) function; SUNLS_PSOLVE_FAIL_UNREC, indicating that the preconditioner solve function psolve failed unrecoverably; SUNLS_GS_FAIL, indicating a failure in the Gram-Schmidt procedure (generated only in SPGMR or SPFGMR); SUNLS_QRSOL_FAIL, indicating that the matrix \(R\) was found to be singular during the QR solve phase (SPGMR and SPFGMR only); or SUN_ERR_EXT_FAIL, indicating an unrecoverable failure in an external iterative linear solver package.

Warning

The previous routines IDADlsGetLastFlag and IDASpilsGetLastFlag are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

char *IDAGetLinReturnFlagName(long int lsflag)

The function IDAGetLinReturnFlagName returns the name of the IDALS constant corresponding to lsflag.

Arguments:

flag – the flag returned by a call to an IDA function

Return value:

char* – the flag name string or if \(1 \leq \mathtt{lsflag} \leq N\) (LU factorization failed), this function returns “NONE”.

Warning

The previous routines IDADlsGetReturnFlagName and IDASpilsGetReturnFlagName are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

6.4.3.13. IDA reinitialization function

The function IDAReInit() reinitializes the main IDA solver for the solution of a new problem, where a prior call to IDAInit() has been made. The new problem must have the same size as the previous one. IDAReInit() performs the same input checking and initializations that IDAInit() does, but does no memory allocation, as it assumes that the existing internal memory is sufficient for the new problem. A call to IDAReInit() deletes the solution history that was stored internally during the previous integration. Following a successful call to IDAReInit(), call IDASolve() again for the solution of the new problem.

The use of IDAReInit() requires that the maximum method order, maxord, is no larger for the new problem than for the problem specified in the last call to IDAInit(). In addition, the same N_Vector module set for the previous problem will be reused for the new problem.

If there are changes to the linear solver specifications, make the appropriate calls to either the linear solver objects themselves, or to the IDALS interface routines, as described in §6.4.3.5.

If there are changes to any optional inputs, make the appropriate IDASet*** calls, as described in §6.4.3.10.1. Otherwise, all solver inputs set previously remain in effect.

One important use of the IDAReInit() function is in the treating of jump discontinuities in the residual function. Except in cases of fairly small jumps, it is usually more efficient to stop at each point of discontinuity and restart the integrator with a readjusted DAE model, using a call to IDAReInit(). To stop when the location of the discontinuity is known, simply make that location a value of \(t_{\text{out}}\). 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 residual function not incorporate the discontinuity, but rather have a smooth extention over the discontinuity, so that the step across it (and subsequent rootfinding, if used) can be done efficiently. Then use a switch within the residual function (communicated through user_data) 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). Similar comments apply if there is to be a jump in the dependent variable vector.

int IDAReInit(void *ida_mem, sunrealtype t0, N_Vector y0, N_Vector yp0)

The function IDAReInit provides required problem specifications and reinitializes IDA.

Arguments:

ida_mem – pointer to the IDA solver object.
t0 – is the initial value of \(t\).
y0 – is the initial value of \(y\).
yp0 – is the initial value of \(\dot{y}\).

Return value:

IDA_SUCCESS – The call to was successful.
IDA_MEM_NULL – The IDA solver object was not initialized through a previous call to IDACreate().
IDA_NO_MALLOC – Memory space for the IDA solver object was not allocated through a previous call to IDAInit().
IDA_ILL_INPUT – An input argument to IDAReInit() has an illegal value.

Notes:

All previously set options are retained but may be updated by calling the appropriate “Set” functions.

If an error occurred, IDAReInit() also sends an error message to the error handler function.

6.4.4. User-supplied functions

The user-supplied functions consist of one function defining the DAE residual, (optionally) a function that handles error and warning messages, (optionally) a function that provides the error weight vector, (optionally) one or two functions that provide Jacobian-related information for the linear solver, and (optionally) one or two functions that define the preconditioner for use in any of the Krylov iteration algorithms.

6.4.4.1. DAE residual function

The user must provide a function of type IDAResFn defined as follows:

typedef int (*IDAResFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)

This function computes the problem residual for given values of the independent variable \(t\), state vector \(y\), and derivative \(\dot{y}\).

Arguments:

tt – is the current value of the independent variable.
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the output residual vector \(F(t,y,\dot{y})\).
user_data – is a pointer to user data, the same as the user_data pointer parameter passed to IDASetUserData().

Return value:

An IDAResFn function type should return a value of \(0\) if successful, a positive value if a recoverable error occurred (e.g., yy has an illegal value), or a negative value if a nonrecoverable error occurred. In the last case, the integrator halts. If a recoverable error occurred, the integrator will attempt to correct and retry.

Notes:

A recoverable failure error return from the IDAResFn is typically used to flag a value of the dependent variable \(y\) that is “illegal” in some way (e.g., negative where only a non-negative value is physically meaningful). If such a return is made, IDA will attempt to recover (possibly repeating the nonlinear solve, or reducing the step size) in order to avoid this recoverable error return.

For efficiency reasons, the DAE residual function is not evaluated at the converged solution of the nonlinear solver. Therefore, in general, a recoverable error in that converged value cannot be corrected. (It may be detected when the residual function is called the first time during the following integration step, but a successful step cannot be undone.)

6.4.4.2. Error weight function

typedef int (*IDAEwtFn)(N_Vector y, N_Vector ewt, void *user_data)

This function computes the WRMS error weights for the vector \(y\).

Arguments:

y – is the value of the dependent variable vector at which the weight vector is to be computed.
ewt – is the output vector containing the error weights.
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

0 – if it the error weights were successfully set.
-1 – if any error occured.

Notes:

Allocation of memory for ewt is handled within IDA.

Warning

The error weight vector must have all components positive. It is the user’s responsiblity to perform this test and return -1 if it is not satisfied.

6.4.4.3. Rootfinding function

If a rootfinding problem is to be solved during the integration of the DAE system, the user must supply a function of type IDARootFn, defined as follows:

typedef int (*IDARootFn)(sunrealtype t, N_Vector y, N_Vector yp, sunrealtype *gout, void *user_data)

This function computes a vector-valued function \(g(t,y,\dot{y})\) such that the roots of the nrtfn components \(g_i(t,y,\dot{y})\) are to be found during the integration.

Arguments:

t – is the current value of the independent variable.
y – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\), the \(t-\text{derivative}\) of \(y\).
gout – is the output array, of length nrtfn, with components \(g_i(t,y,\dot{y})\).
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

0 if successful or non-zero if an error occured (in which case the integration is halted and IDASolve() returs IDA_RTFUNC_FAIL).

Notes:

Allocation of memory for gout is handled within IDA.

6.4.4.4. Jacobian construction (matrix-based linear solvers)

If a matrix-based linear solver module is used (i.e. a non-NULL SUNMatrix object was supplied to IDASetLinearSolver()), the user may provide a function of type IDALsJacFn defined as follows:

typedef int (*IDALsJacFn)(sunrealtype t, sunrealtype c_j, N_Vector y, N_Vector yp, N_Vector r, SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)

This function computes the Jacobian matrix \(J\) of the DAE system (or an approximation to it), defined by (6.6).

Arguments:

tt – is the current value of the independent variable \(t\).
cj – is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the current value of the residual vector \(F(t,y,\dot{y})\).
Jac – is the output (approximate) Jacobian matrix (of type SUNMatrix), \(J = \dfrac{\partial{F}}{\partial{y}} + cj ~ \dfrac{\partial{F}}{\partial{\dot{y}}}\).
user_data - is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().
tmp1, tmp2, and tmp3 – are pointers to memory allocated for variables of type N_Vector which can be used by IDALsJacFn function as temporary storage or work space.

Return value:

An IDALsJacFn should return \(0\) if successful, a positive value if a recoverable error occurred, or a negative value if a nonrecoverable error occurred.

In the case of a recoverable eror return, the integrator will attempt to recover by reducing the stepsize, and hence changing \(\alpha\) in (6.6).

Notes:

Information regarding the structure of the specific SUNMatrix structure (e.g., number of rows, upper/lower bandwidth, sparsity type) may be obtained through using the implementation-specific SUNMatrix interface functions (see Chapter §10 for details).

With direct linear solvers (i.e., linear solvers with type SUNLINEARSOLVER_DIRECT), the Jacobian matrix \(J(t,y,\dot{y})\) is zeroed out prior to calling the user-supplied Jacobian function so only nonzero elements need to be loaded into Jac.

With the default nonlinear solver (the native SUNDIALS Newton method), each call to the user’s IDALsJacFn function is preceded by a call to the IDAResFn user function with the same (tt, yy, yp) arguments. Thus the Jacobian function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual. In the case of a user-supplied or external nonlinear solver, this is also true if the residual function is evaluated prior to calling the linear solver setup function (see §12.1.4 for more information).

If the user’s IDALsJacFn function uses difference quotient approximations, it may need to access quantities not in the call list. These quantities may include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to ida_mem to user_data and then use the IDAGet* functions described in §6.4.3.12.1. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h.

dense:

A user-supplied dense Jacobian function must load the Neq \(\times\) Neq dense matrix Jac with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt, yy, yp). The accessor macros SM_ELEMENT_D and SM_COLUMN_D allow the user to read and write dense matrix elements without making explicit references to the underlying representation of the SUNMATRIX_DENSE type. SM_ELEMENT_D(J, i, j) references the (i, j)-th element of the dense matrix Jac (with i, j\(= 0\ldots \texttt{N}-1\)). This macro is meant for small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from \(1\) to \(N\), the Jacobian element \(J_{m,n}\) can be set using the statement SM_ELEMENT_D(J, m-1, n-1) = \(J_{m,n}\). Alternatively, SM_COLUMN_D(J, j) returns a pointer to the first element of the j-th column of Jac (with j\(= 0\ldots \texttt{N}-1\)), and the elements of the j-th column can then be accessed using ordinary array indexing. Consequently, \(J_{m,n}\) can be loaded using the statements col_n = SM_COLUMN_D(J, n-1); col_n[m-1] = \(J_{m,n}\). For large problems, it is more efficient to use SM_COLUMN_D than to use SM_ELEMENT_D. Note that both of these macros number rows and columns starting from \(0\). The SUNMATRIX_DENSE type and accessor macros are documented in §10.9.

banded:

A user-supplied banded Jacobian function must load the Neq \(\times\) Neq banded matrix Jac with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt, yy, yp). The accessor macros SM_ELEMENT_B, SM_COLUMN_B, and SM_COLUMN_ELEMENT_B allow the user to read and write banded matrix elements without making specific references to the underlying representation of the SUNMATRIX_BAND type. SM_ELEMENT_B(J, i, j) references the (i, j)-th element of the banded matrix Jac, counting from \(0\). This macro is meant for use in small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from \(1\) to \(\texttt{N}\) with \((m,n)\) within the band defined by mupper and mlower, the Jacobian element \(J_{m,n}\) can be loaded using the statement SM_ELEMENT_B(J, m-1, n-1) = \(J_{m,n}\). The elements within the band are those with -mupper \(\le\) m-n \(\le\) mlower. Alternatively, SM_COLUMN_B(J, j) returns a pointer to the diagonal element of the j-th column of Jac, and if we assign this address to sunrealtype *col_j, then the i-th element of the j-th column is given by SM_COLUMN_ELEMENT_B(col_j, i, j), counting from \(0\). Thus, for \((m,n)\) within the band, \(J_{m,n}\) can be loaded by setting col_n = SM_COLUMN_B(J, n-1); and SM_COLUMN_ELEMENT_B(col_n, m-1, n-1) = \(J_{m,n}\). The elements of the j-th column can also be accessed via ordinary array indexing, but this approach requires knowledge of the underlying storage for a band matrix of type SUNMATRIX_BAND. The array col_n can be indexed from \(-\)mupper to mlower. For large problems, it is more efficient to use SM_COLUMN_B and SM_COLUMN_ELEMENT_B than to use the SM_ELEMENT_B macro. As in the dense case, these macros all number rows and columns starting from \(0\). The SUNMATRIX_BAND type and accessor macros are documented in §10.12.

sparse:

A user-supplied sparse Jacobian function must load the Neq \(\times\) Neq compressed-sparse-column or compressed-sparse-row matrix Jac with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt, yy, yp). Storage for Jac already exists on entry to this function, although the user should ensure that sufficient space is allocated in Jac to hold the nonzero values to be set; if the existing space is insufficient the user may reallocate the data and index arrays as needed. The amount of allocated space in a SUNMATRIX_SPARSE object may be accessed using the macro SM_NNZ_S or the routine SUNSparseMatrix_NNZ. The SUNMATRIX_SPARSE type and accessor macros are documented in §10.14.

Warning

The previous function type IDADlsJacFn is identical to IDALsJacFn, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

6.4.4.5. Jacobian-vector product (matrix-free linear solvers)

If a matrix-free linear solver is to be used (i.e., a NULL-valued SUNMatrix was supplied to IDASetLinearSolver()), the user may provide a function of type IDALsJacTimesVecFn in the following form, to compute matrix-vector products \(Jv\). If such a function is not supplied, the default is a difference quotient approximation to these products.

typedef int (*IDALsJacTimesVecFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, sunrealtype cj, void *user_data, N_Vector tmp1, N_Vector tmp2)

This function computes the product \(Jv\) of the DAE system Jacobian \(J\) (or an approximation to it) and a given vector v, where \(J\) is defined by (6.6).

Arguments:

tt – is the current value of the independent variable.
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the current value of the residual vector \(F(t,y,\dot{y})\).
v – is the vector by which the Jacobian must be multiplied to the right.
Jv – is the computed output vector.
cj – is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().
tmp1 and tmp2 – are pointers to memory allocated for variables of type N_Vector which can be used by IDALsJacTimesVecFn as temporary storage or work space.

Return value:

The value returned by the Jacobian-times-vector function should be 0 if successful. A nonzero value indicates that a nonrecoverable error occurred.

Notes:

This function must return a value of \(Jv\) that uses an approximation to the current value of \(J\), i.e. as evaluated at the current \((t,y,\dot{y})\).

If the user’s IDALsJacTimesVecFn function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to ida_mem to user_data and then use the IDAGet* functions described in §6.4.3.12.1. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h.

Warning

The previous function type IDASpilsJacTimesVecFn is identical to IDALsJacTimesVecFn, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

6.4.4.6. Jacobian-vector product setup (matrix-free linear solvers)

If the user’s Jacobian-vector product function requires that any Jacobian-related data be preprocessed or evaluated, then this needs to be done in a user-supplied function of type IDALsJacTimesSetupFn, defined as follows:

typedef int (*IDALsJacTimesSetupFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, ealtype cj, void *user_data);

This function setups any data needed by \(Jv\) product function (see IDALsJacTimesVecFn).

Arguments:

tt – is the current value of the independent variable.
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the current value of the residual vector \(F(t,y,\dot{y})\).
cj – is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

The value returned by the Jacobian-vector setup function should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes:

Each call to the Jacobian-vector product setup function is preceded by a call to the IDAResFn user function with the same \((t,y,\dot{y})\) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual.

If the user’s IDALsJacTimesVecFn function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to ida_mem to user_data and then use the IDAGet* functions described in §6.4.3.12.1. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h.

Warning

The previous function type IDASpilsJacTimesSetupFn is identical to IDALsJacTimesSetupFn, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

6.4.4.7. Preconditioner solve (iterative linear solvers)

If a user-supplied preconditioner is to be used with a SUNLinearSolver solver module, then the user must provide a function to solve the linear system \(Pz = r\) where \(P\) is a left preconditioner matrix which approximates (at least crudely) the Jacobian matrix \(J = \dfrac{\partial{F}}{\partial{y}} + cj ~ \dfrac{\partial{F}}{\partial{\dot{y}}}\). This function must be of type IDALsPrecSolveFn, defined as follows:

typedef int (*IDALsPrecSolveFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector rvec, N_Vector zvec, sunrealtype cj, sunrealtype delta, void *user_data)

This function solves the preconditioning system \(Pz = r\).

Arguments:

tt – is the current value of the independent variable.
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the current value of the residual vector \(F(t,y,\dot{y})\).
rvec – is the right-hand side vector \(r\) of the linear system to be solved.
zvec – is the computed output vector.
cj – is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).
delta – is an input tolerance to be used if an iterative method is employed in the solution. In that case, the residual vector \(Res = r - P z\) of the system should be made less than delta in weighted \(l_2\) norm, i.e., \(\sqrt{\displaystyle\sum_i (Res_i \cdot ewt_i)^2 } <\) delta. To obtain the N_Vector ewt, call IDAGetErrWeights().
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

The value returned by the preconditioner solve function should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

6.4.4.8. Preconditioner setup (iterative linear solvers)

If the user’s preconditioner requires that any Jacobian-related data be evaluated or preprocessed, then this needs to be done in a user-supplied function of type IDALsPrecSetupFn, defined as follows:

typedef int (*IDALsPrecSetupFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, sunrealtype cj, void *user_data)

This function solves the preconditioning system \(Pz = r\).

Arguments:

tt – is the current value of the independent variable.
yy – is the current value of the dependent variable vector, \(y(t)\).
yp – is the current value of \(\dot{y}(t)\).
rr – is the current value of the residual vector \(F(t,y,\dot{y})\).
cj – is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

The value returned by the preconditioner setup function should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes:

With the default nonlinear solver (the native SUNDIALS Newton method), each call to the preconditioner setup function is preceded by a call to the IDAResFn user function with the same \((t,y,\dot{y})\) arguments. Thus the preconditioner setup function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual. In the case of a user-supplied or external nonlinear solver, this is also true if the residual function is evaluated prior to calling the linear solver setup function (see §12.1.4 for more information).

This function is not called in advance of every call to the preconditioner solve function, but rather is called only as often as needed to achieve convergence in the nonlinear solver.

If the user’s IDALsPrecSetupFn function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to ida_mem to user_data and then use the IDAGet* functions described in §6.4.3.12.1. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h.

6.4.5. Preconditioner modules

A principal reason for using a parallel DAE solver such as IDA lies in the solution of partial differential equations (PDEs). Moreover, the use of a Krylov iterative method for the solution of many such problems is motivated by the nature of the underlying linear system of equations (6.5) that must be solved at each time step. The linear algebraic system is large, sparse, and structured. However, if a Krylov iterative method is to be effective in this setting, then a nontrivial preconditioner needs to be used. Otherwise, the rate of convergence of the Krylov iterative method is usually unacceptably slow. Unfortunately, an effective preconditioner tends to be problem-specific.

However, we have developed one type of preconditioner that treats a rather broad class of PDE-based problems. It has been successfully used for several realistic, large-scale problems [70] and is included in a software module within the IDA package. This module works with the parallel vector module NVECTOR_PARALLEL and generates a preconditioner that is a block-diagonal matrix with each block being a band matrix. The blocks need not have the same number of super- and sub-diagonals, and these numbers may vary from block to block. This Band-Block-Diagonal Preconditioner module is called IDABBDPRE.

One way to envision these preconditioners is to think of the domain of the computational PDE problem as being subdivided into \(M\) non-overlapping sub-domains. Each of these sub-domains is then assigned to one of the \(M\) processors to be used to solve the DAE system. The basic idea is to isolate the preconditioning so that it is local to each processor, and also to use a (possibly cheaper) approximate residual function. This requires the definition of a new function \(G(t,y,\dot{y})\) which approximates the function \(F(t, y, \dot{y})\) in the definition of the DAE system (6.1). However, the user may set \(G = F\). Corresponding to the domain decomposition, there is a decomposition of the solution vectors \(y\) and \(\dot{y}\) into \(M\) disjoint blocks \(y_m\) and \(\dot{y}_m\), and a decomposition of \(G\) into blocks \(G_m\). The block \(G_m\) depends on \(y_m\) and \(\dot{y}_m\), and also on components of \(y_{m'}\) and \(\dot{y}_{m'}\) associated with neighboring sub-domains (so-called ghost-cell data). Let \(\bar{y}_m\) and \(\bar{\dot{y}}_m\) denote \(y_m\) and \(\dot{y}_m\) (respectively) augmented with those other components on which \(G_m\) depends. Then we have

\[G(t,y,\dot{y}) = [G_1(t,\bar{y}_1,\bar{\dot{y}}_1), G_2(t,\bar{y}_2,\bar{\dot{y}}_2), \ldots, G_M(t,\bar{y}_M,\bar{\dot{y}}_M)]^T ~,\]

and each of the blocks \(G_m(t,\bar{y}_m,\bar{\dot{y}}_m)\) is uncoupled from the others.

The preconditioner associated with this decomposition has the form

\[\begin{split}P= \begin{bmatrix} P_1 & & &\\ & P_2 & &\\ & & \ddots & \\ & & & P_M\end{bmatrix}\end{split}\]

where

\[P_m \approx \frac{\partial G_m}{\partial y_m} + \alpha \frac{\partial G_m}{\partial \dot{y}_m}\]

This matrix is taken to be banded, with upper and lower half-bandwidths mudq and mldq defined as the number of non-zero diagonals above and below the main diagonal, respectively. The difference quotient approximation is computed using mudq \(+\) mldq \(+ 2\) evaluations of \(G_m\), but only a matrix of bandwidth mukeep \(+\) mlkeep \(+ 1\) is retained.

Neither pair of parameters need be the true half-bandwidths of the Jacobians of the local block of \(G\), if smaller values provide a more efficient preconditioner. Such an efficiency gain may occur if the couplings in the DAE system outside a certain bandwidth are considerably weaker than those within the band. Reducing mukeep and mlkeep while keeping mudq and mldq at their true values, discards the elements outside the narrower band. Reducing both pairs has the additional effect of lumping the outer Jacobian elements into the computed elements within the band, and requires more caution and experimentation.

The solution of the complete linear system

\[Px = b\]

reduces to solving each of the equations

\[P_m x_m = b_m\]

and this is done by banded LU factorization of \(P_m\) followed by a banded backsolve.

Similar block-diagonal preconditioners could be considered with different treatment of the blocks \(P_m\). For example, incomplete LU factorization or an iterative method could be used instead of banded LU factorization.

6.4.5.1. A parallel band-block-diagonal preconditioner module

The IDABBDPRE module calls two user-provided functions to construct \(P\): a required function Gres (of type IDABBDLocalFn) which approximates the residual function \(G(t,y,\dot{y}) \approx F(t,y,\dot{y})\) and which is computed locally, and an optional function Gcomm (of type IDABBDCommFn) which performs all inter-process communication necessary to evaluate the approximate residual \(G\). These are in addition to the user-supplied residual function res. Both functions take as input the same pointer user_data as passed by the user to IDASetUserData() and passed to the user’s function res. The user is responsible for providing space (presumably within user_data) for components of yy and yp that are communicated by Gcomm from the other processors, and that are then used by Gres, which should not do any communication.

typedef int (*IDABBDLocalFn)(sunindextype Nlocal, sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector gval, void *user_data)

This Gres function computes \(G(t,y,\dot{y})\). It loads the vector gval as a function of tt, yy, and yp.

Arguments:

Nlocal – is the local vector length.
tt – is the value of the independent variable.
yy – is the dependent variable.
yp – is the derivative of the dependent variable.
gval – is the output vector.
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

An IDABBDLocalFn function type should return 0 to indicate success, 1 for a recoverable error, or -1 for a non-recoverable error.

Notes:

This function must assume that all inter-processor communication of data needed to calculate gval has already been done, and this data is accessible within user_data.

The case where \(G\) is mathematically identical to \(F\) is allowed.

typedef int (*IDABBDCommFn)(sunindextype Nlocal, sunrealtype tt, N_Vector yy, N_Vector yp, void *user_data)

This Gcomm function performs all inter-processor communications necessary for the execution of the Gres function above, using the input vectors yy and yp.

Arguments:

Nlocal – is the local vector length.
tt – is the value of the independent variable.
yy – is the dependent variable.
yp – is the derivative of the dependent variable.
gval – is the output vector.
user_data – is a pointer to user data, the same as the user_data parameter passed to IDASetUserData().

Return value:

An IDABBDCommFn function type should return 0 to indicate success, 1 for a recoverable error, or -1 for a non-recoverable error.

Notes:

The Gcomm function is expected to save communicated data in space defined within the structure user_data.

Each call to the Gcomm function is preceded by a call to the residual function res with the same \((t,y,\dot{y})\) arguments. Thus Gcomm can omit any communications done by res if relevant to the evaluation of Gres. If all necessary communication was done in res, then Gcomm = NULL can be passed in the call to IDABBDPrecInit().

Besides the header files required for the integration of the DAE problem (see §6.4.1), to use the IDABBDPRE module, the main program must include the header file ida_bbdpre.h which declares the needed function prototypes.

The following is a summary of the usage of this module and describes the sequence of calls in the user main program. Steps that are unchanged from the user main program presented in §6.4.2 are not bold.

Initialize parallel or multi-threaded environment (if appropriate)
Create the vector of initial values
Create matrix object (if appropriate)
Create linear solver object (if appropriate)

When creating the iterative linear solver object, specify the use of left preconditioning (SUN_PREC_LEFT) as IDA only supports left preconditioning.
Create nonlinear solver object (if appropriate)
Create IDA object
Initialize IDA solver
Specify integration tolerances
Attach the linear solver (if appropriate)
Set linear solver optional inputs (if appropriate)

Note that the user should not overwrite the preconditioner setup function or solve function through calls to IDASetPreconditioner() optional input function.
Initialize the IDABBDPRE preconditioner module

Call IDABBDPrecInit() to allocate memory and initialize the internal preconditioner data. The last two arguments of IDABBDPrecInit() are the two user-supplied functions described above.
Attach nonlinear solver module (if appropriate)
Set nonlinear solver optional inputs (if appropriate)
Specify rootfinding problem (optional)
Set optional inputs
Advance solution in time
Get optional outputs

Additional optional outputs associated with IDABBDPRE are available by way of two routines described below, IDABBDPrecGetWorkSpace() and IDABBDPrecGetNumGfnEvals().
Deallocate memory
Finalize MPI, if used

The user-callable functions that initialize or re-initialize the IDABBDPRE preconditioner module are described next.

int IDABBDPrecInit(void *ida_mem, sunindextype Nlocal, sunindextype mudq, sunindextype mldq, sunindextype mukeep, sunindextype mlkeep, sunrealtype dq_rel_yy, IDABBDLocalFn Gres, IDABBDCommFn Gcomm);

The function IDABBDPrecInit initializes and allocates (internal) memory for the IDABBDPRE preconditioner.

Arguments:

ida_mem – pointer to the IDA solver object.
Nlocal – local vector dimension.
mudq – upper half-bandwidth to be used in the difference-quotient Jacobian approximation.
mldq – lower half-bandwidth to be used in the difference-quotient Jacobian approximation.
mukeep – upper half-bandwidth of the retained banded approximate Jacobian block.
mlkeep – lower half-bandwidth of the retained banded approximate Jacobian block.
dq_rel_yy – the relative increment in components of y used in the difference quotient approximations. The default is \(\mathtt{dq\_rel\_yy} = \sqrt{\text{unit roundoff}}\) , which can be specified by passing \(\mathtt{dq\_rel\_yy} = 0.0\).
Gres – the function which computes the local residual approximation \(G(t,y,\dot{y})\).
Gcomm – the optional function which performs all inter-process communication required for the computation of \(G(t,y,\dot{y})\).

Return value:

IDALS_SUCCESS – The call was successful.
IDALS_MEM_NULL – The ida_mem pointer was NULL.
IDALS_MEM_FAIL – A memory allocation request has failed.
IDALS_LMEM_NULL – An IDALS linear solver memory was not attached.
IDALS_ILL_INPUT – The supplied vector implementation was not compatible with the block band preconditioner.

Notes:

If one of the half-bandwidths mudq or mldq to be used in the difference-quotient calculation of the approximate Jacobian is negative or exceeds the value Nlocal-1, it is replaced by 0 or Nlocal-1 accordingly. The half-bandwidths mudq and mldq need not be the true half-bandwidths of the Jacobian of the local block of \(G\), when smaller values may provide a greater efficiency. Also, the half-bandwidths mukeep and mlkeep of the retained banded approximate Jacobian block may be even smaller, to reduce storage and computation costs further. For all four half-bandwidths, the values need not be the same on every processor.

The IDABBDPRE module also provides a reinitialization function to allow for a sequence of problems of the same size, with the same linear solver choice, provided there is no change in local_N, mukeep, or mlkeep. After solving one problem, and after calling IDAReInit() to re-initialize IDA for a subsequent problem, a call to IDABBDPrecReInit() can be made to change any of the following: the half-bandwidths mudq and mldq used in the difference-quotient Jacobian approximations, the relative increment dq_rel_yy, or one of the user-supplied functions Gres and Gcomm. If there is a change in any of the linear solver inputs, an additional call to the “Set”routines provided by the SUNLinearSolver object, and/or one or more of the corresponding IDASet*** functions, must also be made (in the proper order).

int IDABBDPrecReInit(void *ida_mem, sunindextype mudq, sunindextype mldq, sunrealtype dq_rel_yy)

The function IDABBDPrecReInit reinitializes the IDABBDPRE preconditioner.

Arguments:

ida_mem – pointer to the IDA solver object.
mudq – upper half-bandwidth to be used in the difference-quotient Jacobian approximation.
Mldq – lower half-bandwidth to be used in the difference-quotient Jacobian approximation.
dq_rel_yy – the relative increment in components of y used in the difference quotient approximations. The default is \(\mathtt{dq\_rel\_yy} = \sqrt{\text{unit roundoff}}\) , which can be specified by passing \(\mathtt{dq\_rel\_yy} = 0.0\).

Return value:

IDALS_SUCCESS – The call was successful.
IDALS_MEM_NULL – The ida_mem pointer was NULL.
IDALS_LMEM_NULL – An IDALS linear solver memory was not attached.
IDALS_PMEM_NULL – The function IDABBDPrecInit() was not previously called.

Notes:

If one of the half-bandwidths mudq or mldq is negative or exceeds the value Nlocal - 1, it is replaced by 0 or Nlocal - 1, accordingly.

The following two optional output functions are available for use with the IDABBDPRE module:

int IDABBDPrecGetWorkSpace(void *ida_mem, long int *lenrwBBDP, long int *leniwBBDP)

The function IDABBDPrecGetWorkSpace returns the local sizes of the IDABBDPRE real and integer workspaces.

Arguments:

ida_mem – pointer to the IDA solver object.
lenrwBBDP – local number of real values in the IDABBDPRE workspace.
leniwBBDP – local number of integer values in the IDABBDPRE workspace.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer was NULL.
IDALS_PMEM_NULL – The IDABBDPRE preconditioner has not been initialized.

Notes:

The workspace requirements reported by this routine correspond only to memory allocated within the IDABBDPRE module (the banded matrix approximation, banded SUNLinearSolver object, temporary vectors). These values are local to each process. The workspaces referred to here exist in addition to those given by the corresponding function IDAGetLinWorkSpace().

int IDABBDPrecGetNumGfnEvals(void *ida_mem, long int *ngevalsBBDP)

The function IDABBDPrecGetNumGfnEvals returns the cumulative number of calls to the user Gres function due to the finite difference approximation of the Jacobian blocks used within IDABBDPRE’s preconditioner setup function.

Arguments:

ida_mem – pointer to the IDA solver object.
ngevalsBBDP – the cumulative number of calls to the user Gres function.

Return value:

IDALS_SUCCESS – The optional output value has been successfully set.
IDALS_MEM_NULL – The ida_mem pointer was NULL.
IDALS_PMEM_NULL – The IDABBDPRE preconditioner has not been initialized.

In addition to the ngevalsBBDP evaluations of Gres, the costs associated with IDABBDPRE also includes nlinsetups LU factorizations, nlinsetups calls to Gcomm, npsolves banded backsolve calls, and nrevalsLS residual function evaluations, where nlinsetups is an optional IDA output (see §6.4.3.12.1), and npsolves and nrevalsLS are linear solver optional outputs (see §6.4.3.12.4).