Creating a ODESolver from DifferentialEquations.jl

[santiagobadia]

Hi @fverdugo

I have been taking a look at DifferentialEquations.jl, in order to create ODESolvers as wrappers of this library solvers.

In our current version of the GridapTimeStepper package, I think that we can do it.

For instance, looking at the DAE solvers, what one requires is a residual, which is a method that we have.

https://docs.sciml.ai/stable/tutorials/dae_example/
https://docs.sciml.ai/stable/types/dae_types/

On the other hand, you can also provide a Jacobian, with a coefficient that comes from the solver, as we have. So, we could easily provide this method too.

https://docs.sciml.ai/stable/features/performance_overloads/

We can also provide our linear solvers (it says nonlinear too, but I cannot find an example)

https://docs.sciml.ai/stable/features/linear_nonlinear/#Linear-Solvers:-linsolve-Specification-1

We could create a solver from DifferentialEquation at every time step, from t_n to t_n+dt. I cannot see how to use the methods in the library as iterators, so they cannot be used for the whole period. Furthermore, it seems they are storing the unknown at all time steps (it can be tuned), which is not what we want.

It is unclear to me how to create a wrapper for their solvers such that we can easily provide all the options of the library in our creational methods for the wrapper.

What do you think?

[santiagobadia]

I am trying to combine the ODE solvers in DifferentialEquations.jl with the PDE solvers in the Gridap.jl library.

It is mandatory for our applications to use the inplace version of the DifferentialEquations.jl package. This way, we can pre-allocate the required arrays in our Gridap solvers.

On the other hand, even when considering ODE systems, we want to use an implicit ODE expression of the solver, since our mass matrices are not just identity matrices, and can be fairly complicated or even nonlinear.

As a result, I think that the way to go is to consider an IIP problem in DifferentialEquations.jl.

I have written the following silly linear example

tspan = (0.0,1.0)
u0 = zeros(1)
u0[1] = 1/2

function residual(res,du,u,p,t)
  res[1] = - 1.01u[1] + du[1]
end

function jacobian(jac,du,u,p,gamma,t)
  jac[1,1] = gamma - 1.01
  nothing
end

f_iip = ODEFunction{false}(residual;jac=jacobian)
prob_iip = ODEProblem{false}(f_iip,u0,tspan)
sol_iip = solve(prob_inplace, ImplicitEuler(), reltol=1e-8, abstol=1e-8)

but it does not work, because it seems the library is looking for an out-of-place version of jacobian, based on the error:

MethodError: no method matching jacobian(::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64)
Closest candidates are:
  jacobian(::Any, ::Any, ::Any, !Matched::Any, !Matched::Any, !Matched::Any) at /DifEqWrappers.jl:24
...

In fact, even without considering the Jacobian, the same problem appears with the residual when doing this:

prob_iip = ODEProblem{false}(residual,u0,tspan)
sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)

Is it possible to do what I want in DifferentialEquations.jl? If Yes, What am I doing wrong? I could not find an example in the documentation covering this (probably my fault).

[santiagobadia]

Instead, when using the DAE machinery

f_iip = DAEFunction{false}(residual;jac=jacobian)
prob_iip = DAEProblem{false}(f_iip,u0,tspan)
sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)

it returns an error saying that ImplicitEuler is an ODE solver, not a DAE solver

[fverdugo]

We could create a solver from DifferentialEquation at every time step, from t_n to t_n+dt. I cannot see how to use the methods in the library as iterators, so they cannot be used for the whole period. Furthermore, it seems they are storing the unknown at all time steps (it can be tuned), which is not what we want.

Yes, we can always implement the function solve_step! with the DifferentialEquations machinery and use our lazy GenericODESolution to represent all time steps in a lazy way. By the way, we only need to create the DifferentialEquations solver once for all time steps. EDIT: I am not sure if this is true since the inital condition and the time span are inside the ODESolver. We can always create a new solver at each time step as you suggested.

It is unclear to me how to create a wrapper for their solvers such that we can easily provide all the options of the library in our creational methods for the wrapper.

I need to have a look at the DifferentialEquations library to see their API before to be able to answer this question. If they have a "solver object". A flexible way is that the user builds this object directly in the driver using the DifferentialEquations API and passes it to our wrapper.

[fverdugo]

Is it possible to do what I want in DifferentialEquations.jl? If Yes, What am I doing wrong? I could not find an example in the documentation covering this (probably my fault).

It is possible to use in-place functions (at least for the residual) for systems of ODEs. See here:
https://docs.sciml.ai/latest/tutorials/ode_example/#Example-2:-Solving-Systems-of-Equations-1

In the link, they write ODEProblem instead of ODEProblem{false}. I don't know if this is the problem.

[santiagobadia]

The issue is a little bit old, please take a look at this post

