DifferentialEquations.jl inplace implicit ODE problem

General Usage
1

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, i.e., A(du,u,p,t)=0, since our mass matrices are not just identity matrices, and can be fairly complicated or even nonlinear.

Finally, what we would provide are a residual(res,du,u,p,t) and jacobian(J,du,u,p,t) functions.

As a result, I think that the way to go is to consider an IIP problem in DifferentialEquations.jl. But it seems that only the DAE constructor allows such syntax.

However, for implicit ODE systems I cannot use an ODE algorithm, e.g., Backward-Euler, even though it is acceptable for our problems when the mass matrix is non-singular (thus not a DAE). And it seems we must provide an initial value for du that does not have mathematical sense in our problems.

Summarizing, is there a way to create an ODE solver that takes residual!(res,du,u,p,t) and jacobian!(J,du,u,p,t) inplace functions, and allows one to use an ODE algorithm?

I have written the following silly linear example with what I would like to have

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{true}(residual;jac=jacobian)
prob_iip = ODEProblem{true}(f_iip,u0,tspan)
sol_iip = solve(prob_iip, ImplicitEuler(), reltol=1e-8, abstol=1e-8)

but it does not work, because it seems the library is looking for an ODE-style inplace version of residual, i.e., residual(du,u,p,t), based on the error:

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

If you can write it as a matrix, you could use the state-dependent mass matrix part of the interface, and it can be more efficient.

function update_func(A,u,p,t)
    A[1,1] = cos(t)
    A[2,1] = sin(t)*u[1]
    A[3,1] = t^2
    A[1,2] = cos(t)*sin(t)
    A[2,2] = cos(t)^2 + u[3]
    A[3,2] = sin(t)*u[2]
    A[1,3] = sin(t)
    A[2,3] = t^2
    A[3,3] = t*cos(t) + 1
end
dependent_M = DiffEqArrayOperator(ones(3,3),update_func=update_func)
prob = ODEProblem(ODEFunction(f, mass_matrix=mm_A), u0, tspan))

It’s just used in the consistent initialization algorithms. You can set it to all zero if you want. The interface for DAEProblem should probably make that du0 optional in the future, but for now it’s required but sometimes ignored by solvers/initialization algorithms.

No. You can use a state-dependent mass matrix with an algorithm for ODEs, or you can use a DAEProblem with “implicit ODE” DAEs.

3

Gridap.jl is nice! I would like to use it with PseudoArcLengthContinuation.jl. I will see how you expose your jacobian in this thread

2 Likes
4

Thanks for the quick reply, I will consider this option.

In any case, following

https://docs.sciml.ai/latest/features/performance_overloads/#

I am creating a DAE problem

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)
sol_iip = solve(prob_iip, DImplicitEuler(), reltol=1e-8, abstol=1e-8)

but it returns an error,

type ODEIntegrator has no field duprev

which is weird, since I would expect a DAEIntegrator being created.

It would be great if you could tell me what I am doing wrong.

5

It’s us that’s doing something wrong. The new DAE solvers don’t support user Jacobians quite yet:

https://github.com/SciML/OrdinaryDiffEq.jl/issues/1109

Bare with us: these new integrators are missing the polish right now, so I’d still recommend using IDA() for production for the time being (though these are improving daily!).

1 Like