Previous  Up  Next

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
type ('d, 'k) session = ('d, 'k) session

A session with the CVODE solver.

An example session with Cvode (cvode_skel.ml):

open Sundials

(* 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 = 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
                    (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 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
    
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.

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.

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.

module Diag: sig .. end

Diagonal approximation of Jacobians by difference quotients.

module Dls: sig .. end

Direct Linear Solvers operating on dense, banded, and sparse matrices.

module Spils: sig .. end

Scaled Preconditioned Iterative Linear Solvers.

val matrix_embedded_solver : (unit, 'data, 'kind, [> `MatE ]) Sundials.LinearSolver.t ->
       ('data, 'kind) linear_solver

Create a CVode-specific linear solver from a generic matrix embedded solver.

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 lmm =
| Adams (*

Adams-Moulton formulas (non-stiff systems).

*)
| BDF (*

Backward Differentiation Formulas (stiff systems).

*)

Choice of linear multistep method.

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.

type 'd proj_fn = float -> 'd -> 'd -> float -> 'd option -> unit

A function to compute the projection of the solution and, if enabled, the error on the constraint manifold.

Such functions take the following arguments:

(The err argument is None if error projection is not enabled.)

Raising Sundials.RecoverableFailure indicates a recoverable error. Any other exception is treated as an unrecoverable error. The integrator will, in most cases, try to correct and reattempt the step.

The solve should stop when the WRMS norm of the current iterate update is less than eps. The projection routine can access the error weight vector with Cvode.get_err_weights.

val init : ?context:Sundials.Context.t ->
       lmm ->
       ('data, 'kind) tolerance ->
       ?nlsolver:('data, 'kind, ('data, 'kind) session, [ `Nvec ])
                 Sundials_NonlinearSolver.t ->
       ?nlsrhsfn:'data rhsfn ->
       ?lsolver:('data, 'kind) linear_solver ->
       'data rhsfn ->
       ?roots:int * 'data rootsfn ->
       ?projfn:'data proj_fn ->
       float -> ('data, 'kind) Nvector.t -> ('data, 'kind) session

Creates and initializes a session with the solver. The call

init lmm tol ~nlsolver ~nlsrhsfn ~lsolver f ~roots:(nroots, g) ~projfn 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.

If an nlsolver is not specified, then the Newton module is used by default. In this case only, lsolver defaults to Cvode.Diag.solver if not otherwise specified. Specifying an nlsolver that requires a linear solver without specifying an lsolver results in a Cvode.NonlinearInitFailure (or Cvode.IllInput for Sundials < 4.0.0) exception on the first call to Cvode.solve_normal or Cvode.solve_one_step.

The projection feature is only supported for Sundials >= 5.3.0 and the BDF method.

The alternative right-hand-side function for nonlinear system function evaluations is only supported for Sundials >= 5.8.0.

By default, the session is created using the context returned by Sundials.Context.default, but this can be overridden by passing an optional context argument.

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.

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.

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.

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.

val reinit : ('d, 'k) session ->
       ?nlsolver:('d, 'k, ('d, 'k) session, [ `Nvec ])
                 Sundials_NonlinearSolver.t ->
       ?nlsrhsfn:'d rhsfn ->
       ?lsolver:('d, 'k) linear_solver ->
       ?roots:int * 'd rootsfn ->
       ?rhsfn:'d rhsfn -> float -> ('d, 'k) 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, nlsolver specifies a nonlinear solver, nlsrhsfn specifies an alternative rhs function for nonlinear system function evaluations,lsolver specifies a linear solver, roots specifies a new root finding function, and rhsfn specifies a new rhs function; all default to unchanged.

If the new problem does not have a constraint equation, but the old one did, then Cvode.set_proj_frequency must a zero argument to disable projection.

Modifying the solver (optional input functions)

val set_tolerances : ('d, 'k) session -> ('d, 'k) tolerance -> unit

Sets the integration tolerances.

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 Logfile.stderr.

val set_err_handler_fn : ('d, 'k) session -> (Sundials.Util.error_details -> unit) -> unit

Specifies a custom function for handling error messages. The handler must not fail: any exceptions are trapped and discarded.

val clear_err_handler_fn : ('d, 'k) session -> unit

Restores the default error handling function.

val set_monitor_fn : ('d, 'k) session -> int -> (('d, 'k) session -> unit) -> unit

Specifies a function to be called after the given number of successful steps.

The solver solution may be read by the monitoring function, but it should not be changed.

This function requires that the underlying library was explicitly built with support for monitoring (see Sundials_Config.monitoring_enabled).

val set_monitor_frequency : ('d, 'k) session -> int -> unit

Sets the number of successful steps between calls to the monitoring function.

val clear_monitor_fn : ('d, 'k) session -> unit

Turns monitoring off.

val set_max_ord : ('d, 'k) session -> int -> unit

Specifies the maximum order of the linear multistep method.

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.

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.

val set_stab_lim_det : ('d, 'k) session -> bool -> unit

Indicates whether the BDF stability limit detection algorithm should be used.

val set_init_step : ('d, 'k) session -> float -> unit

Specifies the initial step size.

val set_min_step : ('d, 'k) session -> float -> unit

Specifies a lower bound on the magnitude of the step size.

val set_max_step : ('d, 'k) session -> float -> unit

Specifies an upper bound on the magnitude of the step size.

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.

val set_max_err_test_fails : ('d, 'k) session -> int -> unit

Specifies the maximum number of error test failures permitted in attempting one step.

val set_max_nonlin_iters : ('d, 'k) session -> int -> unit

Specifies the maximum number of nonlinear solver iterations permitted per step.

val set_max_conv_fails : ('d, 'k) session -> int -> unit

Specifies the maximum number of nonlinear solver convergence failures permitted during one step.

val set_nonlin_conv_coef : ('d, 'k) session -> float -> unit

Specifies the safety factor used in the nonlinear convergence test.

val set_constraints : ('d, 'k) session -> ('d, 'k) Nvector.t -> unit

Specifies a vector defining inequality constraints for each component of the solution vector y. See Sundials.Constraint.

val clear_constraints : ('d, 'k) session -> unit

Disables constraint checking.

val set_proj_err_est : ('d, 'k) session -> bool -> unit

Enables or disables projection of the error estimate by the projection function.

val set_proj_frequency : ('d, 'k) session -> int -> unit

Set the frequency with which the projection is performed. The default is 1, that is, every time step. A value of 0 disables projection and a value less than zero restores the default.

val set_max_num_proj_fails : ('d, 'k) session -> int -> unit

Set the maximum number of projection failures in a step attempt before an unrecoverable error is returned. The default is 10. A value less than 1 restores the default.

val set_eps_proj : ('d, 'k) session -> float -> unit

Set the tolerance for the nonlinear-constrained least-squares problem solved by the projection function. The default is 0.1. A value less than or equal to zero restores the default.

val set_proj_fail_eta : ('d, 'k) session -> float -> unit

Sets the time-step reduction factor to apply on a projection function failure. The default is 0.25. A value less than or equal to 1, or greater than 1 restores the default.

Querying the solver (optional output functions)

val get_work_space : ('d, 'k) session -> int * int

Returns the real and integer workspace sizes.

val get_num_steps : ('d, 'k) session -> int

Returns the cumulative number of internal steps taken by the solver.

val get_num_rhs_evals : ('d, 'k) session -> int

Returns the number of calls to the right-hand side function.

val get_num_lin_solv_setups : ('d, 'k) session -> int

Returns the number of calls made to the linear solver's setup function.

val get_num_err_test_fails : ('d, 'k) session -> int

Returns the number of local error test failures that have occurred.

val get_last_order : ('d, 'k) session -> int

Returns the integration method order used during the last internal step.

val get_current_order : ('d, 'k) session -> int

Returns the integration method order to be used on the next internal step.

val get_current_state : ('d, 'k) session -> 'd

Returns the current state vector. This vector provides direct access to the data within the integrator.

type 'd nonlin_system_data = {
   tn : float; (*

Independent variable value $t_n$ .

*)
   ypred : 'd; (*

Predicted state vector $y_{\mathit{pred}}$ at $t_n$ . This data must not be changed.

*)
   yn : 'd; (*

State vector $y^n$ . This data may not be current and may need to be filled.

*)
   fn : 'd; (*

The right-hand side function evaluated at the current time and state, $f(t_n, y^n)$ . * This data may not be current and may need to be filled.

*)
   gamma : float; (*

Current value of $\gamma$ .

*)
   rl1 : float; (*

A scaling factor used to compute $\tilde{a}_n = \mathtt{rl1}\cdot\mathtt{zn1}$ .

*)
   zn1 : 'd; (*

A vector used to compute $\tilde{a}_n = \mathtt{rl1}\cdot\mathtt{zn1}$ .

*)
}

Internal data required to construct the current nonlinear implicit system within a nonlinear solver.

val get_nonlin_system_data : ('d, 'k) session -> 'd nonlin_system_data

Gives direct access to the internal data required to construct the current nonlinear system within a nonlinear solver. This function should be called inside the nonlinear system function. If the nonlinear solver uses the lsetup or lsolve functions, then the nonlinear solver system function must fill the zi and fi vectors with, respectively, the current state and corresponding evaluation of the right-hand-side function: $y^n = y_{\mathit{pred}} + y_{\mathit{cor}}$ and $f_n = f(t_n, y^n)$ where $y_{\mathit{cor}}$ is the first argument of the nonlinear solver system function. Within a custom linear solver, then the vectors yn and fn are only current after an evaluation of the nonlinear system function.

val compute_state : ('d, 'k) session -> ('d, 'k) Nvector.t -> ('d, 'k) Nvector.t -> unit

Computes the current stage state vector using the stored prediction and the supplied correction from the nonlinear solver. The call compute_state s ycor yn computes $y^n = y_{\mathit{pred}} + y_{\mathit{cor}}$ .

val get_current_gamma : ('d, 'k) session -> float

Returns the current value of $\gamma$ . This scalar appears in the internal Newton equation, $M = I - \gamma J$ .

val get_last_step : ('d, 'k) session -> float

Returns the integration step size taken on the last internal step.

val get_current_step : ('d, 'k) session -> float

Returns the integration step size to be attempted on the next internal step.

val get_actual_init_step : ('d, 'k) session -> float

Returns the the value of the integration step size used on the first step.

val get_current_time : ('d, 'k) session -> float

Returns the the current internal time reached by the solver.

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.

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.

val get_err_weights : ('d, 'k) session -> ('d, 'k) Nvector.t -> unit

Returns the solution error weights at the current time.

val get_est_local_errors : ('d, 'k) session -> ('d, 'k) Nvector.t -> unit

Returns the vector of estimated local errors.

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.

val print_integrator_stats : ('d, 'k) session -> Stdlib.out_channel -> unit

Prints the integrator statistics on the given channel.

type linear_solver_stats = {
   jac_evals : int; (*

Number of calls made by a linear solver to the Jacobian approximation function.

*)
   lin_rhs_evals : int; (*

Number of calls to the right-hand side callback due to the finite difference Jacobian approximation.

*)
   lin_iters : int; (*

The cumulative number of linear iterations.

*)
   lin_conv_fails : int; (*

The cumulative number of linear convergence failures.

*)
   prec_evals : int; (*

The cumulative number of calls to the setup function.

*)
   prec_solves : int; (*

The cumulative number of calls to the solve function.

*)
   jtsetup_evals : int; (*

The cumulative number of calls to the Jacobian-vector setup function.

*)
   jtimes_evals : int; (*

The cumulative number of calls to the Jacobian-vector function.

*)
}

Summaries of linear solver statistics.

val get_linear_solver_stats : ('d, 'k) session -> linear_solver_stats

Returns linear solver statistics as a group.

val get_num_nonlin_solv_iters : ('d, 'k) session -> int

Returns the cumulative number of nonlinear (functional or Newton) iterations.

val get_num_nonlin_solv_conv_fails : ('d, 'k) session -> int

Returns the cumulative number of nonlinear convergence failures.

val get_nonlin_solv_stats : ('d, 'k) session -> int * int

Returns both the numbers of nonlinear iterations performed nniters and nonlinear convergence failures nncfails.

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.

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.

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.

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.

val get_num_g_evals : ('d, 'k) session -> int

Returns the cumulative number of calls made to the user-supplied root function g.

val get_num_proj_evals : ('d, 'k) session -> int

Returns the current total number of projection evaluations.

val get_num_proj_fails : ('d, 'k) session -> int

Returns the current total number of projection evaluation failures.

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.

exception TooClose

The initial and final times are too close to each other and an initial step size was not specified.

exception TooMuchWork

The requested time could not be reached in mxstep internal steps. See Cvode.set_max_num_steps

exception TooMuchAccuracy

The requested accuracy could not be satisfied.

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.

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.

exception LinearInitFailure

Linear solver initialization failed.

exception LinearSetupFailure of exn option

Linear solver setup failed in an unrecoverable manner. If possible, the exception in the underlying linear solver is specified. It is typically one of Sundials_LinearSolver.ZeroInDiagonal, Sundials_LinearSolver.PSetFailure, or Sundials_LinearSolver.PackageFailure.

exception LinearSolveFailure of exn option

Linear solver solution failed in an unrecoverable manner. If possible, the exception in the underlying linear solver is specified. It is typically one of Sundials_LinearSolver.ZeroInDiagonal, Sundials_LinearSolver.ATimesFailure, Sundials_LinearSolver.PSolveFailure, Sundials_LinearSolver.GSFailure, Sundials_LinearSolver.QRSolFailure, or Sundials_LinearSolver.PackageFailure.

exception NonlinearSolverFailure

The nonlinear solver failed in a general way.

exception NonlinearInitFailure

Nonlinear solver initialization failed.

exception NonlinearSetupFailure

Nonlinear solver setup failed in an unrecoverable manner.

exception RhsFuncFailure

The right-hand side function failed in an unrecoverable manner.

exception FirstRhsFuncFailure

The right-hand side function had a recoverable error when first called.

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.

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.

exception RootFuncFailure

The rootfinding function failed.

exception ConstraintFailure

No solution satisfying the inequality constraints could be found.

exception BadK

Raised by Cvode.get_dky for invalid order values.

exception BadT

Raised by Cvode.get_dky for invalid time values.

exception VectorOpErr

A fused vector operation failed.

exception ProjFuncFailure

The projection function failed.

exception RepeatedProjFuncError

The projection function failed repeatedly.

exception ProjectionNotEnabled

The project functionality is not enabled. A projection function must be given in the call to Cvode.init and the last call to Cvode.set_proj_frequency must not have set the frequency to zero.