https://discourse.julialang.org/t/differentialequations-jl-inplace-implicit-ode-problem/37495

It must be ODEProblem{true} where true means inplace but the tutorial does not apply to us, because there is no mass matrix. But you can consider a non-trivial mass matrix, as replied in the post. I knew this option but still, it is not all we need.

But still, we cannot provide a jacobian so I consider this is not enough. And it assumes a linear mass matrix, which is not always the case. So, up to now, I don't know how to make it work for the most general case.

I like the DAE definition in DifferentialEquations.jl, it really matches our needs, but when using this constructor combined with a ODE algorithm, e.g., ImplicitEuler(), it return an error, because we are trying to solve a DAE with an ODE algorithm. However, this is not always the case, when jacobian_t part is non-singular...

[ChrisRackauckas]

Hey, just found this. I can help you out here.

I cannot see how to use the methods in the library as iterators, so they cannot be used for the whole period

https://docs.sciml.ai/latest/basics/integrator/

Furthermore, it seems they are storing the unknown at all time steps (it can be tuned), which is not what we want.

https://docs.sciml.ai/latest/basics/common_solver_opts/#Output-Control-1 save_everystep=false,save_start=false saves nothing except the final timepoint, so you can save as much or as little as you want.

it says nonlinear too, but I cannot find an example

You probably never want to supply that, which is why we haven't updated support for that in awhile. Most standard nonlinear solver schemes are highly highly inefficient for time-stepping problems because you always want to do some kind of fancy quasi-Newton with Jacobian reuse. We'll re-add support here, but the built-in ones of quasi-Newton, Anderson acceleration, and fixed point are probably what you want anyways (and they do the proper reuse stuff).

It must be ODEProblem{true} where true means inplace but the tutorial does not apply to us, because there is no mass matrix. But you can consider a non-trivial mass matrix, as replied in the post. I knew this option but still, it is not all we need.

Yes, you probably always want to do ODEProblem{true}. I'm not sure why you were setting ODEProblem{false}, or where you even found that.

But still, we cannot provide a jacobian so I consider this is not enough. And it assumes a linear mass matrix, which is not always the case. So, up to now, I don't know how to make it work for the most general case.

Why not set the Jacobian? There's an example in the documentation: https://docs.sciml.ai/latest/features/performance_overloads/#Declaring-Explicit-Jacobians-for-DAEs-1

But still, we cannot provide a jacobian so I consider this is not enough. And it assumes a linear mass matrix, which is not always the case. So, up to now, I don't know how to make it work for the most general case.

It doesn't assume a linear mass matrix. It mentions in the documentation which solvers can handle a state-dependent mass matrix: M(u,p,t)u' = f(u,p,t), and ImplicitEuler is one of the ones in the list.

https://docs.sciml.ai/latest/solvers/dae_solve/#OrdinaryDiffEq.jl-(Mass-Matrix)-1

There's some examples of using it here:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/test/interface/mass_matrix_tests.jl#L36-L63

But you can consider a non-trivial mass matrix, as replied in the post. I knew this option but still, it is not all we need.

I'm curious it isn't, because it can express anything that can be expressed in residual form.

---

Path Forward

I think the easiest path forward is to first target the DAEProblem, defining the residual form and the Jacobian, and using IDA() . In this style, your example is just:

using Sundials
tspan = (0.0,1.0)
u0 = zeros(1)
u0[1] = 1/2

function residual(res,du,u,p,t)
  res[1] = - 1.01u[1] + du[1]
end

function jacobian(jac,du,u,p,gamma,t)
  jac[1,1] = gamma - 1.01
  nothing
end

f_iip = DAEFunction{true}(residual;jac=jacobian)
prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
sol_iip = solve(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)

or in the integrator interface:

using Sundials
tspan = (0.0,1.0)
u0 = zeros(1)
u0[1] = 1/2

function residual(res,du,u,p,t)
  res[1] = - 1.01u[1] + du[1]
end

function jacobian(jac,du,u,p,gamma,t)
  jac[1,1] = gamma - 1.01
  nothing
end

f_iip = DAEFunction{true}(residual;jac=jacobian)
prob_iip = DAEProblem{true}(f_iip,u0,u0,tspan,differential_vars=[true])
integ = init(prob_iip, IDA(), reltol=1e-8, abstol=1e-8)
step!(integ)
step!(integ)

# Show of using integrators as iterators
using Base.Iterators
for i in take(integ,100)
      @show integ.u
end

The fully implicit form will get you working, and we can then discuss the mass matrix forms to open up new solver possibilities, but this can at least get a version 1 working. The downside with using the Sundials solver here is it doesn't have the full flexibility yet, but we can get it to what you need quite easily. Do you have an example of the residual form I could test on?

[santiagobadia]

