Module Arkode

module Arkode: sig .. end

Adaptive-step time integration for stiff, nonstiff, and multi-rate systems of ODE initial value problems with zero-crossing detection.

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

This documented interface is structured as follows.

  1. Linear and mass matrix 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 ARKODE solver.

An example session with Arkode (arkode_skel.ml):

open Sundials

(* 1. Define right-hand-side functions. *)
let f_e t y yd = yd.{0} <- y.{1}
let f_i t y yd = 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 determine whether the problem is purely explicit, purely
      implicit, or both. Initialize a specific linear solver if there is an
      implicit component, and root-finding if necessary. It is also possible to
      specify a mass matrix solver. *)
let m = Matrix.dense 2
let s = Arkode.(
  init
    (ImEx { explicit = f_e;
            implicit = (f_i, Newton Arkode.Dls.(solver (dense y m)),
                        Linear true); })
    (SStolerances (1e-4, 1e-9))
    ~roots:(1, g)
    0.0
    y);;

(* 5. Set optional inputs, e.g.,
      call [set_*] functions to change solver parameters. *)
Arkode.set_stop_time s 10.0;;
Arkode.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
  | Arkode.Success -> go (Arkode.solve_normal s (t +. 0.5) y)
  | Arkode.RootsFound -> begin
        yd.{1} <- -0.8 *. yd.{1};
        Arkode.reinit s t y;
        go (t, Arkode.Success)
      end
  | Arkode.StopTimeReached -> ();;

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

(* 7. Get optional outputs,
      call the [get_*] functions to examine solver statistics. *)
let ns = Arkode.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 and mass matrix solvers

type ('data, 'kind) session_linear_solver = ('data, 'kind) session_linear_solver 

Linear solvers used by Arkode.

type [> Nvector_serial.kind ] serial_session_linear_solver = (Nvector_serial.data, [> Nvector_serial.kind ] as 'a)
session_linear_solver

Alias for linear solvers that are restricted to serial nvectors.

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 $f_I(t, y)$.

*)
   jac_tmp : 't; (*

Workspace data.

*)
}

Arguments common to Jacobian callback functions.

module Dls: sig .. end

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

module Spils: sig .. end

Scaled Preconditioned Iterative Linear Solvers.

module Alternate: sig .. end

Alternate Linear Solvers.

module Mass: sig .. end

Mass Matrix 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-9.

type 'data res_weight_fun = 'data -> 'data -> unit 

Functions that compute the weighted RMS residual weights. The call rfun y rwt takes the dependent variable vector y and fills the residual-weight vector rwt with positive values or raises Sundials.NonPositiveEwt. Other exceptions are eventually propagated, but should be avoided (ffun is not allowed to abort the solver).

type ('data, 'kind) res_tolerance = 
| ResStolerance of float (*

abs : scalar absolute residual tolerance.

*)
| ResVtolerance of ('data, 'kind) Nvector.t (*

abs : vector of absolute residual tolerances.

*)
| ResFtolerance of 'data res_weight_fun (*

Compute the residual weight vector.

*)

Solver initialization and use

type ('d, 'kind) iter = 
| Newton of ('d, 'kind) session_linear_solver (*

Modified Newton iteration with a given linear solver.

*)
| FixedPoint of int (*

Accelerated fixed-point solver. Specifies the number of vectors to store within the Anderson acceleration subspace.

*)

Choice of method for solving non-linear systems that arise in the solution of implicit systems.

type linearity = 
| Linear of bool (*

Implicit portion is linear. Specifies whether $f_I(t, y)$ or the preconditioner is time-dependent.

*)
| Nonlinear (*

Implicit portion is nonlinear.

*)

The linearity of the implicit portion of the problem.

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.

type ('d, 'k) imex = {
   implicit : 'd rhsfn * ('d, 'k) iter * linearity;
   explicit : 'd rhsfn;
}
type ('d, 'k) problem = 
| Implicit of 'd rhsfn * ('d, 'k) iter * linearity (*

Diagonally Implicit Runge-Kutta (DIRK) solution of stiff problem.

*)
| Explicit of 'd rhsfn (*

Explicit Runge-Kutta (ERK) solution of non-stiff problem.

*)
| ImEx of ('d, 'k) imex (*

Additive Runge-Kutta (ARK) solution of multi-rate problem.

*)

The form of the initial value problem.

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.

val no_roots : int * 'd rootsfn

