Module Cvode

module Cvode: sig .. end
Variable-step solution of ODE initial value problems with zero-crossing detection.

This module solves numerically problems of the form $\dot{y} = f(t, y)$, $y(t_0) = y_0$.

This documented interface is structured as follows.

  1. Linear solvers
  2. Tolerances
  3. Solver initialization and use
  4. Modifying the solver
  5. Querying the solver
  6. Additional root finding functions
  7. Exceptions

Author(s): Timothy Bourke (Inria/ENS), Jun Inoue (Inria/ENS), Marc Pouzet (UPMC/ENS/Inria)
Version: 2.7.0

type ('d, 'k) session = ('d, 'k) session 
A session with the CVODE solver.

An example session with Cvode (cvode_skel.ml):

(* 1. Define a right-hand-side function. *)
let f t y yd = yd.{0} <- y.{1}; yd.{1} <- -9.81

(* 2. Optionally define a root function. *)
let g t y gout = gout.{0} <- 1.0 -. y.{0}

(* 3. Set vector of initial values.
      The length of this vector determines the problem size. *)
let yd = Sundials.RealArray.of_list [ 10.0; 0.0 ]
let y = Nvector_serial.wrap yd

(* 4. Create and initialize a solver session.
      This will initialize a specific linear solver and the root-finding
      mechanism, if necessary. *)
let s = Cvode.(init Adams Functional
                    (SStolerances (1e-4, 1e-8))
                    f ~roots:(1, g) 0.0 y);;

(* 5. Set optional inputs, e.g.,
      call [set_*] functions to change solver parameters. *)
Cvode.set_stop_time s 10.0;;
Cvode.set_all_root_directions s Sundials.RootDirs.Increasing;;

(* 6. Advance the solution in time,
      by repeatedly calling [solve_normal] or [solve_one_step]. *)
let rec go (t, r) =
  Printf.printf "% .10e\t% .10e\t% .10e\n" t yd.{0} yd.{1};
  match r with
  | Cvode.Success -> go (Cvode.solve_normal s (t +. 0.5) y)
  | Cvode.RootsFound -> begin
        yd.{1} <- -0.8 *. yd.{1};
        Cvode.reinit s t y;
        go (t, Cvode.Success)
      end
  | Cvode.StopTimeReached -> ();;

Printf.printf "time\ty\ty'\n";;
go (0.0, Cvode.Success);;

(* 7. Get optional outputs,
      call the [get_*] functions to examine solver statistics. *)
let ns = Cvode.get_num_steps s
    

See sundials: Skeleton of main program
type [> Nvector_serial.kind ] serial_session = (Nvector_serial.data, [> Nvector_serial.kind ] as 'a) session 
Alias for sessions based on serial nvectors.

Linear solvers


type ('data, 'kind) linear_solver = ('data, 'kind) linear_solver 
Linear solvers used by Cvode.
See sundials: Linear Solver Specification Functions
type [> Nvector_serial.kind ] serial_linear_solver = (Nvector_serial.data, [> Nvector_serial.kind ] as 'a) linear_solver 
Alias for linear solvers that are restricted to serial nvectors.
type 'd double = 'd * 'd 
Workspaces with two temporary vectors.
type 'd triple = 'd * 'd * 'd 
Workspaces with three temporary vectors.
type ('t, 'd) jacobian_arg = ('t, 'd) jacobian_arg = {
   jac_t : float; (*
The independent variable.
*)
   jac_y : 'd; (*
The dependent variable vector.
*)
   jac_fy : 'd; (*
The derivative vector (i.e., $\frac{\mathrm{d}y}{\mathrm{d}t}$).
*)
   jac_tmp : 't; (*
Workspace data.
*)
}
Arguments common to Jacobian callback functions.
See sundials: CVDlsDenseJacFn
See sundials: CVDlsBandJacFn
See sundials: CVSpilsJacTimesVecFn
See sundials: CVSpilsPrecSolveFn
See sundials: CVSpilsPrecSetupFn
type bandrange = bandrange = {
   mupper : int; (*
The upper half-bandwidth.
*)
   mlower : int; (*
The lower half-bandwidth.
*)
}
The range of nonzero entries in a band matrix.
module Diag: sig .. end
Diagonal approximation of Jacobians by difference quotients.
module Dls: sig .. end
Direct Linear Solvers operating on dense and banded matrices.
module Sls: sig .. end
Sparse Linear Solvers.
module Spils: sig .. end
Scaled Preconditioned Iterative Linear Solvers.
module Alternate: sig .. end
Alternate Linear Solvers.

