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
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
The main inner loop can be summarized by the
@inbounds while !isempty(integrator.opts.tstops) while integrator.tdir*integrator.t < integrator.tdir*top(integrator.opts.tstops) loopheader!(integrator) @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.
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.
If the method is a FSAL method then it needs to be set via
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
rateTypes and your
uTypes 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$).