The DiffEq Internals

# The DiffEq Internals

The DiffEq solvers, OrdineryDiffEq, StochasticDiffEq, FiniteElementDiffEq, etc. all follow a similar scheme which leads to rapid development and high performance. This portion of the documentation explains how the algorithms are written.

## Developing New Solver Algorithms

The easiest way to get started would be to add new solver algorithms. This is a pretty simple task as there are tools which put you right into the "hot loop". For example, take a look at the ODE solver code. The mode `solve(::ODEProblem,::OrdinaryDiffEqAlgorithm)` is glue code to a bunch of solver algorithms. The algorithms which are coded in DifferentialEquations.jl can be found in src/integrators.jl. The actual step is denoted by the `perform_step!(integrator)` function. For example, take a look at the Midpoint method's implementation:

``````@inline function perform_step!(integrator,cache::MidpointConstantCache,f=integrator.f)
@unpack t,dt,uprev,u,k = integrator
halfdt = dt/2
k = integrator.fsalfirst
k = f(t+halfdt,uprev+halfdt*k)
u = uprev + dt*k
integrator.fsallast = f(u,p,t+dt) # For interpolation, then FSAL'd
integrator.k = integrator.fsalfirst
integrator.k = integrator.fsallast
@pack integrator = t,dt,u
end``````

The available items are all unloaded from the `integrator` in the first line. `fsalfirst` inherits the value of `fsallast` on the last line. The algorithm is written in this form so that way the derivative of both endpints is defined, allowing the vector `integrator.k` to determine a Hermite interpolation polynomial (in general, the `k` values for each algorithm form the interpolating polynomial). Other than that, the algorithm is as basic as it gets for the Midpoint method, making sure to set `fsallast` at the last line. The results are then packaged back into the integrator to mutate the state.

Note the name `ConstantCache`. In OrdinaryDiffEq.jl, the `Algorithm` types are used for holding type information about which solver to choose, and its the "cache" types which then hold the internal caches and state. `ConstantCache` types are for not inplace calculations `f(t,u)` and `Cache` types (like `MidpointCache`) include all of the internal arrays. Their constructor is specified in the `cache.jl` file. The cache (and the first `fsalfirst`) is initialized in the `initialize!` function next to the cache's `perform_step!` function.

The main inner loop can be summarized by the `solve!` command:

``````@inbounds while !isempty(integrator.opts.tstops)
while integrator.tdir*integrator.t < integrator.tdir*top(integrator.opts.tstops)
@ode_exit_conditions
perform_step!(integrator,integrator.cache)
loopfooter!(integrator)
if isempty(integrator.opts.tstops)
break
end
end
handle_tstop!(integrator)
end
postamble!(integrator)``````

The algorithm runs until `tstop` is empty. It hits the `loopheader!` in order to accept/reject the previous step and choose a new `dt`. This is done at the top so that way the iterator interface happens "mid-step" instead of "post-step". In here the algorithm is also chosen for the switching algorithms.

Then the `perform_step!` is called. (The exit conditions throw an error if necessary.) Afterwards, the `loopfooter!` is used to calculate new timesteps, save, and apply the callbacks. If a value of `tstops` is hit, the algorithm breaks out of the inner-most loop to save the value.

Adding algorithms to the other problems is very similar, just in a different pacakge.

## Extras

If the method is a FSAL method then it needs to be set via `isfsal` and `fsalfirst` should be defined before the loop, with `fsallast` what's pushed up to `fsalfirst` upon a successful step. See `:DP5` for an example.

If tests fail due to units (i.e. Unitful), don't worry. I would be willing to fix that up. To do so, you have to make sure you keep separate your `rateType`s and your `uType`s since the rates from `f` will have units of `u` but divided by a unit of time. If you simply try to write these into `u`, the units part will fail (normally you have to multiply by a \$dt\$).