Tolerances


type 'data error_weight_fun = 'data -> 'data -> unit 
Functions that set the multiplicative error weights for use in the weighted RMS norm. The call efun y ewt takes the dependent variable vector y and fills the error-weight vector ewt with positive values or raises Sundials.NonPositiveEwt. Other exceptions are eventually propagated, but should be avoided (efun is not allowed to abort the solver).
type ('data, 'kind) tolerance = 
| SStolerances of float * float (*
(rel, abs) : scalar relative and absolute tolerances.
*)
| SVtolerances of float * ('data, 'kind) Nvector.t (*
(rel, abs) : scalar relative and vector absolute tolerances.
*)
| WFtolerances of 'data error_weight_fun (*
Set the multiplicative error weights for the weighted RMS norm.
*)
Tolerance specifications.
val default_tolerances : ('data, 'kind) tolerance
A default relative tolerance of 1.0e-4 and absolute tolerance of 1.0e-8.

Solver initialization and use


type ('d, 'kind) iter = 
| Newton of ('d, 'kind) linear_solver (*
Newton iteration with a given linear solver.
*)
| Functional (*
Functional iteration (non-stiff systems only).
*)
Choice of method for solving non-linear systems that arise in solver formulas.
See sundials: IVP Solution
See sundials: CVodeCreate
type lmm = 
| Adams (*
Adams-Moulton formulas (non-stiff systems).
*)
| BDF (*
Backward Differentiation Formulas (stiff systems).
*)
Choice of linear multistep method.
See sundials: IVP Solution
See sundials: CVodeCreate
type 'd rhsfn = float -> 'd -> 'd -> unit 
Right-hand side functions for calculating ODE derivatives. They are passed three arguments: Within the function, raising a Sundials.RecoverableFailure exception indicates a recoverable error. Any other exception is treated as an unrecoverable error.

y and y' should not be accessed after the function returns.

See sundials: CVRhsFn
type 'd rootsfn = float -> 'd -> Sundials.RealArray.t -> unit 
Called by the solver to calculate the values of root functions. These ‘zero-crossings’ are used to detect significant events. The function is passed three arguments:
y and gout should not be accessed after the function has returned.

See sundials: cvRootFn
val init : lmm ->
('data, 'kind) iter ->
('data, 'kind) tolerance ->
'data rhsfn ->
?roots:int * 'data rootsfn ->
float -> ('data, 'kind) Nvector.t -> ('data, 'kind) session
Creates and initializes a session with the solver. The call
init lmm iter tol f ~roots:(nroots, g) t0 y0
has as arguments: This function does everything necessary to initialize a session, i.e., it makes the calls referenced below. The Cvode.solve_normal and Cvode.solve_one_step functions may be called directly.
See sundials: CVodeCreate/CVodeInit
See sundials: CVodeRootInit
See sundials: Linear solvers
See sundials: CVodeSStolerances
See sundials: CVodeSVtolerances
See sundials: CVodeWFtolerances
See sundials: CVEwtFn
val no_roots : int * 'd rootsfn
A convenience value for signalling that there are no roots to monitor.
type solver_result = 
| Success (*
The solution was advanced. (CV_SUCCESS)
*)
| RootsFound (*
A root was found. See Cvode.get_root_info. (CV_ROOT_RETURN)
*)
| StopTimeReached (*
The stop time was reached. See Cvode.set_stop_time. (CV_TSTOP_RETURN)
*)
Values returned by the step functions. Failures are indicated by exceptions.
See sundials: CVode
val solve_normal : ('d, 'k) session ->
float -> ('d, 'k) Nvector.t -> float * solver_result
Integrates an ODE system over an interval. The call tret, r = solve_normal s tout yout has as arguments It returns tret, the time reached by the solver, which will be equal to tout if no errors occur, and, r, a Cvode.solver_result.
Raises
See sundials: CVode (CV_NORMAL)
val solve_one_step : ('d, 'k) session ->
float -> ('d, 'k) Nvector.t -> float * solver_result
Like Cvode.solve_normal but returns after one internal solver step.
See sundials: CVode (CV_ONE_STEP)
val get_dky : ('d, 'k) session -> ('d, 'k) Nvector.t -> float -> int -> unit
Returns the interpolated solution or derivatives. get_dky s dky t k computes the kth derivative of the function at time t, i.e., $\frac{d^\mathtt{k}y(\mathtt{t})}{\mathit{dt}^\mathtt{k}}$, and stores it in dky. The arguments must satisfy $t_n - h_u \leq \mathtt{t} \leq t_n$—where $t_n$ denotes Cvode.get_current_time and $h_u$ denotes Cvode.get_last_step,— and $0 \leq \mathtt{k} \leq q_u$—where $q_u$ denotes Cvode.get_last_order.

