Merge pull request #41 from martinjrobins/faer-doc

docs: faer containers, ode_equations and ops

src/lib.rs

@@ -70,6 +70,7 @@
 //!
 //! The provided linear solvers are:
 //! - [NalgebraLU]: a direct solver that uses the LU decomposition implemented in the [nalgebra](https://nalgebra.org) library.
+//! - [FaerLU]: a direct solver that uses the LU decomposition implemented in the [faer](https://github.com/sarah-ek/faer-rs) library.
 //! - [SundialsLinearSolver]: a linear solver that uses the [sundials](https://computation.llnl.gov/projects/sundials) library (requires the `sundials` feature).
 //!
 //! The provided nonlinear solvers are:
@@ -90,6 +91,7 @@
 //!
 //! When solving ODEs, you will need to choose a matrix and vector type to use. DiffSol uses the following types:
 //! - [nalgebra::DMatrix] and [nalgebra::DVector] from the [nalgebra](https://nalgebra.org) library.
+//! - [faer::Mat] and [faer::Col] from the [faer](https://github.com/sarah-ek/faer-rs) library.
 //! - [SundialsMatrix] and [SundialsVector] from the [sundials](https://computation.llnl.gov/projects/sundials) library (requires the `sundials` feature).
 //!
 //! If you wish to use your own matrix and vector types, you will need to implement the following traits:

src/ode_solver/bdf/mod.rs

