Open MartinOtter opened 5 years ago
There's tons of other things that can be done here as well. Lots of allocations to get rid of, overloads for user-Jacobians, etc. But I don't plan on working on it at all because this is somewhat of a dead-end development-wise. Essentially, we'd have to add every feature of OrdinaryDiffEq.jl here, or just add one function to OrdinaryDiffEq.jl and it'll get all of the features, so we're doing the later. So this code has just been kept pretty much in tact from its original creation (way before I started using Julia, this is "inherited code"), and sooner rather than later (a GSoC is slated for it) we'll have OrdinaryDiffEq.jl handle it, and this code will be deprecated. If you want to work on it a bit that's fine, but I don't plan to use of my own time improving this codebase myself.
[...] sooner rather than later (a GSoC is slated for it) we'll have OrdinaryDiffEq.jl handle it
Hm. Just to make sure that I understand correctly: DASSL.jl will be replaced by a new implementation within OrdinaryDiffEq.jl? So, it is not worth to put effort in improving DASSL.jl?
Indeed. OrdinaryDiffEq already does a few mass-matrix ODEs, and it just needs a small tweak to allow DAEProblems. Once that's done, someone just needs to make a perform_step!
overload for how DAE BDF works, and then it will have all of the events, GPU, Newton-Krylov, adaptivity, etc. support for free, which is why we do it that way. And we need the consistent initial condition finder from DASSL.
So I think that a better use of time would be to just find out why QNDF is not performing as well as we'd like and maybe implement a divided differences form as well. But for example, the whole code for QNDF is https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/bdf_perform_step.jl#L739-L878
which is much easier to maintain and give features to rather than having a library per method (we are at 300 methods, so...)
I'll open a new issue in OrdinaryDiffEq.jl to discuss the BDF differences. @YingboMa and I have been talking about it a lot internally but we should start writing things down.
File src/DASSL.jl line 686 has the function:
This (central) function should be improved:
In the for loop
f(y)
is called in every iteration, althoughy
does not change. This means that there aren-1
unnecessary calls of the right hand side function. The call off(y)
should be moved out of the for loop and the result stored in a vector.The statement
jac = Array{eltype(y)}(undef,n,n)
generates the Jacobian matrix whenever this function is called. So, when this function is called 100 times, then memory for 100 Jacobians are allocated. This should be improved: At the start of the simulation, memory for one Jacobian should be allocated that is then updated within this function.