Attempt at ForwardDiff Jacobian preparation when explicit Jacobian is provided yields error

[SouthEndMusic]

I isolated this M(N)WE:

using DifferentialEquations
using LinearAlgebra

function f!(du, u, p, t)
    du .= u
end

N = 10
jac_prototype = Matrix{Float64}(I, N, N)
function jac(J, u, p, t)
    J .= 0
    for i = 1:N
        J[i, i] = 1
    end
end

u0 = ones(N)
timespan = (0, 10)

RHS = ODEFunction(f!; jac_prototype, jac)
prob = ODEProblem(RHS, u0, timespan)
solve(
    prob,
    QNDF(;
        nlsolve = NonlinearSolveAlg(
            OrdinaryDiffEq.OrdinaryDiffEqNonlinearSolve.NewtonRaphson(;),
        ),
    ),
)

Stack trace:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearFunction{…}, Float64}, Float64, 10})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it. 

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::Float32)
   @ Base float.jl:341
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{NonlinearFunction{…}, Float64}, Float64, 10})        
    @ Base .\number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 10}, i::Int64)
    @ Base .\array.jl:987
  [3] macro expansion
    @ C:\Users\konin_bt\.julia\packages\FastBroadcast\wfdTr\src\FastBroadcast.jl:163 [inlined]
  [4] macro expansion
    @ .\simdloop.jl:77 [inlined]
  [5] __fast_materialize!
    @ C:\Users\konin_bt\.julia\packages\FastBroadcast\wfdTr\src\FastBroadcast.jl:162 [inlined]
  [6] _fast_materialize!
    @ C:\Users\konin_bt\.julia\packages\FastBroadcast\wfdTr\src\FastBroadcast.jl:196 [inlined]
  [7] fast_materialize!
    @ C:\Users\konin_bt\.julia\packages\FastBroadcast\wfdTr\src\FastBroadcast.jl:277 [inlined]
  [8] compute_ustep!
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\newton.jl:296 [inlined]
  [9] _compute_rhs!(tmp::Vector{…}, ztmp::Vector{…}, ustep::Vector{…}, γ::Rational{…}, α::Float64, tstep::Float64, k::Vector{…}, invγdt::Float64, method::OrdinaryDiffEqCore.MethodType, p::SciMLBase.NullParameters, dt::Float64, f::ODEFunction{…}, z::Vector{…})
    @ OrdinaryDiffEqNonlinearSolve C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\newton.jl:373
 [10] (::OrdinaryDiffEqNonlinearSolve.var"#6#9")(ztmp::Vector{…}, z::Vector{…}, p::Tuple{…})
    @ OrdinaryDiffEqNonlinearSolve C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\utils.jl:227
 [11] NonlinearFunction
    @ C:\Users\konin_bt\.julia\packages\SciMLBase\2TYas\src\scimlfunctions.jl:2471 [inlined]
 [12] FixTail
    @ C:\Users\konin_bt\.julia\packages\DifferentiationInterface\7eD1K\src\utils\context.jl:169 [inlined]
 [13] vector_mode_dual_eval!
    @ C:\Users\konin_bt\.julia\packages\ForwardDiff\UBbGT\src\apiutils.jl:31 [inlined]
 [14] vector_mode_jacobian(f!::DifferentiationInterface.FixTail{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff C:\Users\konin_bt\.julia\packages\ForwardDiff\UBbGT\src\jacobian.jl:139
 [15] jacobian
    @ C:\Users\konin_bt\.julia\packages\ForwardDiff\UBbGT\src\jacobian.jl:40 [inlined]
 [16] jacobian(f!::NonlinearFunction{…}, y::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{…}, backend::AutoForwardDiff{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceForwardDiffExt C:\Users\konin_bt\.julia\packages\DifferentiationInterface\7eD1K\ext\DifferentiationInterfaceForwardDiffExt\twoarg.jl:470
 [17] construct_jacobian_cache(prob::NonlinearProblem{…}, alg::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}, f::NonlinearFunction{…}, fu::Vector{…}, u::Vector{…}, p::Tuple{…}; stats::SciMLBase.NLStats, autodiff::AutoForwardDiff{…}, vjp_autodiff::AutoFiniteDiff{…}, jvp_autodiff::AutoForwardDiff{…}, linsolve::Nothing)
    @ NonlinearSolveBase C:\Users\konin_bt\.julia\packages\NonlinearSolveBase\xDyV6\src\jacobian.jl:76
 [18] construct_jacobian_cache
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolveBase\xDyV6\src\jacobian.jl:34 [inlined]
 [19] __init(::NonlinearProblem{…}, ::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveFirstOrder C:\Users\konin_bt\.julia\packages\NonlinearSolveFirstOrder\o91ZE\src\solve.jl:161        
 [20] __init
    @ C:\Users\konin_bt\.julia\packages\NonlinearSolveFirstOrder\o91ZE\src\solve.jl:121 [inlined]
 [21] #init_call#31
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:545 [inlined]
 [22] init_call
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:518 [inlined]
 [23] #init_up#34
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:571 [inlined]
 [24] init_up
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:566 [inlined]
 [25] #init#32
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:559 [inlined]
 [26] init(prob::NonlinearProblem{…}, args::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…})
    @ DiffEqBase C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:549
 [27] build_nlsolver(alg::QNDF{…}, nlalg::NonlinearSolveAlg{…}, u::Vector{…}, uprev::Vector{…}, p::SciMLBase.NullParameters, t::Float64, dt::Float64, f::ODEFunction{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, γ::Rational{…}, c::Int64, α::Int64, ::Val{…})
    @ OrdinaryDiffEqNonlinearSolve C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\utils.jl:233
 [28] build_nlsolver
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\utils.jl:152 [inlined]
 [29] build_nlsolver
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqNonlinearSolve\e4hoO\src\utils.jl:142 [inlined]
 [30] alg_cache(alg::QNDF{…}, u::Vector{…}, rate_prototype::Vector{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Vector{…}, uprev2::Vector{…}, f::ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, ::Val{…})
    @ OrdinaryDiffEqBDF C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqBDF\Ga8aD\src\bdf_caches.jl:430
 [31] __init(prob::ODEProblem{…}, alg::QNDF{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{…}, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqCore\bMOsj\src\solve.jl:409
 [32] __init (repeats 5 times)
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqCore\bMOsj\src\solve.jl:11 [inlined]
 [33] __solve(::ODEProblem{…}, ::QNDF{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqCore\bMOsj\src\solve.jl:6
 [34] __solve
    @ C:\Users\konin_bt\.julia\packages\OrdinaryDiffEqCore\bMOsj\src\solve.jl:1 [inlined]
 [35] #solve_call#35
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:635 [inlined]
 [36] solve_call(_prob::ODEProblem{…}, args::QNDF{…})
    @ DiffEqBase C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:592
 [37] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::QNDF{…}; kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:1142
 [38] solve_up
    @ C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:1120 [inlined]
 [39] solve(prob::ODEProblem{…}, args::QNDF{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:1057
 [40] solve(prob::ODEProblem{…}, args::QNDF{…})
    @ DiffEqBase C:\Users\konin_bt\.julia\packages\DiffEqBase\zYZst\src\solve.jl:1047
 [41] top-level scope
    @ c:\Users\konin_bt\OneDrive - Stichting Deltares\Documents\Ribasim_development\runners\runner.jl:38
Some type information was truncated. Use `show(err)` to see complete types.

I presume this problem was introduced by #2567.

Ping @gdalle @jClugstor

[gdalle]

I'm not familiar enough with NonlinearSolve or OrdinaryDiffEq to decipher where there should be Dual numbers and where there shouldn't be. This is probably a job for @avik-pal

[SouthEndMusic]

In this case though there shouldn't be Dual numbers anywhere (at least not in the OrdinaryDiffEq internals). I think the default AD method for implicit solvers is ForwardDiff, but there should be no AD Jacobian preparation at all when a Jacobian function is provided.

[SouthEndMusic]

It looks like this should be checked around here:

https://github.com/jClugstor/OrdinaryDiffEq.jl/blob/57ea0bc8f6b6f0ab17534aa6b1063c16f0169812/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl#L53

[jClugstor]

Am I missing something or does this just set the bottom right hand corner of the Jacobian to be 1?

function jac(J, u, p, t)
    J .= 0
    for i = 1:N
        J[N, N] = 1
    end
end

did you mean

    J .= 0
    for i = 1:N
        J[i, i] = 1
    end
end

?

Regardless, it is strange that there are Dual numbers showing up here, let me take a look.

[SouthEndMusic]

Ah you're right, but the specifics of the ODE problem aren't relevant.

[jClugstor]

Yeah, I don't think we have it set up to use the user provided jac function when you specifically request to use a NonlinearSolveAlg.

@oscardssmith do you think this is something with a quick fix , like just putting the user jac function into the NonlinearFunction here?

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl

Lines 316 to 318 in 2239fc8

	prob = NonlinearProblem(NonlinearFunction(nlf), copy(ztmp), nlp_params)
	cache = init(prob, nlalg.alg)
	nlcache = NonlinearSolveCache(

I'm not sure that the Jacobian for nlf would be exactly the same as the users jac function though.

[ChrisRackauckas]

NonlinearSolveAlg is just incomplete. Don't worry about that. That's my next big project, intertwined with the ImplicitDiscreteSolve stuff. We currently don't document this interface because it's almost guaranteed to be slower and is missing stuff right now, but I have big plans for how this + some special sauce in MTK will make the stiff ODE solves about 10x faster. That's the JuliaCon talk I have lined up an the fun paper I am working towards, and if anyone else wants to join in this work the more the marrier and I'd be happy to chat in more detail but it's not just the simple bugfix route.