Inspiration from PETSc ODE/DAE solvers for API and algorithms

[mauro3]

This is to propose to support the API PETSc uses for its ODE/DAE solvers.
I think the PETSc framework is the nicest around and the most modern (afaik).
(Other issues related are #20, #18, JuliaLang/julia#75)

The PETSc ODE/DAE solvers have been implemented over the last few years with modern developments in mind (they are citing papers from 2003, 2009), which is more than can be said for most other suites of ODE solvers which carry decades old ballast around. It supports, explicit, implicit and IMEX methods for ODEs and DAEs all within the same API. The API is based around formulating the problem like

F(t, u, u') = G(t, u)

For problems using implicit methods this is used as F(t, u, u') =0 for
explicit ones to u'=G(t, u) (i.e. F=u'), IMEX mixes the two at once.

Best have a look at the PETSc manual which gives a nice concise intro (Chapter 6):
http://www.mcs.anl.gov/petsc/petsc-current/docs/manual.pdf

PETSc has a liberal licence so code could probably be translated
from PETSc into Julia. The problem with that approach is that
one needs to understand the whole PETSc framework, otherwise
the code looks quite obfuscated.

Further, some cool solvers PETSc has:

• General linear methods:
  http://www.mcs.anl.gov/petsc/petsc-3.4/docs/manualpages/TS/TSGL.html#TSGL
  These methods contain Runge-Kutta and multistep
  schemes (e.g. BDF) as special cases. As far as I can tell with
  my limited ODE/DAE solver knowledge, these are the best
  currently available.
• IMEX: implicit explicit methods
  http://www.mcs.anl.gov/petsc/petsc-3.4/docs/manualpages/TS/TSARKIMEX.html
  "These methods are intended for problems with well-separated
  time scales, especially when a slow scale is strongly nonlinear
  such that it is expensive to solve with a fully implicit
  method."

Disclaimer: I have not actually used the PETSc solvers, just read
the docs and wished I could use them in Matlab...

[fsaporito]

I read over the project list in the GSOC page on the julia site (julialang.org/gsoc/2014/)
There is this project proposal:

"Project: PETSc integration for scalable technical computing

PETSc is a widely used framework of data structures and computational routines suitable for massively scaling scientific computations. Many of these algorithms are also ideally suited for big data applications such as computing principal components of very large sparse matrices and solving complicated forecasting models with distributed methods for solving partial differential equations.

This project proposal is to develop a new Julia package to interface with PETsc, thus allowing users access to state of the art scalable algorithms for optimization, eigenproblem solvers, finite element mesh computations, and hyperbolic partial differential equation solvers. The more mathematically oriented student may choose to study the performance of these various algorithms as compared to other libraries and naïve implementations. Alternatively, students may also be interested in working on the LLVM BlueGene port for deploying Julia with PetSc integration in an actual supercomputing environment.

Expected Results: Wrapper package for the PETSc suite of libraries."

It's a good idea to recreate the API in native Julia when, if someone does this project, the PETCs framework can be called directly from Julia?
In my opinion it's worthwhile to do so, for two main reasons, but others may think otherwise:

1. Speed
   While surely the PETSc implementation is very optmised, at such level that it would be difficult to reach within a reasonable time, calling over an external library over and over is computationally heavy, and for Julia would means to grow slower.

2. Glue Language
   If we start from creating interfaces to other libraries instead of coding ours implementations, and the process continues, in a few time Julia will slowly became a glue languages for linking external libraries in a common framework. While I don't believe that a language like that is a total garbage, since it could have a few uses, I don't think that this is the goal for the Julia language that anybody had in mind.

What is your opinion on this?

[jiahao]

You're quite right to be concerned at the possibility of overhead because of calling an external library. In fact, benchmarking the PETSc interface vs one or a few particular algorithms reimplemented in native Julia could make for a very nice component of a summer project. The usual counterargument is that computing the functions and its derivatives will probably be expensive anyway, but it's an assertion that can be tested.

[stevengj]

The main reason to implement this in pure Julia is to support ODE integration with arbitrary types; basically, you want to support integrating ODE functions over any normed vector space, i.e. any type that supports at least +, -, * by scalar, and a norm. (Not to mention support for different scalar types, e.g. BigFloat.)

[ivarne]

Another reason to implement in Julia, is that many more people will be able to dive into the source code for the solvers they are using if it is written in Julia.

This topic for this issue is anyway to draw inspiration from the API's for the PETSc solvers, into our own API discussions, not reimplementing everything in Julia. If someone are going to create wrappers for PETSc in Julia, it would be nice if they could make it effortlessly as a new backend for the solvers in ODE.jl

[mauro3]

Here a quick comparison of the PETSc and the matlab API for stiff solvers, e.g. Matlab's ode15s and ode15i.

ode15s can be applied to problems of the form M(t,u) u'=G(t,u) where M is called the mass matrix. Thus it can be applied to nearly as general problems as above PETCs setup. The mass matrix needs to be set in the options structure with Mass, MStateDependence, MvPattern and MassSingular. This is a bit clunky, considering that the mass matrix is a central part of the equation and not an option.

ode15i can be applied to problems of the form F(t, u, u') = 0, i.e. as the LHS of the PETSc setup. Here there is no need to set the mass matrix. (however, ode15i seems to be a lot less efficient that ode15s)

The PETSc interface unifies the two matlab interfaces quite nicely. Also solvers which only support the F or the G term would only implement that. Also in PETSc the difference between F and G is that for F a Jacobian is needed (either provided of by finite difference) whereas for G non Jacobian is needed as an explicit method is applied to that.

(PS: anther reason to re-implement some of the PETSc-solvers in Julia is that PETSc is a big dependency which involves a lot of compilation. Too much to ask of someone just wanting to solve a ODE.)

[tomasaschan]

If we decided to generalize the API to the form PETSc uses, i.e. ode(F, G, y0, tspan; kwargs...), the solvers we use today would only be able to use one of F and G to define the system. If I understand the PETSc API correctly, that's the case in their solvers as well.

What is a good way to handle the case when the user calls a solver with both F and G specified, but the solver cannot handle that type of systems? How do we detect it, and how do we let the user know they're doing something that won't work?

I think the idea of the API is really nice, and perhaps it's possible to create a wrapper function which, much like e.g. the \ operator, inspects the input arguments and chooses a solver based on the properties of the problem. The user would simply call ode and let the library (try to) work it out. If one needs more direct control over what solver method is used (e.g. specific performance/accuracy requirements), the specific solver methods are of course also available.

(This comment somehow disappeared when I posted it about an hour ago...)

[stevengj]

Regarding API, I think we should prescribe both:

ts,ys = odeXX(G, y0, tspan; kws...)  # solve y' = G(t,y)
ts,ys = odeXX(F, G, y0, tspan; kws...) # solve F(t,y,y') = G(t,y)

Not all solvers will support the second form, which is fine; they will just throw a MethodError.

All solvers which support the second form should also define the first form. Most simply, they may of course define odeXX(G, y0, tspan; kws...) = odeXX((t,y,yp) -> yp, G, y0, tspan; kws...), but I expect that they may want to have special-case optimizations for F = y'.

[mauro3]

This sounds good and would also be compatible with the API already defined.

[stevengj]

On second thought, although we should definitely have a simple interface odeXX(G, y0, tspan), I'm not sure that odeXX(F, G, y0, tspan) is sufficient for the complex interface, because you potentially want to specify other pieces of information, such as the Jacobian of F (and for efficiency you really want to allow a single function to compute and return both F and the Jacobian dF/dy' at the same time), whether F or G is linear, and so on.

(PETSc has you create a problem "object" and provides various methods to modify its properties. I suppose we could have something in a similar spirit but more Julian.)

[mauro3]

Matlab also sets the Jacobian or its sparsity structure in the ode-options.

We could have a keyword argument, say FJac&GJac, to either (1) pass a Jacobian function to, or (2) to provide a sparsity pattern (for finite differencing), or (3) to indicate to use automatic differentiation, or (4) to indicate that F or G return both F or Gand its Jacobian. If we wanted to use dispatch instead, then maybe a trick like #31 may work (until the Function-type contains its signature...if ever).

[acroy]

We have actually introduced a keyword jacobian in #14 for ode4s (non-adaptive Rosenbrock). Currently, we only check if jacobian is a function and otherwise use a (crude?) finite-difference approximation. We could extend that to support other options as well.

I like the idea of having a "problem type" and using dispatch to select the most appropriate solver, but IMO we should also support a low-level interface like odeXX(F, G, y0, tspan) and odeXX(F, y0, tspan).

[stevengj]

@acroy, a problem type is a lower level interface than F, G, etcetera. It will allow one to avoid anonymous functions, allow the Jacobian to be computed at the same time as the function, and so on. But agree that we should still have the odeXX(F, G, y0, tspan) and odeXX(F, y0, tspan) high-level simple interfaces.

[tomasaschan]

Rather than throw a MethodError when calling e.g. ode45(F, G, y0, tspan) (since ode45 is an explicit solver and therefore can't handle any system with F != y') I think we should only define this level of API with methods that actually make sense for the solvers in question. Implicit solvers for systems like F(t,y,y') benefit from specifying F(t,y,y') (and can, possibly, solve F(t,y,y') = G(t,y) simply by rewriting the system as F - G = 0 if necessary, but maybe with a performance hit), but for explicit solvers that can't handle that type of systems I think it makes more sense to just not define a method with that signature. ode45 et al solve y' = G(t,y), so then you only specify G.

I'm unsure about the problem specific types, and whether there is a place for them in this package. At the moment, ODE.jl focuses on very general methods, usable on a few very wide classes of systems. There could absolutely be a case for implementing specific solvers for specific problems, if there are such solvers that are common enough (and good enough) to be worth the effort, but then it might make more sense to just define a solver method with a name that indicates what problem we're solving, and that instead of general system functions take specific parameters for the problem.

What I was talking about in case of a general wrapper function, was something that would look like this (in pseudocode):

function odesol(F, G, y0, tspan)
    # if F == y' and G is nonstiff
        return our_normal_solver(G, y0, tspan)
    # else if F == y' and G is stiff
        return our_stiff_solver(G, y0, tspan) # with a kwarg for the jacobian, if specified
    # else if F != y' and G == 0
        return our_implicit_solver(F, y0, tspan) # with jacobian as above
   # etc for various other cases
end

our_normal_solver et al are reasonable choices for the type of problems, e.g. ode45 for non-stiff, explicit problems, and for whatever properties we can extract about F and G we can choose an appropriate solver. That way, the user will first try to use odesol to solve the problem, and if that works (s)he doesn't even have to care what method was used. If more detailed control is needed, ode45 and the lot are of course still available, but to "just do it" you don't have to to even choose a solver.

[stevengj]

@tlycken, only defining the API for methods that make sense is what throws a MethodError. You don't need to explicitly throw a MethodError, you just don't define the method.

The problem with your general wrapper is that there is no easy way to implement a test like F == y'.

[pwl]

Speaking from my experience I think it would be nice to support four types of equations with increasing generality:

1. y'=F(t,y)
2. M(t)*y'=F(t,y)
3. M(y,t)*y'=F(t,y)
4. F(t,y,y')=0

Optionally, the second case could be split into M=M(t) and M=const.

I think these four types strike a balance between providing enough information on the equation and being easy enough to handle them through the API. Unfortunately, this complicates the process of passing the user defined jacobian as for each type of equation we have to provide a different version of jacobian. For example in the first case we should provide dF/dy, but in the third case we have to additionally provide dM/dy. In the fourth case we would need dF/dy and dF/dy'. The simplest solution that comes into my mind would be to accept F as the required argument (despite F from 1. being different from F from 2.) and use keyword arguments for M and dM/dy, dF/dy, dF/dy'. Neither of the jacobians should be required for the code to work and they should be independent of each other, e.g. user should be able to pass dM/dy without passing dF/dy. Additionally to distinguish the fourth case from the first case user would set a flag DAE=true or something alike. Alternatively, instead of the keyword DAE we could use type="DAE", which would also give us the option to explicitly specify other types of equations type="MassODE" for 2. or type="MassYODE" for 3.

The downside of this approach is that it would need a lot of code inside the functions ode23 etc. to distinguish between all different cases. The remedy could be to implement some kind of dispatcher which would be called like ode(..., method="ode23"), in the likeness of NDSolve from Mathematica, which would filter invalid argument combinations so that ode23 could just assume it got correct arguments.

Please let me know what do you think about these ideas.

[stevengj]

@pwl, I don't think that anything but the simplest interfaces, possibly just case (1), should have the user passing a function directly. For everything else, the user should just pass an opaque ODEProblem subtype with appropriate methods defined on it to compute M, Jacobians etcetera, possibly in-place. We can then use method_exists to determine on this object to determine what methods the user has provided.

[pwl]

@stevengj I think this is an even better option for practical purposes: you could seamlessly switch between various methods and apply them to the same ODEProblem. Although it definitely adds complexity to the use cases.

What about returning additional data from the solver? Data such as the number of function/jacobian evaluations, local error estimates, steps accepted vs. steps rejected etc. These results can be very solver-dependent, but on the other hand they might be useful to compare the performance of different solvers.

[stevengj]

I'm not sure what "this" you are referring to.

For benchmarking, all you really want is to count either wall-clock time or number of function/Jacobian evaluations. Both of these can be implemented by the user (e.g. by adding a counter into their function) without modifying the solver.

[pwl]

Sorry for being ambiguous. By "this" I meant your proposal to enclose all the equation-related data into one structure. You are also correct about the counters, some of them can be definitely removed from the implementation of the solvers.

[pwl]

What about outputting derivatives as well as the solution? For explicit systems like y'=F(t,y) it is trivial to obtain the derivative of the solution by just evaluating F, but with implicit systems F(t,y,y')=0 it is not that easy. In both cases the values of y' are computed in each step anyway, so it makes sense not to waste this information and return a tuple (t,y,dy) as a result.

[ivarne]

Currently doc/api.md does not describe the API for implicit systems. For large problems, it might be problematic to hold on to more memory for a longer time, so we might want to somehow make it a choice. Maybe this can be a property of the ODEProblem as suggested in #31 and #33 (comment). Eg. you could define return_derivative(a::YourProblem) = true, to keep and return the derivative

[jedbrown]

Note that when integrating F(t,y,y') = G(t,y), you end up needing to solve with J[a] = dF/dy + a*dF/dy' (where a is defined by the method and step size). It is desirable for your Jacobian evaluation function to return this matrix directly in order to save memory and because J might be implemented matrix-free or using some other structured representation (possibly needed for preconditioning). IDA in Sundials has a similar interface (where the user builds the preconditioner for a given shift a instead of providing pieces and having the implementation put them together). I designed the PETSc interface the way I did because of a number of applications where special representations of matrices provide integer-factor to orders of magnitude benefit in terms of memory and time-to-solution in the implicit solver.

[delta-pi-1701]

https://github.com/hplgit/odespy