A convenience value for signalling that there are no roots to monitor.

val init : ('data, 'kind) problem ->
('data, 'kind) tolerance ->
?restol:('data, 'kind) res_tolerance ->
?order:int ->
?mass:('data, 'kind) Mass.solver ->
?roots:int * 'data rootsfn ->
float -> ('data, 'kind) Nvector.t -> ('data, 'kind) session

Creates and initializes a session with the solver. The call

init problem tol ~restol ~order:ord ~mass:msolver ~roots:(nroots, g) t0 y0

has as arguments:

The allowed values for ord are:

This function does everything necessary to initialize a session, i.e., it makes the calls referenced below. The Arkode.solve_normal and Arkode.solve_one_step functions may be called directly.

type solver_result = 
| Success (*

The solution was advanced. (ARK_SUCCESS)

*)
| RootsFound (*

A root was found. See Arkode.get_root_info. (ARK_ROOT_RETURN)

*)
| StopTimeReached (*

The stop time was reached. See Arkode.set_stop_time. (ARK_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 Arkode.solver_result.

val solve_one_step : ('d, 'k) session ->
float -> ('d, 'k) Nvector.t -> float * solver_result

Like Arkode.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}}{\mathit{dt}^\mathtt{k}}y(\mathtt{t})$, and stores it in dky. The arguments must satisfy $t_n - h_n \leq \mathtt{t} \leq t_n$—where $t_n$ denotes Arkode.get_current_time and $h_n$ denotes Arkode.get_last_step,— and $0 \leq \mathtt{k} \leq 3$.

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

val reinit : ('d, 'kind) session ->
?problem:('d, 'kind) problem ->
?order:int ->
?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. The new problem must have the same size as the previous one. Only the functions in the given problem description have any effect; the iter and linear fields are ignored.

The allowed values for order are:

type 'd resize_fn = 'd -> 'd -> unit 

Called to resize a vector to match the dimensions of another. The call resizefn y ytemplate must resize y to match the size of ytemplate.

y and ytemplate should not be accessed after the function has returned.

val resize : ('d, 'kind) session ->
?resize_nvec:'d resize_fn ->
?linsolv:('d, 'kind) session_linear_solver ->
('d, 'kind) tolerance ->
?restol:('d, 'kind) res_tolerance ->
float -> ('d, 'kind) Nvector.t -> float -> unit

Change the number of equations and unknowns between integrator steps. The call resize s ~resize_nvec:rfn ~linsolv:ls tol ~restol hscale ynew t0 has as arguments:

If a new linear solver is not specified, any existing linear solver in use is destroyed and reinitialized; settings must be reconfigured after resizing.

If the mass matrix residual tolerance was previously set to ResVtolerance and restol is not given, then it is reset to ResStolerance with the default value.

The tol argument is ignored in versions 2.6.1 and 2.6.2 since it may cause a segmentation error.

Modifying the solver (optional input functions)

Optional inputs for ARKode

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

Sets the integration tolerances.

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

Sets the residual tolerance.

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

Resets all optional input parameters to their default values. Neither the problem-defining functions nor the root-finding functions are changed.

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

Specifies the order of accuracy for the polynomial interpolant used for dense output. The interpolant is used both for solution output values and implicit method predictors.

val set_diagnostics : ('d, 'k) session -> Sundials.Logfile.t -> unit

Write step adaptivity and solver diagnostics to the given file.

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

Do not write step adaptivity or solver diagnostics of a file.

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_init_step : ('d, 'k) session -> float -> unit

Specifies the initial step size.

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

Disables time step adaptivity and fix the step size for all internal steps. Pass None to restore time step adaptivity. This function is not recommended since there is no assurance of the validity of the computed solutions. It is provided primarily for code-to-code verification testing. Use in conjunction with Arkode.set_min_step or Arkode.set_max_step.

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_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_err_test_fails : ('d, 'k) session -> int -> unit

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

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_optimal_params : ('d, 'k) session -> unit

Sets all adaptivity and solver parameters to ‘best guess’ values. This routine takes into account the integration method (ERK, DIRK, or ARK) and a given method order; it should only be called after these have been set.

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.

Optional inputs for IVP method selection

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

Enables both the implicit and explicit portions of a problem.

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

Disables the implicit portion of a problem.

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

Disables the explicit portion of a problem.