This function may only be called after a successful return from either Cvode.solve_normal or Cvode.solve_one_step.

Raises
See sundials: CVodeGetDky
val reinit : ('d, 'kind) session ->
?iter_type:('d, 'kind) iter ->
?roots:int * 'd rootsfn -> float -> ('d, 'kind) Nvector.t -> unit
Reinitializes the solver with new parameters and state values. The values of the independent variable, i.e., the simulation time, and the state variables must be given. If given, iter_type specifies a new iteration method, and roots specifies a new root finding function; both default to unchanged.
See sundials: CVodeReInit

Modifying the solver (optional input functions)


val set_tolerances : ('d, 'k) session -> ('d, 'k) tolerance -> unit
Sets the integration tolerances.
See sundials: CVodeSStolerances
See sundials: CVodeSVtolerances
See sundials: CVodeWFtolerances
See sundials: CVEwtFn
val set_error_file : ('d, 'k) session -> Sundials.Logfile.t -> unit
Configure the default error handler to write messages to a file. By default it writes to Sundials.Logfile.stderr.
See sundials: CVodeSetErrFile
val set_err_handler_fn : ('d, 'k) session -> (Sundials.error_details -> unit) -> unit
Specifies a custom function for handling error messages. The handler must not fail: any exceptions are trapped and discarded.
See sundials: CVodeSetErrHandlerFn
See sundials: CVErrHandlerFn
val clear_err_handler_fn : ('d, 'k) session -> unit
Restores the default error handling function.
See sundials: CVodeSetErrHandlerFn
val set_max_ord : ('d, 'k) session -> int -> unit
Specifies the maximum order of the linear multistep method.
See sundials: CVodeSetMaxOrd
val set_max_num_steps : ('d, 'k) session -> int -> unit
Specifies the maximum number of steps taken in attempting to reach a given output time.
See sundials: CVodeSetMaxNumSteps
val set_max_hnil_warns : ('d, 'k) session -> int -> unit
Specifies the maximum number of messages warning that t + h = t on the next internal step.
See sundials: CVodeSetMaxHnilWarns
val set_stab_lim_det : ('d, 'k) session -> bool -> unit
Indicates whether the BDF stability limit detection algorithm should be used.
See sundials: CVodeSetStabLimDet
See sundials: BDF Stability Limit Detection
val set_init_step : ('d, 'k) session -> float -> unit
Specifies the initial step size.
See sundials: CVodeSetInitStep
val set_min_step : ('d, 'k) session -> float -> unit
Specifies a lower bound on the magnitude of the step size.
See sundials: CVodeSetMinStep
val set_max_step : ('d, 'k) session -> float -> unit
Specifies an upper bound on the magnitude of the step size.
See sundials: CVodeSetMaxStep
val set_stop_time : ('d, 'k) session -> float -> unit
Limits the value of the independent variable t when solving. By default no stop time is imposed.
See sundials: CVodeSetStopTime
val set_max_err_test_fails : ('d, 'k) session -> int -> unit
Specifies the maximum number of error test failures permitted in attempting one step.
See sundials: CVodeSetMaxErrTestFails
val set_max_nonlin_iters : ('d, 'k) session -> int -> unit
Specifies the maximum number of nonlinear solver iterations permitted per step.
See sundials: CVodeSetMaxNonlinIters
val set_max_conv_fails : ('d, 'k) session -> int -> unit
Specifies the maximum number of nonlinear solver convergence failures permitted during one step.
See sundials: CVodeSetMaxConvFails
val set_nonlin_conv_coef : ('d, 'k) session -> float -> unit
Specifies the safety factor used in the nonlinear convergence test.
See sundials: CVodeSetNonlinConvCoef
See sundials: IVP Solution

