Explicit Jacobians for DAEs

[ChrisRackauckas]

MWE:

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

The issue is threefold.

1. The Jacobian references duprev which doesn't exist. This should be fsalfirst
2. fsalfirst isn't defined in the first step, so that will be undefined until the first Jacobian is calculated, which gives a causality issue...
3. jac for a DAE is actually the W-matrix, so that needs to move to the calc_W! section.

[oxinabox]

I am just adding the error message that the above gives, so I can find this again with search

julia> sol_iip = solve(prob_iip, DImplicitEuler(), reltol=1e-8, abstol=1e-8)

ERROR: type Nothing has no field α
Stacktrace:
  [1] setproperty!(x::Any, f::Symbol, v::Any)
    @ Base ./Base.jl:39 [inlined]
  [2] calc_W!(W::Any, integrator::Any, nlsolver::Union{…}, cache::Any, dtgamma::Any, repeat_step::Any, W_transform::Any, newJW::Any)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/derivative_utils.jl:668 [inlined]
  [3] update_W!(nlsolver::OrdinaryDiffEq.AbstractNLSolver, integrator::SciMLBase.DEIntegrator{…}, cache::Any, dtgamma::Any, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/derivative_utils.jl:796 [inlined]
  [4] update_W!(nlsolver::OrdinaryDiffEq.AbstractNLSolver, integrator::SciMLBase.DEIntegrator{…}, cache::Any, dtgamma::Any, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/derivative_utils.jl:795 [inlined]
  [5] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.DImplicitEulerCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/nlsolve/nlsolve.jl:25
  [6] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.DImplicitEulerCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/perform_step/dae_perform_step.jl:67
  [7] perform_step!(integrator::Any, cache::OrdinaryDiffEq.DImplicitEulerCache, repeat_step::Any)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/perform_step/dae_perform_step.jl:60 [inlined]
  [8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:520
  [9] __solve(::DAEProblem{…}, ::DImplicitEuler{…}; kwargs::Base.Pairs{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/87BwC/src/solve.jl:6
 [10] solve_call(_prob::DAEProblem{…}, args::DImplicitEuler{…}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:511 [inlined]
 [11] solve_up(prob::DAEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::DImplicitEuler{…}; kwargs::Base.Pairs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:972
 [12] solve_up(prob::SciMLBase.AbstractDEProblem, sensealg::Any, u0::Any, p::Any, args::Vararg{Any})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:945 [inlined]
 [13] solve(prob::DAEProblem{…}, args::DImplicitEuler{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::Base.Pairs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HoOGI/src/solve.jl:882
 [14] top-level scope
    @ REPL[10]:1
Some type information was truncated. Use `show(err)` to see complete types.

We run into this problem often in the thing I am working on.
And the solution is to use IDA from Sundials.jl instead of any DAESolver from OrdinaryDiffEq.jl