type rk_method = {
   stages : int; (*

Number of stages ($s$).

*)
   global_order : int; (*

Global order of accuracy ($q$).

*)
   global_embedded_order : int; (*

Global order of accuracy for the embedded RK method ($p$).

*)
}

Parameters for the RK method.

type rk_timescoefs = {
   stage_times : Sundials.RealArray.t; (*

Array (of length stages) of stage times ($c$).

*)
   coefficients : Sundials.RealArray.t; (*

Array (of length stages) of solution coefficients ($b$).

*)
   bembed : Sundials.RealArray.t option; (*

Optional array (of length stages) of embedding coefficients ($b2$).

*)
}

Coefficients for the RK method.

val set_ark_tables : ('d, 'k) session ->
rk_method ->
Sundials.RealArray.t ->
Sundials.RealArray.t -> rk_timescoefs -> rk_timescoefs -> unit

Specifies a customized Butcher table pair for the additive RK method. The call set_ark_tables rkm ai ae specifies:

In versions of Sundials prior to 2.7.0, the e parameter is ignored; only the i parameter is used.

If the i.bembed or the e.bembed field is None then the solver will run in fixed step mode and the step size must be set, or have been set, using either Arkode.set_fixed_step or Arkode.set_init_step. This feature is not available for Sundials versions prior to 2.7.0 (the Config.NotImplementedBySundialsVersion exception is raised).

val set_erk_table : ('d, 'k) session ->
rk_method -> Sundials.RealArray.t -> rk_timescoefs -> unit

Specifies a customized Butcher table pair for the explicit portion of the system. The call set_erk_table rkm ae specifies:

If the e.bembed field is None then the solver will run in fixed step mode and the step size must be set, or have been set, using either Arkode.set_fixed_step or Arkode.set_init_step. This feature is not available for Sundials versions prior to 2.7.0 (the Config.NotImplementedBySundialsVersion exception is raised).

val set_irk_table : ('d, 'k) session ->
rk_method -> Sundials.RealArray.t -> rk_timescoefs -> unit

Specifies a customized Butcher table pair for the implicit portion of the system. The call set_irk_table rkm ai specifies:

If the i.bembed field is None then the solver will run in fixed step mode and the step size must be set, or have been set, using either Arkode.set_fixed_step or Arkode.set_init_step. This feature is not available for Sundials versions prior to 2.7.0 (the Config.NotImplementedBySundialsVersion exception is raised).

type erk_table = 
| HeunEuler_2_1_2 (*

Default 2nd order explicit method (table 0).

*)
| BogackiShampine_4_2_3 (*

Default 3rd order explicit method (table 1).

*)
| ARK_4_2_3_Explicit (*

Explicit portion of default 3rd order additive method (table 2).

*)
| Zonneveld_5_3_4 (*

Default 4th order explicit method (table 3).

*)
| ARK_6_3_4_Explicit (*

Explicit portion of default 3rd order additive method (table 4).

*)
| SayfyAburub_6_3_4 (*

Butcher table number 5.

*)
| CashKarp_6_4_5 (*

Default 5th order explicit method (table 6).

*)
| Fehlberg_6_4_5 (*

Butcher table number 7.

*)
| DormandPrince_7_4_5 (*

Butcher table number 8.

*)
| ARK_8_4_5_Explicit (*

Explicit portion of default 5th order additive method (table 9).

*)
| Verner_8_5_6 (*

Default 6th order explicit method (table 10).

*)
| Fehlberg_13_7_8 (*

Default 8th order explicit method (table 11).

*)

Explicit Butcher tables

type irk_table = 
| SDIRK_2_1_2 (*

Default 2nd order implicit method (table 12).

*)
| Billington_3_2_3 (*

Butcher table number 13.

*)
| TRBDF2_3_2_3 (*

Butcher table number 14.

*)
| Kvaerno_4_2_3 (*

Butcher table number 15.

*)
| ARK_4_2_3_Implicit (*

Default 3rd order implicit method and the implicit portion of the default 3rd order additive method (table 16).

*)
| Cash_5_2_4 (*

Butcher table number 17.

*)
| Cash_5_3_4 (*

Butcher table number 18.

*)
| SDIRK_5_3_4 (*

Default 4th order implicit method (table 19).

*)
| Kvaerno_5_3_4 (*

Butcher table number 20.

*)
| ARK_6_3_4_Implicit (*

Implicit portion of the default 4th order additive method (table 21).

*)
| Kvaerno_7_4_5 (*

Butcher table number 22.

*)
| ARK_8_4_5_Implicit (*

Default 5th order method and the implicit portion of the default 5th order additive method (table 23).

*)

