Merge pull request #24 from martinjrobins/docs

#19 more docs on the front-page and odeequations

Author: martinjrobins

src/lib.rs

@@ -2,7 +2,7 @@
 //!
 //! DiffSol is a library for solving differential equations. It provides a simple interface to solve ODEs and semi-explicit DAEs.
 //!
-//! ## Getting Started
+//! ## Solving ODEs
 //!
 //! To create a new problem, use the [OdeBuilder] struct. You can set the initial time, initial step size, relative tolerance, absolute tolerance, and parameters,
 //! or leave them at their default values. Then, call the [OdeBuilder::build_ode] method with the ODE equations, or the [OdeBuilder::build_ode_with_mass] method
@@ -11,8 +11,11 @@
 //! You will also need to choose a matrix type to use. DiffSol can use the [nalgebra](https://nalgebra.org) `DMatrix` type, or any other type that implements the
 //! [Matrix] trait. You can also use the [sundials](https://computation.llnl.gov/projects/sundials) library for the matrix and vector types (see [SundialsMatrix]).
 //!
-//! To solve the problem, you need to choose a solver. DiffSol provides a pure rust [Bdf] solver, or you can use the [SundialsIda] solver from the sundials library (requires the `sundials` feature).
-//! See the [OdeSolverMethod] trait for a more detailed description of the available methods on the solver.
+//! To solve the problem, you need to choose a solver. DiffSol provides the following solvers:
+//! - A pure rust Backwards Difference Formulae [Bdf] solver, suitable for stiff problems and singular mass matrices.
+//! - The IDA solver solver from the sundials library ([SundialsIda], requires the `sundials` feature).
+//!
+//! See the [OdeSolverMethod] trait for a more detailed description of the available methods on each solver.
 //!
 //! ```rust
 //! use diffsol::{OdeBuilder, Bdf, OdeSolverState, OdeSolverMethod};
@@ -45,6 +48,53 @@
 //! }
 //! let y = solver.interpolate(&state, t);
 //! ```
+//!
+//! ## DiffSL
+//!
+//! DiffSL is a domain-specific language for specifying differential equations <https://github.com/martinjrobins/diffsl>. It uses the LLVM compiler framwork
+//! to compile the equations to efficient machine code and uses the EnzymeAD library to compute the jacobian. You can use DiffSL with DiffSol by enabling one of the `diffsl-llvm*` features
+//! corresponding to the version of LLVM you have installed, and using the [OdeBuilder::build_diffsl] method.
+//!
+//! For more information on the DiffSL language, see the [DiffSL documentation](https://martinjrobins.github.io/diffsl/)
+//!
+//! ## Custom ODE problems
+//!
+//! The [OdeBuilder] struct can be used to create an ODE problem from a set of closures.
+//! If this is not suitable for your problem or you want more control over how your equations are implemented, you can also implement the [OdeEquations] trait manually.
+//!
+//! ## Nonlinear and linear solvers
+//!
+//! DiffSol provides generic nonlinear and linear solvers that are used internally by the ODE solver. You can use the solvers provided by DiffSol, or implement your own following the provided traits.
+//! The linear solver trait is [LinearSolver], and the nonlinear solver trait is [NonLinearSolver]. The [SolverProblem] struct is used to define the problem to solve.
+//!
+//! The provided linear solvers are:
+//! - [LU]: a direct solver that uses the LU decomposition implemented in the [nalgebra](https://nalgebra.org) 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:
+//! - [NewtonNonlinearSolver]: a nonlinear solver that uses the Newton method.
+//!
+//! ## Jacobian and Mass matrix calculation
+//!
+//! Via [OdeEquations], the user provides the action of the jacobian on a vector `J(x) v`. By default DiffSol uses this to generate a jacobian matrix for the ODE solver.
+//! Generally this requires `n` evaluations of the jacobian action for a system of size `n`, so it is often more efficient if the user can provide the jacobian matrix directly
+//! by implementing the [OdeEquations::jacobian_matrix] and the [OdeEquations::mass_matrix] (is applicable) functions.
+//!
+//! DiffSol also provides an experimental feature to calculate sparse jacobians more efficiently by automatically detecting the sparsity pattern of the jacobian and using
+//! colouring \[1\] to reduce the number of jacobian evaluations. You can enable this feature by enabling [OdeBuilder::use_coloring()] option when building the ODE problem.
+//!
+//! \[1\] Gebremedhin, A. H., Manne, F., & Pothen, A. (2005). What color is your Jacobian? Graph coloring for computing derivatives. SIAM review, 47(4), 629-705.
+//!
+//! ## Matrix and vector types
+//!
+//! 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.
+//! - [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:
+//! - For matrices: [Matrix], [MatrixView], [MatrixViewMut], [DenseMatrix], and [MatrixCommon].
+//! - For vectors: [Vector], [VectorIndex], [VectorView], [VectorViewMut], and [VectorCommon].
+//!
 
 #[cfg(feature = "diffsl-llvm10")]
 pub extern crate diffsl10_0 as diffsl;