@@ -372,11 +372,6 @@ where
     }
 
     fn step(&mut self) -> Result<()> {
-        // we will try and use the old jacobian unless convergence of newton iteration
-        // fails
-        // tells callable to update rhs jacobian if the jacobian is requested (by nonlinear solver)
-        // initialise step size and try to make the step,
-        // iterate, reducing step size until error is in bounds
         let mut d: Eqn::V;
         let mut safety: Eqn::T;
         let mut error_norm: Eqn::T;

src/ode_solver/equations.rs

@@ -37,11 +37,14 @@ impl OdeEquationsStatistics {
 /// this is the trait that defines the ODE equations of the form
 ///
 /// $$
-///  M \frac{dy}{dt} = F(t, y, p)
-///  y(t_0) = y_0(t_0, p)
+///  M \frac{dy}{dt} = F(t, y)
+///  y(t_0) = y_0(t_0)
 /// $$
 ///
-/// The ODE equations are defined by the right-hand side function $F(t, y, p)$, the initial condition $y_0(t_0, p)$, and the mass matrix $M$.
+/// The ODE equations are defined by:
+/// - the right-hand side function `F(t, y)`, which is given as a [NonLinearOp] using the `Rhs` associated type and [Self::rhs] function,
+/// - the mass matrix `M` which is given as a [LinearOp] using the `Mass` associated type and the [Self::mass] function,
+/// - the initial condition `y_0(t_0)`, which is given using the [Self::init] function.
 pub trait OdeEquations {
     type T: Scalar;
     type V: Vector<T = Self::T>;
@@ -53,12 +56,16 @@ pub trait OdeEquations {
     /// Note that `set_params` must always be called before calling any of the other functions in this trait.
     fn set_params(&mut self, p: Self::V);
 
+    /// returns the right-hand side function `F(t, y)` as a [NonLinearOp]
     fn rhs(&self) -> &Rc<Self::Rhs>;
+
+    /// returns the mass matrix `M` as a [LinearOp]
     fn mass(&self) -> &Rc<Self::Mass>;
 
-    /// returns the initial condition, i.e. $y(t_0, p)$
+    /// returns the initial condition, i.e. `y(t)`, where `t` is the initial time
     fn init(&self, t: Self::T) -> Self::V;
 
+    /// returns true if the mass matrix is constant over time
     fn is_mass_constant(&self) -> bool {
         true
     }

src/op/mod.rs

@@ -13,19 +13,29 @@ pub mod matrix;
 pub mod sdirk;
 pub mod unit;
 
+/// Op is a trait for operators that, given a paramter vector `p`, operates on an input vector `x` to produce an output vector `y`.
+/// It defines the number of states (i.e. length of `x`), the number of outputs (i.e. length of `y`), and number of parameters (i.e. length of `p`) of the operator.
+/// It also defines the type of the scalar, vector, and matrices used in the operator.
 pub trait Op {
     type T: Scalar;
     type V: Vector<T = Self::T>;
     type M: Matrix<T = Self::T, V = Self::V>;
+
+    /// Return the number of input states of the operator.
     fn nstates(&self) -> usize;
+
+    /// Return the number of outputs of the operator.
     fn nout(&self) -> usize;
+
+    /// Return the number of parameters of the operator.
     fn nparams(&self) -> usize;
 
     /// Return sparsity information for the jacobian or matrix (if available)
     fn sparsity(&self) -> Option<&<Self::M as Matrix>::Sparsity> {
         None
     }
 
+    /// Return statistics about the operator (e.g. how many times it was called, how many times the jacobian was computed, etc.)
     fn statistics(&self) -> OpStatistics {
         OpStatistics::default()
     }
@@ -70,11 +80,15 @@ pub trait NonLinearOp: Op {
 
     /// Compute the product of the Jacobian with a given vector.
     fn jac_mul_inplace(&self, x: &Self::V, t: Self::T, v: &Self::V, y: &mut Self::V);
+
+    /// Compute the operator at a given state and time, and return the result.
     fn call(&self, x: &Self::V, t: Self::T) -> Self::V {
         let mut y = Self::V::zeros(self.nout());
         self.call_inplace(x, t, &mut y);
         y
     }
+
+    /// Compute the product of the Jacobian with a given vector, and return the result.
     fn jac_mul(&self, x: &Self::V, t: Self::T, v: &Self::V) -> Self::V {
         let mut y = Self::V::zeros(self.nstates());
         self.jac_mul_inplace(x, t, v, &mut y);
@@ -87,6 +101,7 @@ pub trait NonLinearOp: Op {
         self._default_jacobian_inplace(x, t, y);
     }
 
+    /// Default implementation of the Jacobian computation.
     fn _default_jacobian_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::M) {
         let mut v = Self::V::zeros(self.nstates());
         let mut col = Self::V::zeros(self.nout());
@@ -98,6 +113,7 @@ pub trait NonLinearOp: Op {
         }
     }
 
+    /// Compute the Jacobian of the operator and return it.
     fn jacobian(&self, x: &Self::V, t: Self::T) -> Self::M {
         let n = self.nstates();
         let mut y = Self::M::new_from_sparsity(n, n, self.sparsity());
@@ -106,24 +122,31 @@ pub trait NonLinearOp: Op {
     }
 }
 
+/// LinearOp is a trait for linear operators (i.e. they only depend linearly on the input `x`). It extends the Op trait with methods for
+/// calling the operator via a GEMV operation (i.e. `y = t * A * x + beta * y`), and for computing the matrix representation of the operator.
 pub trait LinearOp: Op {
+    /// Compute the operator at a given state and time
     fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V) {
         let beta = Self::T::zero();
         self.gemv_inplace(x, t, beta, y);
     }
 
+    /// Computer the operator via a GEMV operation (i.e. `y = t * A * x + beta * y`)
     fn gemv_inplace(&self, x: &Self::V, t: Self::T, beta: Self::T, y: &mut Self::V);
 
+    /// Compute the matrix representation of the operator and return it.
     fn matrix(&self, t: Self::T) -> Self::M {
         let mut y = Self::M::new_from_sparsity(self.nstates(), self.nstates(), self.sparsity());
         self.matrix_inplace(t, &mut y);
         y
     }
 
+    /// Compute the matrix representation of the operator and store it in the matrix `y`.
     fn matrix_inplace(&self, t: Self::T, y: &mut Self::M) {
         self._default_matrix_inplace(t, y);
     }
 
+    /// Default implementation of the matrix computation.
     fn _default_matrix_inplace(&self, t: Self::T, y: &mut Self::M) {
         let mut v = Self::V::zeros(self.nstates());
         let mut col = Self::V::zeros(self.nout());