Querying the solver (optional output functions)


val get_work_space : ('d, 'k) session -> int * int
Returns the real and integer workspace sizes.
Returns (real_size, integer_size)
See sundials: CVodeGetWorkSpace
val get_num_steps : ('d, 'k) session -> int
Returns the cumulative number of internal steps taken by the solver.
See sundials: CVodeGetNumSteps
val get_num_rhs_evals : ('d, 'k) session -> int
Returns the number of calls to the right-hand side function.
See sundials: CVodeGetNumRhsEvals
val get_num_lin_solv_setups : ('d, 'k) session -> int
Returns the number of calls made to the linear solver's setup function.
See sundials: CVodeGetNumLinSolvSetups
val get_num_err_test_fails : ('d, 'k) session -> int
Returns the number of local error test failures that have occurred.
See sundials: CVodeGetNumErrTestFails
val get_last_order : ('d, 'k) session -> int
Returns the integration method order used during the last internal step.
See sundials: CVodeGetLastOrder
val get_current_order : ('d, 'k) session -> int
Returns the integration method order to be used on the next internal step.
See sundials: CVodeGetCurrentOrder
val get_last_step : ('d, 'k) session -> float
Returns the integration step size taken on the last internal step.
See sundials: CVodeGetLastStep
val get_current_step : ('d, 'k) session -> float
Returns the integration step size to be attempted on the next internal step.
See sundials: CVodeGetCurrentStep
val get_actual_init_step : ('d, 'k) session -> float
Returns the the value of the integration step size used on the first step.
See sundials: CVodeGetActualInitStep
val get_current_time : ('d, 'k) session -> float
Returns the the current internal time reached by the solver.
See sundials: CVodeGetCurrentTime
val get_num_stab_lim_order_reds : ('d, 'k) session -> int
Returns the number of order reductions dictated by the BDF stability limit detection algorithm.
See sundials: CVodeGetNumStabLimOrderReds
See sundials: BDF stability limit detection
val get_tol_scale_factor : ('d, 'k) session -> float
Returns a suggested factor by which the user's tolerances should be scaled when too much accuracy has been requested for some internal step.
See sundials: CVodeGetTolScaleFactor
val get_err_weights : ('d, 'k) session -> ('d, 'k) Nvector.t -> unit
Returns the solution error weights at the current time.
See sundials: CVodeGetErrWeights
See sundials: IVP solution (W_i)
val get_est_local_errors : ('d, 'k) session -> ('d, 'k) Nvector.t -> unit
Returns the vector of estimated local errors.
See sundials: CVodeGetEstLocalErrors
type integrator_stats = {
   num_steps : int; (*
Cumulative number of internal solver steps.
*)
   num_rhs_evals : int; (*
Number of calls to the right-hand side function.
*)
   num_lin_solv_setups : int; (*
Number of setups calls to the linear solver.
*)
   num_err_test_fails : int; (*
Number of local error test failures.
*)
   last_order : int; (*
Integration method order used in the last internal step.
*)
   current_order : int; (*
Integration method order to be used in the next internal step.
*)
   actual_init_step : float; (*
Integration step sized used on the first step.
*)
   last_step : float; (*
Integration step size of the last internal step.
*)
   current_step : float; (*
Integration step size to attempt on the next internal step.
*)
   current_time : float; (*
Current internal time reached by the solver.
*)
}
Summaries of integrator statistics.
val get_integrator_stats : ('d, 'k) session -> integrator_stats
Returns the integrator statistics as a group.
See sundials: CVodeGetIntegratorStats
val print_integrator_stats : ('d, 'k) session -> Pervasives.out_channel -> unit
Prints the integrator statistics on the given channel.
See sundials: CVodeGetIntegratorStats
val get_num_nonlin_solv_iters : ('d, 'k) session -> int
Returns the number of nonlinear (functional or Newton) iterations performed.
See sundials: CVodeGetNumNonlinSolvIters
val get_num_nonlin_solv_conv_fails : ('d, 'k) session -> int
Returns the number of nonlinear convergence failures that have occurred.
See sundials: CVodeGetNumNonlinSolvConvFails
val get_nonlin_solv_stats : ('d, 'k) session -> int * int
Returns both the numbers of nonlinear iterations performed nniters and nonlinear convergence failures nncfails.
Returns (nniters, nncfails)
See sundials: CVodeGetNonlinSolvStats