Implicit Butcher tables

type ark_table = 
| ARK_4_2_3 (*

3rd-order pair combining tables 2 and 15.

*)
| ARK_6_3_4 (*

4th-order pair combining tables 4 and 20.

*)
| ARK_8_4_5 (*

5th-order pair combining tables 9 and 22.

*)

Additive Butcher tables

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

Use specific built-in Butcher tables for an ImEx system.

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

Use specific built-in Butcher tables for an explicit integration of the problem.

The Fehlberg_13_7_8 method is not available prior to Sundials 2.7.0.

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

Use specific built-in Butcher tables for an implicit integration of the problem.

Optional inputs for time step adaptivity

type adaptivity_args = {
   h1 : float; (*

the current step size, $t_m - t_{m-1}$.

*)
   h2 : float; (*

the previous step size, $t_{m-1} - t_{m-2}$.

*)
   h3 : float; (*

the step size $t_{m-2} - t_{m-3}$.

*)
   e1 : float; (*

the error estimate from the current step, $m$.

*)
   e2 : float; (*

the error estimate from the previous step, $m-1$.

*)
   e3 : float; (*

the error estimate from the step $m-2$.

*)
   q : int; (*

the global order of accuracy for the integration method.

*)
   p : int; (*

the global order of accuracy for the embedding.

*)
}
type 'd adaptivity_fn = float -> 'd -> adaptivity_args -> float 

A function implementing a time step adaptivity algorithm that chooses an $h$ that satisfies the error tolerances. The call hnew = adapt_fn t y args has as arguments

This function should focus on accuracy-based time step estimation; for stability based time steps, Arkode.set_stability_fn should be used.

type adaptivity_params = {
   ks : (float * float * float) option;
   method_order : bool;
}

Parameters for the standard adaptivity algorithms. There are two:

type 'd adaptivity_method = 
| PIDcontroller of adaptivity_params (*

The default time adaptivity controller.

*)
| PIcontroller of adaptivity_params (*

Uses the two most recent step sizes.

*)
| Icontroller of adaptivity_params (*

Standard time adaptivity control algorithm.

*)
| ExplicitGustafsson of adaptivity_params (*

Primarily used with explicit RK methods.

*)
| ImplicitGustafsson of adaptivity_params (*

Primarily used with implicit RK methods.

*)
| ImExGustafsson of adaptivity_params (*

An ImEx version of the two preceding controllers.

*)
| AdaptivityFn of 'd adaptivity_fn (*

A custom time-step adaptivity function.

*)

Asymptotic error control algorithms.

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

Specifies the method and associated parameters used for time step adaptivity.

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

Specifies the fraction of the estimated explicitly stable step to use. Any non-positive argument resets to the default value (0.5).

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

Specifies the bias to apply to the error estimates within accuracy-based adaptivity strategies.

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

Specifies the step growth interval in which the step size will remain unchanged. In the call set_fixed_step_bounds s lb ub, lb specifies a lower bound on the window to leave the step size fixed and ub specifies an upper bound. Any interval not containing 1.0 resets to default values.

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

Specifies the maximum step size growth factor upon a convergence failure on a stage solve within a step. Any value outside the interval $(0, 1]$ resets to the default value.

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

Specifies the maximum step size growth factor upon multiple successive accuracy-based error failures in the solver. Any value outside the interval $(0, 1]$ resets to the default value.

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

Specifies the maximum allowed step size change following the very first integration step. Any value $\le 1$ resets to the default value.

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

Specifies the maximum growth of the step size between consecutive time steps. Any value $\le 1$ resets to the default value.

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

Specifies the safety factor to be applied to the accuracy-based estimated step. Any non-positive value resets to the default value.

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

Specifies the threshold for “multiple” successive error failures before the factor from Arkode.set_max_efail_growth is applied. Any non-positive value resets to the default value.

type 'd stability_fn = float -> 'd -> float 

A function that predicts the maximum stable step size for the explicit portions of an ImEx ODE system. The call hstab = stab_fn t y has as arguments

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

Sets a problem-dependent function to estimate a stable time step size for the explicit portion of the ODE system.

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

Clears the problem-dependent function that estimates a stable time step size for the explicit portion of the ODE system.

Optional inputs for implicit stage solves

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

