gabrielgellner / DiffEQ.jl

Staging Area for a Differential Solver suite in Julia (DEAD)
Other
0 stars 0 forks source link

Convert to internal in place evaluation of the RHS #13

Open gabrielgellner opened 9 years ago

gabrielgellner commented 9 years ago

So currently we have implemented a type system that is similar to how the fortran and C codes pass around memory for the solvers, currently we call this an ODESystem.

Looking at https://github.com/JuliaLang/ODE.jl/pull/33 and https://github.com/JuliaLang/ODE.jl/issues/30 a related issue is how to think about how the callback function is used internally. They key is to have the ability to change the function to work inplace instead of having:

dydt = f(t, y)

to have instead:

f!(t, y, dydt)

That way we don't incur the memory allocation on each call. The way that the above suggests to to have a special type for the function and then make a high level wrapper that takes a regular dydt = f(t, y) and place it into a dummy version that does the call inplace. This way the user has the option to make their own efficient version if they choose. This is similar to how I think the Optim.jl package uses DifferentiableFunction and TwiceDifferentiableFunction. Looking at the way the Optim.jl uses this to also allow for efficient calculation of both the function and the jacobian in one function call will be good to think about for when we try to implement stiff solvers.

So in summary it is worth looking into having as part of the ODESystem type some kind of ODEFunction or some such type to wrap the callback. Then have all internal calls use the inplace version, dummy or not.

In the issue 30 discussion there is a more subtle issue of how to model the form of the ODE in general citing the package PETSc which uses the model

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

for the general case. This seems interesting, but I have so little experience with this kind of issue I am not sure how to deal with it.