Additional root-finding functions


val set_root_direction : ('d, 'k) session -> Sundials.RootDirs.d array -> unit
set_root_direction s dir specifies the direction of zero-crossings to be located and returned. dir may contain one entry for each root function.
See sundials: CVodeSetRootDirection
val set_all_root_directions : ('d, 'k) session -> Sundials.RootDirs.d -> unit
Like Cvode.set_root_direction but specifies a single direction for all root functions.
See sundials: CVodeSetRootDirection
val set_no_inactive_root_warn : ('d, 'k) session -> unit
Disables issuing a warning if some root function appears to be identically zero at the beginning of the integration.
See sundials: CVodeSetNoInactiveRootWarn
val get_num_roots : ('d, 'k) session -> int
Returns the number of root functions.
val get_root_info : ('d, 'k) session -> Sundials.Roots.t -> unit
Fills an array showing which functions were found to have a root.
See sundials: CVodeGetRootInfo
val get_num_g_evals : ('d, 'k) session -> int
Returns the cumulative number of calls made to the user-supplied root function g.
See sundials: CVodeGetNumGEvals

Exceptions


exception IllInput
Raised on missing or illegal solver inputs. Also raised if an element of the error weight vector becomes zero during time stepping, or the linear solver initialization function failed, or a root was found both at t and very near t.
See sundials: CV_ILL_INPUT
exception TooClose
The initial and final times are too close to each other and an initial step size was not specified.
See sundials: CV_TOO_CLOSE
exception TooMuchWork
The requested time could not be reached in mxstep internal steps. See Cvode.set_max_num_steps
See sundials: CV_TOO_MUCH_WORK
exception TooMuchAccuracy
The requested accuracy could not be satisfied.
See sundials: CV_TOO_MUCH_ACC
exception ErrFailure
Too many error test failures within a step or at the minimum step size. See Cvode.set_max_err_test_fails and Cvode.set_min_step.
See sundials: CV_ERR_FAILURE
exception ConvergenceFailure
Too many convergence test failures within a step or at the minimum step size. See Cvode.set_max_conv_fails and Cvode.set_min_step.
See sundials: CV_CONV_FAILURE
exception LinearInitFailure
Linear solver initialization failed.
See sundials: CV_LINIT_FAIL
exception LinearSetupFailure
Linear solver setup failed in an unrecoverable manner.
See sundials: CV_LSETUP_FAIL
exception LinearSolveFailure
Linear solver solution failed in an unrecoverable manner.
See sundials: CV_LSOLVE_FAIL
exception RhsFuncFailure
The right-hand side function failed in an unrecoverable manner.
See sundials: CV_RHSFUNC_FAIL
exception FirstRhsFuncFailure
The right-hand side function had a recoverable error when first called.
See sundials: CV_FIRST_RHSFUNC_ERR
exception RepeatedRhsFuncFailure
Too many convergence test failures, or unable to estimate the initial step size, due to repeated recoverable errors in the right-hand side function.
See sundials: CV_REPTD_RHSFUNC_ERR
exception UnrecoverableRhsFuncFailure
The right-hand side function had a recoverable error, but no recovery was possible. This error can only occur after an error test failure at order one.
See sundials: CV_UNREC_RHSFUNC_ERR
exception RootFuncFailure
The rootfinding function failed.
See sundials: CV_RTFUNC_FAIL
exception BadK
Raised by Cvode.get_dky for invalid order values.
See sundials: CVodeGetDky (CV_BAD_K)
exception BadT
Raised by Cvode.get_dky for invalid time values.
See sundials: CVodeGetDky (CV_BAD_T)