@@ -100,7 +150,7 @@ pub use linear_solver::sundials::SundialsLinearSolver;
 #[cfg(feature = "sundials")]
 pub use ode_solver::sundials::SundialsIda;
 
-use matrix::{DenseMatrix, Matrix, MatrixViewMut};
+use matrix::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut};
 pub use nonlinear_solver::newton::NewtonNonlinearSolver;
 use nonlinear_solver::NonLinearSolver;
 pub use ode_solver::{
@@ -110,7 +160,7 @@ pub use ode_solver::{
 use op::NonLinearOp;
 use scalar::{IndexType, Scalar, Scale};
 use solver::SolverProblem;
-use vector::{Vector, VectorIndex, VectorRef, VectorView, VectorViewMut};
+use vector::{Vector, VectorCommon, VectorIndex, VectorRef, VectorView, VectorViewMut};
 
 pub use scalar::scale;
 

src/matrix/dense_serial.rs

@@ -3,7 +3,7 @@ use nalgebra::{DMatrix, DMatrixView, DMatrixViewMut, DVector, DVectorView, DVect
 
 use crate::{IndexType, Scalar};
 
-use super::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut};
+use crate::{DenseMatrix, Matrix, MatrixCommon, MatrixView, MatrixViewMut};
 
 impl<'a, T: Scalar> MatrixCommon for DMatrixViewMut<'a, T> {
     type V = DVector<T>;

src/ode_solver/diffsl.rs

@@ -220,6 +220,19 @@ mod tests {
 
     use super::DiffSl;
 
+    #[test]
+    fn diffsl_expontential_decay() {
+        // TODO: put this example into the docs once https://github.com/rust-lang/cargo/pull/13490 makes it to stable
+        let code = "
+            u { y = 1 }
+            F { -y }
+            out { y }
+        ";
+        let problem = OdeBuilder::new().build_diffsl(code).unwrap();
+        let mut solver = Bdf::default();
+        let _y = solver.solve(&problem, 1.0).unwrap();
+    }
+
     #[test]
     fn diffsl_logistic_growth() {
         let text = "

src/ode_solver/equations.rs

@@ -36,32 +36,47 @@ 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)
+/// $$
+///
+/// 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$.
 pub trait OdeEquations: Op {
-    /// This must be called first
+    /// The parameters of the ODE equations are assumed to be constant. This function sets the parameters to the given value before solving the ODE.
+    /// 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);
 
-    /// calculates $F(y)$ where $y$ is given in `y` and stores the result in `rhs_y`
+    /// calculates $F(t, y, p)$ where $y$ is given in `y` and stores the result in `rhs_y`. Note that the parameter vector $p$ is assumed to be
+    /// already provided via [Self::set_params()]
     fn rhs_inplace(&self, t: Self::T, y: &Self::V, rhs_y: &mut Self::V);
 
-    /// calculates $y = J(x)v$
+    /// calculates $y = J(x)v$, where $J(x)$ is the Jacobian matrix of the right-hand side function $F(t, y, p)$ at $y = x$. The result is stored in `y`.
     fn rhs_jac_inplace(&self, t: Self::T, x: &Self::V, v: &Self::V, y: &mut Self::V);
 
-    /// initializes `y` with the initial condition
+    /// returns the initial condition, i.e. $y(t_0, p)$
     fn init(&self, t: Self::T) -> Self::V;
 
+    /// calculates the right-hand side function $F(t, y, p)$ where $y$ is given in `y`. The result is allocated and returned.
+    /// The default implementation calls [Self::rhs_inplace()] and allocates a new vector for the result.
     fn rhs(&self, t: Self::T, y: &Self::V) -> Self::V {
         let mut rhs_y = Self::V::zeros(self.nstates());
         self.rhs_inplace(t, y, &mut rhs_y);
         rhs_y
     }
 
+    /// calculates $y = J(x)v$, where $J(x)$ is the Jacobian matrix of the right-hand side function $F(t, y, p)$ at $y = x$. The result is allocated and returned.
+    /// The default implementation calls [Self::rhs_jac_inplace()] and allocates a new vector for the result.
     fn jac_mul(&self, t: Self::T, x: &Self::V, v: &Self::V) -> Self::V {
         let mut rhs_jac_y = Self::V::zeros(self.nstates());
         self.rhs_jac_inplace(t, x, v, &mut rhs_jac_y);
         rhs_jac_y
     }
 
-    /// rhs jacobian matrix J(x), re-use jacobian calculation from NonLinearOp
+    /// calculate and return the jacobian matrix $J(x)$ of the right-hand side function $F(t, y, p)$ at $y = x$.
+    /// The default implementation calls [Self::rhs_jac_inplace()] and uses the jacobian calculation in [NonLinearOp].
     fn jacobian_matrix(&self, x: &Self::V, t: Self::T) -> Self::M {
         let rhs_inplace = |x: &Self::V, _p: &Self::V, t: Self::T, y_rhs: &mut Self::V| {
             self.rhs_inplace(t, x, y_rhs);
@@ -81,28 +96,39 @@ pub trait OdeEquations: Op {
         closure.jacobian(x, t)
     }
 
-    /// mass matrix action: y = M(t)
+    /// calculate the action of the mass matrix $M$ on the vector $v$ at time $t$, i,e. $y = M(t)v$.
+    /// The default implementation assumes that the mass matrix is the identity matrix and returns $y = v$.
     fn mass_inplace(&self, _t: Self::T, v: &Self::V, y: &mut Self::V) {
         // assume identity mass matrix
         y.copy_from(v);
     }
 
-    /// returns the indices of the algebraic state variables
+    /// For semi-explicit DAEs (with zeros on the diagonal of the mass matrix), this function
+    /// returns the indices of the algebraic state variables. This is used to determine which
+    /// components of the solution vector are algebraic, and therefore must be solved for
+    /// when calculating the initial condition to make sure that the algebraic constraints are satisfied.
+    /// The default implementation returns an empty vector, assuming that the mass matrix is the identity matrix.
     fn algebraic_indices(&self) -> <Self::V as Vector>::Index {
         // assume identity mass matrix
         <Self::V as Vector>::Index::zeros(0)
     }
 
+    /// calculate and return the mass matrix $M(t)$ at time $t$.
+    /// The default implementation assumes that the mass matrix is the identity matrix and returns the identity matrix.
     fn mass_matrix(&self, _t: Self::T) -> Self::M {
         // assume identity mass matrix
         Self::M::from_diagonal(&Self::V::from_element(self.nstates(), Self::T::one()))
     }
 
+    /// calculate and return the statistics of the ODE equation object (i.e. how many times the right-hand side function was evaluated, how many times the jacobian was multiplied, etc.)
+    /// The default implementation returns an empty statistics object.
     fn get_statistics(&self) -> OdeEquationsStatistics {
         OdeEquationsStatistics::new()
     }
 }
 
+/// This struct implements the ODE equation trait [OdeEquations] for a given right-hand side function, jacobian function, mass matrix function, and initial condition function.
+/// These functions are provided as closures, and the parameters are assumed to be constant.
 pub struct OdeSolverEquations<M, F, G, H, I>
 where
     M: Matrix,
@@ -179,7 +205,6 @@ where
     }
 }
 
-// impl Op
 impl<M, F, G, H, I> Op for OdeSolverEquations<M, F, G, H, I>
 where
     M: Matrix,
@@ -288,6 +313,9 @@ where
     }
 }
 
+/// This struct implements the ODE equation trait [OdeEquations] for a given right-hand side function, jacobian function, and initial condition function.
+/// These functions are provided as closures, and the parameters are assumed to be constant.
+/// The mass matrix is assumed to be the identity matrix.
 pub struct OdeSolverEquationsMassI<M, F, G, I>
 where
     M: Matrix,

src/solver/mod.rs

@@ -10,6 +10,7 @@ pub struct SolverStatistics {
     pub nmaxiter: IndexType,
 }
 
+/// A generic linear or nonlinear solver problem, containing the function to solve $f(t, y)$, the current time $t$, and the relative and absolute tolerances.
 pub struct SolverProblem<C: Op> {
     pub f: Rc<C>,
     pub t: C::T,
@@ -46,6 +47,8 @@ impl<C: Op> SolverProblem<C> {
 }
 
 impl<C: NonLinearOp> SolverProblem<C> {
+    /// Create a new solver problem from a nonlinear operator that solves for the linearised operator.
+    /// That is, if the original function is $f(t, y)$, this function creates a new problem $f'$ that solves $f' = J(x) v$, where $J(x)$ is the Jacobian of $f$ at $x$.
     pub fn linearise(&self, x: &C::V) -> SolverProblem<LinearisedOp<C>> {
         let linearised_f = Rc::new(LinearisedOp::new(self.f.clone(), x));
         SolverProblem::new_from_problem(linearised_f, self)

src/vector/faer_serial.rs

@@ -4,7 +4,7 @@ use faer::{unzipped, zipped, Col, ColMut, ColRef};
 
 use crate::{scalar::Scale, IndexType, Scalar};
 
-use super::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut};
+use crate::{Vector, VectorCommon, VectorIndex, VectorView, VectorViewMut};
 
 macro_rules! impl_op_for_faer_struct {
     ($struct:ident, $trait_name:ident, $func_name:ident) => {