Module Ida

A session with the IDA solver.

An example session with Ida (ida_skel.ml):

open Sundials

(* 1. Define a residual function. *)
let resf _t v v' r =
  r.{0} <- v.{2}  -. v'.{0};
  r.{1} <- v.{3}  -. v'.{1};
  r.{2} <- v'.{2} -. v.{4} *. v.{0};
  r.{3} <- v'.{3} -. v.{4} *. v.{1} +. 9.81;
  r.{4} <- v.{2}*.v.{2} +. v'.{2}*.v.{0} +. v.{3}*.v'.{1} +. v'.{3}*.v.{1}

(* 2. Optionally define a root function. *)
let g _t v _v' gout = gout.{0} <- v.{0} -. v.{1} *. 0.5774

(* 3. Set vector of initial values.
      The length of this vector determines the problem size. *)
let vd = RealArray.of_list [ 0.9848; 0.1736; 0.0; 0.0; 0.0 ]
let v  = Nvector_serial.wrap vd
let v' = Nvector_serial.make 5 0.0

(* 4. Create and initialize a solver session.
      This will initialize a specific linear solver and the root-finding
      mechanism, if necessary. *)
let m = Matrix.dense 5
let s = Ida.(init (SStolerances (1e-9, 1e-9))
                  ~lsolver:Dls.(solver (dense v m))
                  resf ~roots:(1, g) 0.0 v v');;

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

(* 6. Correct initial values *)
let vids = Nvector_serial.make 5 Ida.VarId.differential;;
(Nvector_serial.unwrap vids).{4} <- Ida.VarId.algebraic;

Ida.set_suppress_alg s ~varid:vids true;
Ida.calc_ic_ya_yd' s ~y:v ~y':v' 0.1;;

(* 7. 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\t% .10e\t% .10e\t% .10e\n"
                t vd.{0} vd.{1} vd.{2} vd.{3} vd.{4};
  match r with
  | Ida.Success -> go (Ida.solve_normal s (t +. 0.05) v v')
  | Ida.RootsFound -> begin
         vd.{2} <- -. 0.5 *. vd.{2};
         vd.{3} <- -. 0.5 *. vd.{3};
         Ida.reinit s t v v';
         Ida.calc_ic_ya_yd' ~y:v ~y':v' s ~varid:vids (t +. 0.05);
          go (t, Ida.Success)
      end
  | Ida.StopTimeReached -> ();;

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

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

Alias for sessions based on serial nvectors.

Linear solvers used by Ida.

Alias for linear solvers that are restricted to serial nvectors.

Workspaces with two temporary vectors.

Workspaces with three temporary vectors.

Arguments common to Jacobian callback functions.

Residual functions that define a DAE problem. They are passed four * arguments:

• t, the value of the independent variable, i.e., the simulation time,
• y, the vector of dependent-variable values, i.e., $y(t)$,
• y', the vector of dependent-variable derivatives, i.e., $\dot{y} = \frac{\mathrm{d}y}{\mathrm{d}t}$, and,
• r a vector for storing the residual value, $F(t, y, \dot{y})$.

Within the function, raising a Sundials.RecoverableFailure exception indicates a recoverable error. Any other exception is treated as an unrecoverable error.

Create an IDA-specific linear solver from a generic matrix embedded solver.

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

Tolerance specifications.

A default relative tolerance of 1.0e-4 and absolute tolerance of 1.0e-8.

Called by the solver to calculate the values of root functions. These ‘zero-crossings’ are used to detect significant events. The function is passed four arguments:

• t, the value of the independent variable, i.e., the simulation time,
• y, the vector of dependent-variable values, i.e., $y(t)$,
• y', the vector of dependent-variable derivatives, i.e., $\dot{y} = \frac{\mathrm{d}y}{\mathrm{d}t}$, and,
• gout, a vector for storing the value of $g(t, y, \dot{y})$.

Creates and initializes a session with the solver. The call

init linsolv tol ~nlsolver ~nlsresfn ~lsolver f ~varid:varid ~roots:(nroots, g) t0 y0 y'0

has as arguments:

• tol, the integration tolerances,
• nlsolver, the solver to use to calculate integration steps and initial conditions,
• nlsresfn, alternative residual function to use in nonlinear system function evaluations,
• lsolver, used by nlsolvers based on Newton interation,
• f, the DAE residual function,
• varid, optionally classifies variables as algebraic or differential,
• nroots, the number of root functions,
• g, the root function ((nroots, g) defaults to Ida.no_roots),
• t0, the initial value of the independent variable,
• y0, a vector of initial values that also determines the number of equations, and,
• y'0, the initial values for $\dot{y} = \frac{\mathrm{d}y}{\mathrm{d}t}$.

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

If an nlsolver is not specified, then the Newton module is used by default. The nlsolver must be of type RootFind, otherwise an Ida.IllInput exception is raised.

The alternative residual 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.

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

Class components of the state vector as either algebraic or differential. These classifications are required by Ida.calc_ic_ya_yd' and Ida.set_suppress_alg. See also Ida.VarId.

Indicates whether or not to ignore algebraic variables in the local error test. When ignoring algebraic variables (true), a varid vector must be specified either in the call or by a prior call to Ida.init or Ida.set_id. Suppressing local error tests for algebraic variables is discouraged for DAE systems of index 1 and encouraged for systems of index 2 or more.

Computes the initial state vector for certain index-one problems. All components of $y$ are computed, using $\dot{y}$, to satisfy the constraint $F(t_0, y_0, \dot{y}_0) = 0$. If given, the ~y vector is filled with the corrected values. The last argument is the first vale of $t$ at which a solution will be requested. A Ida.reinit is required before calling this function after Ida.solve_normal or Ida.solve_one_step.

Computes the algebraic components of the initial state and derivative vectors for certain index-one problems. The elements of $y$ and $\dot{y}$ marked as algebraic are computed, using $\dot{y}$, to satisfy the constraint $F(t_0, y_0, \dot{y}_0) = 0$. The variable ids must be given in ~varid or by a prior call to Ida.init or Ida.set_id. If given, the ~y and ~y' vectors are filled with the corrected values. The last argument is the first vale of $t$ at which a solution will be requested. A Ida.reinit is required before calling this function after Ida.solve_normal or Ida.solve_one_step.

Specifies the positive constant in the nonlinear convergence test of the initial condition calculation.

Specifies the maximum number of steps taken in attempting to reach a given output time in the initial condition calculation.

Specifies the maximum number of approximate Jacobian or preconditioner evaluations allowed when the Newton iteration appears to be slowly converging.

Specifies the maximum number of Newton iterations allowed in any one attempt to calculate initial conditions.

Specifies the maximum number of linesearch backtracks allowed in any Newton iteration, when solving the initial conditions calculation problem.

Enables (true) or disables (false) the linesearch algorithm in the initial condition calculation.

Specifies a positive lower bound on the Newton step in the initial condition calculation.

Returns the number of backtrack operations during Ida.calc_ic_ya_yd' or Ida.calc_ic_y.

Values returned by the step functions. Failures are indicated by exceptions.

Integrates a DAE system over an interval. The call tret, r = solve_normal s tout yout y'out has as arguments

• s, a solver session,
• tout the next time at which the solution is desired,
• yout, a vector to store the computed solution, and,
• y'out, a vector to store the computed derivative.

It returns tret, the time reached by the solver, which will be equal to tout if no errors occur, and, r, a Ida.solver_result.

Like Ida.solve_normal but returns after one internal solver step.

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 Ida.get_current_time and $h_u$ denotes Ida.get_last_step,— and $0 \leq \mathtt{k} \leq q_u$—where $q_u$ denotes Ida.get_last_order.

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

Reinitializes the solver with new parameters and state values. The values of the independent variable, i.e., the simulation time, the state variables, and the derivatives must be given. If the argument ~linsolv is not given, the current linear solver remains unchanged. The argument ~roots works similarly; pass Ida.no_roots to disable root finding. The ~resfn argument gives the possibility to change the residual function. Similarly for ~nlsresfn and the alternative residual function for nonlinear system function iterations.

Set the integration tolerances.

Configure the default error handler to write messages to a file. By default it writes to Logfile.stderr.

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

Restores the default error handling function.

Specifies the maximum order of the linear multistep method.

Specifies the maximum number of steps taken in attempting to reach a given output time.

Specifies the initial step size.

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

Limits the value of the independent variable t when solving. By default no stop time is imposed.

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

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

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

Specifies the safety factor used in the nonlinear convergence test.

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

Disables constraint checking.

Returns the sizes of the real and integer workspaces.

Returns the cumulative number of internal solver steps.

Returns the number of calls to the residual function.

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

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

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

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

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

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

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

Returns the the current internal time reached by the solver.

Returns a suggested factor by which the user's tolerances should be scaled when too much accuracy has been requested for some internal step.

Returns the solution error weights at the current time.

Returns the vector of estimated local errors.

Summaries of integrator statistics.

Returns the integrator statistics as a group.

Prints the integrator statistics on the given channel.

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

Returns the cumulative number of nonlinear convergence failures.

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

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.

Like Ida.set_root_direction but specifies a single direction for all root functions.

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

Returns the number of root functions.

Fills an array showing which functions were found to have a root.

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

Returns the scalar $c_j$ , which is proportional to the inverse of the step size.

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

Returns the current $\dot{y}$ vector. This vector provides direct access to the data within the integrator.

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

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 yyn, ypn, and res vectors with, respectively: $\mathit{yyn} = y_{\mathit{pred}} + y_{\mathit{cor}}$ , $\mathit{ypn} = \dot{y}_{\mathit{pred}} + \alpha\dot{y}_{\mathit{cor}}$ , and $\mathit{res} = F(t_n, y_n, \dot{y}_n)$ where $y_{\mathit{cor}}$ is the first argument of the nonlinear solver system function. Within a custom linear solver, then the vectors yyn, ypn, and res are only current after an evaluation of the nonlinear system function.

Computes the current $y$ vector from a correction vector.

Computes the current $\dot{y}$ vector from a correction vector.

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.

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

The requested accuracy could not be satisfied.

Too many error test failures within a step. See Ida.set_max_err_test_fails.

Too many convergence test failures within a step, or Newton convergence failed. See Ida.set_max_conv_fails.

Linear solver initialization failed.

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.

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.

The nonlinear solver failed in a general way.

Nonlinear solver initialization failed.

Nonlinear solver setup failed in an unrecoverable manner.

Nonlinear solver setup failed in a recoverable manner.

The residual function failed in an unrecoverable manner.

The residual function had a recoverable error when first called.

Too many convergence test failures, or unable to estimate the initial step size, due to repeated recoverable errors in the residual function.

The rootfinding function failed.

No solution satisfying the inequality constraints could be found.

Linesearch could not find a solution with a step larger than steptol in weighted RMS norm.

A recoverable error occurred in a callback but no recovery was possible.

A component of the error weight vector, either for the input value or a corrected value, is zero.

Raised by Ida.get_dky for invalid order values.

Raised by Ida.get_dky for invalid time values.

Variable ids are required but not set. See Ida.set_id.

A fused vector operation failed.