Solve the implicit portion of the problem using the accelerated fixed-point solver instead of the modified Newton iteration. The integer argument gives the maximum dimension of the accelerated subspace (i.e., the number of vectors to store within the Anderson acceleration subspace).

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

Solve the implicit portion of the problem using the modified Newton solver. Overrides a previous call to Arkode.set_fixed_point.

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

Specifies that the implicit portion of the problem is linear. The flag indicates whether the Jacobian of $f_I(t,y)$, or, when using an iterative linear solver, the preconditioner is time-dependent (true) or not (false).

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

Specifies that the implicit portion of the problem is nonlinear.

type predictor_method = 
| TrivialPredictor (*

Piece-wise constant interpolant $p_0(\tau) = y_{n-1}$ %

*)
| MaximumOrderPredictor (*

An interpolant $p_q(t)$ of polynomial order up to $q=3$.

*)
| VariableOrderPredictor (*

Decrease the polynomial degree for later RK stages.

*)
| CutoffOrderPredictor (*

Maximum order for early RK stages and first-order for later ones.

*)
| BootstrapPredictor (*

Second-order predictor based only on the current step.

*)

Method choices for predicting implicit solutions.

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

Specifies the method for predicting implicit solutions.

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

Specifies the maximum number of nonlinear solver iterations permitted per RK stage at each 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_nonlin_crdown : ('d, 'k) session -> float -> unit

Specifies the constant used in estimating the nonlinear solver convergence rate.

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

Specifies the nonlinear correction threshold beyond which the iteration will be declared divergent.

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

Specifies a scaled step size ratio tolerance beyond which the linear solver setup routine will be signalled.

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

Specifies the frequency of calls to the linear solver setup routine. Positive values specify the number of time steps between setup calls, negative values force recomputation at each Newton step, and zero values reset to the default.

type 'd postprocess_step_fn = float -> 'd -> unit 

A function to process the results of each timestep solution. The arguments are

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

Set a post processing step function.

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

Clear the post processing step function.

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_exp_steps : ('d, 'k) session -> int

Returns the cumulative number of stability-limited steps taken by the solver.

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

Returns the cumulative number of accuracy-limited steps taken by the solver.

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

Returns the cumulative number of steps attempted by the solver.

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

Returns the number of calls to the right-hand side functions. In the call (nfe_evals, nfi_evals) = get_num_rhs_evals s,

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

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

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

Returns the integration step size taken on the last successful 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_current_butcher_tables : ('d, 'k) session ->
rk_method * Sundials.RealArray.t * Sundials.RealArray.t *
rk_timescoefs * rk_timescoefs

Returns the explicit and implicit Butcher tables in use by the solver. In the call (rkm, ai, ae, i, e) = get_current_butcher_tables s,

For versions of Sundials prior to 2.7.0, the i and e results are identical.

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.

*)
   exp_steps : int; (*

Cumulative number of stability-limited solver steps.

*)
   acc_steps : int; (*

Cumulative number of accuracy-limited solver steps.

*)
   step_attempts : int; (*

Cumulative number of steps attempted by the solver.

*)
   num_nfe_evals : int; (*

Number of calls to the explicit right-hand side function.

*)
   num_nfi_evals : int; (*

Number of calls to the implicit 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.

*)
   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.

Implicit solver optional output functions

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_nonlin_solv_iters : ('d, 'k) session -> int

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

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

Returns the number of nonlinear convergence failures that have occurred.

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 Arkode.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.

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 Arkode.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 Arkode.set_max_err_test_fails and Arkode.set_min_step.

exception ConvergenceFailure

Too many convergence test failures within a step or at the minimum step size. See Arkode.set_max_conv_fails and Arkode.set_min_step.

exception LinearInitFailure

Linear solver initialization failed.

exception LinearSetupFailure

Linear solver setup failed in an unrecoverable manner.

exception LinearSolveFailure

Linear solver solution failed in an unrecoverable manner.

exception MassInitFailure

Mass matrix solver initialization failed.

exception MassSetupFailure

Mass matrix solver setup failed in an unrecoverable manner.

exception MassSolveFailure

Mass matrix solver solution failed in an unrecoverable manner.

exception MassMultFailure

Mass matrix-vector multiplication failed.

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 PostprocStepFailure

The postprocess step function failed.

exception BadK

Raised by Arkode.get_dky for invalid order values.

exception BadT

Raised by Arkode.get_dky for invalid time values.