Integrate the dassl DAE solver from DASSL.jl

[ivarne]

See:

• https://github.com/pwl/DASSL.jl
• https://groups.google.com/forum/m/#!topic/julia-users/egfvEo-4hP8

This might be a good excuse to continue the discussion on efficient ways to use different problem specifications in a consistent way.

@pwl

[pwl]

I have been working on DASSL API today after I discovered the issue #27. With relatively minor changes to the source code I was able to implement the iterative version of DASSL solver (see DASSL.jl/README.md) using coroutines and so far I haven't found any downsides of my implementation. I have also implemented some of the keyword arguments from the API.md, now you can use abstol, reltol, maxstep, minstep and initstep. I also modified the algorithm to prevent it from overshooting beyond tspan[end], in the iterative version I used a keyword tstop with default value tstop=Inf.

I have a problem with norm keyword argument as in DASSL the default norm is given by the formula

dassl_norm(v, wt)=norm(v./wt)/sqrt(length(v))

with the default weights depending on the current state of solution wt=reltol*abs(y).+abstol. The error estimate, given by error=dassl_norm(v,wt), is then kept to be less then one during integration. To my understanding this type of condition cannot be simply replaced by a norm with fixed weights. I know that this is not a typical approach to error estimation, but it is a part of DASSL algorithm and I am unsure whether I should modify it.

As for the output (see also #30 (comment) in issue #30), I chose to return a tuple (t,y,dy), adding the derivative dy to the standard output. The reason is that for implicit systems F(t,y,dy)=0 obtaining dy from t and y is non-trivial and may be computationally intensive. At the same time if the iterative solvers are introduced, it would be pretty easy to add an option to ignore the dy when saving all the steps in memory (although in my opinion, if you are constrained by memory dropping dy will bring you only a marginal advantage and you probably should dump your data on a disk anyway).

[ivarne]

Have you tested the overhead of using a task based approach, versus storing the state in a type object and using the normal iterator functions, versus pushing the results into arrays?

[pwl]

I have run some tests and there doesn't seem to be any overhead at all. (but I only compared the iterative and non-iterative versions)

[ChrisRackauckas]

Closing this because it's been completed by the organization layout and the common interface(s).