martinjrobins / diffsol Goto Github PK

• Explicit Runge-Kutta (ERK)
• #30
• fully implicit Runge-Kutta (FRK)

the NonLinearOp and LinearOp have a jacobian function which returns a matrix

fn jacobian(&self, x: &Self::V, t: Self::T) -> Self::M

This requires creating a new matrix everytime the jacobian is recalculated. Need an "inplace" version that can update an existing jacobian.

This will require additions to the Matrix trait to allow for mutation of existing matrices

related to #6 and #7

Add an optional stop time to step to halt the solve at this time point (e.g for time-based events, see #42 )

fn step(&mut self, stop_time: Option<Eqn::T>) -> Result<()>;

OdeEquations trait is duplicating lots of stuff from NonLinearOp and LinearOp, this trait should just return references to objects of these traits

at the moment OdeEquations has the following function for take state:

/// Take the current state of the solver, if it exists, returning it to the user. This is useful if you want to use this
/// state in another solver or problem. Note that this will unset the current problem and solver state, so you will need to call
/// `set_problem` again before calling `step` or `solve`.
fn take_state(&mut self) -> Option<OdeSolverState<Eqn::V>>;

However, when you set the state again with set_problem, there is a lot of unnecessary initialisation done, e.g. setting the initial timestep (see #34 ).

Option A:

I'm thinking of replacing take_step by a set_state, which would allow a user to set the internal OdeSolverState, while keeping the initialisation of other internal variables to a minimum

fn set_state(&mut self, new_state: &Eqn::V);

So user code would be:

let mut new_state = solver.state().clone();
new_state[0] += 1.0;
solver.set_state(&new_state);

Option B:

An alternative is to keep take_step, and add set_state that takes ownership of the vector, like so

fn set_state(&mut self, new_state: Eqn::V);

So then in user code you have:

let mut state = solver.take_state().unwrap();
state[0] += 1.0;
solver.set_state(state);

This allows you to mutate the state in-place so you don't have to have two versions of the state vector

atm the mass matrix is sorta optional in that it defaults to the identity matrix if not provided, but to be more explicit and allow algorithms to do efficiencies in the case of an identity mass matrix this should be truely optional, so the OdeEquations trait goes from

pub trait OdeEquations {
    /// returns the mass matrix `M` as a [LinearOp]
    fn mass(&self) -> &Rc<Self::Mass>;

to

pub trait OdeEquations {
    /// returns the mass matrix `M` as a [LinearOp]
    fn mass(&self) -> Option<&Rc<Self::Mass>>;

• OdeBuilder
• include latex in docs, see https://github.com/victe/rust-latex-doc-minimal-example
• everything else!

When set_problem is called for the solvers they calculate an initial step size. If the solve is a continuation of a previous solve this might not be required, so there should be an option to disable this

• #7
• #35
• sparse jacobian calculation (#6)
• implement sparse ode problem for testing

Add forward sensitivity calculation feature to all solvers

• Op traits need updating in order to supply the gradients wrt parameters
• perhaps use the methods used in sundials? https://sundials.readthedocs.io/en/latest/idas/Mathematics_link.html#forward-sensitivity-analysis

implement sparse jacobian coloring for efficient calculation of the jacobian, using https://github.com/JuliaDiff/SparseDiffTools.jl as a reference

mass matrix currently implemented via an implicit function $Mv$, for a large majority of cases $M=I$ or $M$ is constant, so I think we need to factor in efficiencies for these cases

• Should the solvers be aware if $M=I$ or not and act accordingly? Or should the efficiencies be implemented in OdeEquations?
• for the latter above, should OdeEquations provide a gemv style function for the implicit solvers to use?
• and an $M^{-1}$ function for the explicit solvers?

now that the rhs and mass types have been moved out of OdeEquations Trait, there is no need for these to be separate structs

Some other issues that might be relevant:

• rust-lang/rust#116914
• rust-lang/rust#114891
• rust-lang/rust#108826

at the moment you need to implement both Clone and Default for every type of matrix which leads to a lot of duplicated code. Should be able to clean this up

There needs to be a way for the user to specify events at which to terminate the solve. This could be another function in the OdeEquations trait that writes to a vector given the current state x and time t. The solver will terminate if any of the elements in that vector cross zero

atm the main call function for NonLinearOp is:

pub trait NonLinearOp: Op {
    /// Compute the operator at a given state and time.
    fn call_inplace(&self, x: &Self::V, t: Self::T, y: &mut Self::V);
    ...

To make this more flexible, we should use VectorView and VectorViewMut instead, ie:

pub trait NonLinearOp: Op {
    /// Compute the operator at a given state and time.
    fn call_inplace(&self, x: Self::VView, t: Self::T, y: Self::VViewMut);
    ...

This allows algorithms to call these functions with views rather than references to owned vectors. One usecase for this is to concatentate solves with varying parameters into a larger statevector

To allow Scale to be in the LhS this trait should be implemented.

add sundials solvers behind a feature gate

• #11
• add sundials ida solver with klu sparse linear solver (requires Chris00/sundials-sys#3)
• #15

Its a bit confusing atm who owns the solver state, I propose the API is refactored thus:

pub trait OdeSolverMethod<Eqn: OdeEquations> {

    fn problem(&self) -> Option<&OdeSolverProblem<Eqn>>;

    // solver now takes ownership of the state, which has been initialised by the user
    fn set_problem(&mut self, state: OdeSolverState<Eqn::M>, problem: &OdeSolverProblem<Eqn>);

    // solver can return a ref to the state for reading
    fn state(&self) -> Option<&OdeSolverState<Eqn>>

    // or you can take ownership of the state for editing, or initialising another problem/solver
    // note that `set_problem` must be called again after this, or subsequent `step` calls
    // will fail
    fn take_state(&mut self) -> Option<OdeSolverState<Eqn>>

    fn step(&mut self) -> Result<()>;
    fn interpolate(&self, t: Eqn::T) -> Eqn::V;
    fn solve(&mut self, problem: &OdeSolverProblem<Eqn>, t: Eqn::T) -> Result<Eqn::V> {
        let state = OdeSolverState::new(problem);
        self.set_problem(state, problem);
        while state.t <= t {
            self.step()?;
        }
        Ok(self.interpolate(t))
    }

    fn make_consistent_and_solve<RS: NonLinearSolver<FilterCallable<OdeRhs<Eqn>>>>(
        &mut self,
        problem: &OdeSolverProblem<Eqn>,
        t: Eqn::T,
        root_solver: &mut RS,
    ) -> Result<Eqn::V> {
        let state = OdeSolverState::new_consistent(problem, root_solver)?;
        self.set_problem(state, problem);
        while state.t <= t {
            self.step()?;
        }
        Ok(self.interpolate(t))
    }
}