OrdinaryDiffEq.jl icon indicating copy to clipboard operation
OrdinaryDiffEq.jl copied to clipboard

Implicit Radau5DAE

Open ufechner7 opened this issue 4 years ago • 8 comments

Please make it possible to use Radau5DAE for implicit problems. See: https://jmodelica.org/assimulo/DAE_Radau5DAE.html This would make it much easier to port existing code from Python to Julia.

For my problems it has a much better stability and performance than the sundials IDA solver.

In Assimulo both a pure Python version and a Fortran wrapper for Python are available which should make it easy to create a similar solution for Julia.

ufechner7 avatar May 13 '21 17:05 ufechner7

We have RadauIIA5 and in ODEInterfaceDiffEq we have a wrapper to the Fortran radau5 (and radau). But the original Fortran radau only does linearly implicit, i.e. allows singular mass matrices. It doesn't allow the fully implicit form. So jmodelica has either modified the solver (which, I don't quite see how to change the IRK here), or they must do an automatic conversion of the fully implicit DAE to a linearly-implicit system. @YingboMa have you looked at this?

ChrisRackauckas avatar May 14 '21 16:05 ChrisRackauckas

RADAU5 implicit Runge-Kutta method of order 5 (Radau IIA) for problems of the form My'=f(x,y) with possibly singular matrix M; with dense output (collocation solution). Concerning the linear algebra routines the user has the choice to link the program either with DC_DECSOL and DECSOL or with DC_LAPACK and LAPACK and LAPACKC From: http://www.unige.ch/~hairer/software.html

So it DOES allow a singular matrix M.

ufechner7 avatar May 14 '21 19:05 ufechner7

Yes sorry, that was a typo. Our RadauIIA5, radau, and radau5 methods all allow DAEs through a singular mass matrix.

ChrisRackauckas avatar May 14 '21 19:05 ChrisRackauckas

https://diffeq.sciml.ai/stable/solvers/dae_solve/ has all of the DAE solvers, mentioning fully implicit and mass matrix form.

ChrisRackauckas avatar May 14 '21 19:05 ChrisRackauckas

Any progress on this? What is missing?

ufechner7 avatar Jan 03 '22 03:01 ufechner7

It's missing someone to work on it. RadauIIA5 supports DAEs in mass matrix form and that can be used today. A separate form would require a PR.

ChrisRackauckas avatar Jan 03 '22 03:01 ChrisRackauckas

If I would like to implement it, which files would have to be modified?

ufechner7 avatar Jan 03 '22 09:01 ufechner7

It would just follow the standards of the OrdinaryDiffEq.jl repository. This is documented in the developer documentation http://devdocs.sciml.ai/latest/. http://devdocs.sciml.ai/latest/contributing/adding_algorithms/ describes how new algorithms are added to OrdinaryDiffEq.jl with an example on explicit RK methods. Some of the fine details of Radau and DAEs are not documented there, but that should get you started on what's important (defining a cache and a perform_step! function). I think the right way to do this might be to add a new DRadauIIA5 that uses the same cache and perform_step! dispatches as RadauIIA5 but hooks into the f calls in the right way with a if prob isa DAEProblem to swap over to the other f call etc.

ChrisRackauckas avatar Jan 05 '22 11:01 ChrisRackauckas