Dear @ChrisRackauckas

Thanks a lot for you reply! And sorry for the huge delay answering.

I have followed your suggestions, and created a driver in a branch in which I have used Sundials with Gridap

https://github.com/gridap/GridapODEs.jl/blob/diffeqwrapper/test/TransientFEsTests/DiffEqsWrapperTests.jl

If you have some time, you can take a look and let me know.

At the end of this file, I have listed the main issues I see in the current solution, attached below:

• Issue 1: How do I initialize the jacobian (sparse CSR matrix)? Dense matrices cannot be used in FE codes.

• Issue 2: No control over the (non)linear solver, we would like to be able to provide certainly linear and possibly nonlinear solvers. Let us assume that we just want to consider a fixed point algorithm and we consider an implicit time integration of a nonlinear PDE. Our solvers are efficient since they re-use cache among nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization) and for the transient case in our library, we can also reuse all this between time steps. Could we attain something like this using DifferentialEquations/Sundials?

• Issue 3: Iterator-like version as in GridapODEs. HOWEVER, is this computation lazy? I don't think so, so we need to store
  all time steps, e.g., for computing time-dependent functionals (drag or lift in CFD, etc), which is going to incur in a lot of memory consumption.

It would be great if you can suggest any idea or recommendation.

Cheers

--Santi

[ChrisRackauckas]

Issue 1: How do I initialize the jacobian (sparse CSR matrix)? Dense matrices cannot be used in FE codes.

f_iip = DAEFunction{true}(residual;jac=jacobian,jac_prototype=mysparsematrix)

See https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ for more tutorials on that stuff.

Issue 2: No control over the (non)linear solver, we would like to be able to provide certainly linear and possibly nonlinear solvers. Let us assume that we just want to consider a fixed point algorithm and we consider an implicit time integration of a nonlinear PDE.

https://docs.sciml.ai/latest/tutorials/advanced_ode_example/ shows swapping out linear solvers. https://docs.sciml.ai/latest/features/linear_nonlinear/ is all of the extra details.

Our solvers are efficient since they re-use cache among nonlinear steps (e.g., symbolic factorization, arrays for numerical factorization) and for the transient case in our library, we can also reuse all this between time steps. Could we attain something like this using DifferentialEquations/Sundials?

All of the DiffEq and Sundials solvers already do that.

Issue 3: Iterator-like version as in GridapODEs. HOWEVER, is this computation lazy? I don't think so, so we need to store all time steps, e.g., for computing time-dependent functionals (drag or lift in CFD, etc), which is going to incur in a lot of memory consumption.

Yes, it's lazy. It doesn't store anything if you do save_everystep=false,save_start=false.

[santiagobadia]

I have worked a little bit on the driver that makes use of Sundials in this file

I have found a couple of issues 1) related to jac_prototype (when I pass my own sparse matrix I get an internal error). If I don't pass the Jacobian prototype, 2) the nonlinear solver does not converge if the dimension of the ODE system is > 1, while the problem I am solving is linear (?). I guess I am doing something wrong.

To be investigated.

[kanav99]

Hi @santiagobadia, I am from SciML. I was having a look, do you have an example of the following issues? Any code snippets which I can change to reproduce the error?

[santiagobadia]

Hi @kanav99

Please take a look at this test

https://github.com/gridap/GridapODEs.jl/blob/diffeqwrapper/test/TransientFEsTests/SimpleDiffEqsTest.jl

When I pass jac_prototype it gets stuck.

[kanav99]

Yes I was having a look at that file. What I wanted to know was what you pass as jac_prototype that it gets stuck. If I could get a script that I can run to reproduce the error, that would be very helpful

[santiagobadia]

I pass a SparseArrays.SparseMatrixCSC{Float64,Int64}, which is J in the code.

[kanav99]

I was trying DImplicitEuler instead of IDA, I am not sure if IDA is compatible with Sparse Matrices (cc: @ChrisRackauckas ). There where some changes needed for OrdinaryDiffEq (SciML/OrdinaryDiffEq.jl#1231) but we need the jacobian function to have 5 (J, du, u, p, t) instead of 6 arguments (gamma is extra)

[santiagobadia]

Let us assume that I have a nonlinear problem A(t,u,u') = 0, and I can provide ∂A/∂u and ∂A/∂(u') or ∂A/∂u + γ*∂A/∂(u'), where γ represents d(Δₜ u)/du, with Δₜ being the discretisation of the time derivative in the time integrator. Without γ, I am not sure how to combine the partial derivative of the operator A wrt u and its time derivative u'. Should I assume a particular value for γ based on the solver I pick? E.g., assume that I am going to use Implicit Euler and take γ=1. It is even worse that that, because γ should depend on the time step size, so I cannot see how to build the Jacobian if γ is extra.