Autodiff fails when initial conditions contain dual and non-dual

[hersle]

The simple example

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using ForwardDiff 

@variables x(t) y(t)
@named sys = ODESystem([D(x) ~ 0, D(y) ~ 0], t; defaults = [y => 1.0]) # change y => 1.0 to y => x, and the code runs
sys = structural_simplify(sys)

function y_at_1(x0)
    prob = ODEProblem(sys, [sys.x => x0], (0.0, 1.0), [])
    sol = solve(prob)
    return sol(1.0, idxs=y)
end

ForwardDiff.derivative(y_at_1, 1.0) # expect 0.0

gives

ERROR: LoadError: MethodError: no method matching Union{Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y_at_1), Float64}, Float64, 1}}(::Int64)

Closest candidates are:
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::DynamicQuantities.AbstractQuantity) where T<:Number
   @ DynamicQuantities C:\Users\herma\.julia\packages\DynamicQuantities\M5pvK\src\utils.jl:117
  (::Type{T})(::DynamicQuantities.AbstractGenericQuantity) where T<:Number
   @ DynamicQuantities C:\Users\herma\.julia\packages\DynamicQuantities\M5pvK\src\utils.jl:117
  ...

Stacktrace:
  [1] convert(::Type{Union{Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y_at_1), Float64}, Float64, 1}}}, x::Int64)   
    @ Base .\number.jl:7
  [2] one(::Type{Union{Float64, ForwardDiff.Dual{ForwardDiff.Tag{typeof(y_at_1), Float64}, Float64, 1}}})
    @ Base .\number.jl:347
  [3] recursive_unitless_eltype(a::Type{Union{Float64, ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}}})
    @ RecursiveArrayTools C:\Users\herma\.julia\packages\RecursiveArrayTools\xGKIm\src\utils.jl:286
  [4] recursive_unitless_eltype(a::Vector{Union{Float64, ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}}})
    @ RecursiveArrayTools C:\Users\herma\.julia\packages\RecursiveArrayTools\xGKIm\src\utils.jl:276
  [5] promote_f(f::ODEFunction{…}, ::Val{…}, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64)
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:1243
  [6] get_concrete_problem(prob::ODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:1173
  [7] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::ModelingToolkit.MTKParameters{…}; kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:1070
  [8] solve_up(::ODEProblem{…}, ::Nothing, ::Vector{…}, ::ModelingToolkit.MTKParameters{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:1066
  [9] solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:1003
 [10] solve(::ODEProblem{…})
    @ DiffEqBase C:\Users\herma\.julia\packages\DiffEqBase\V6SCE\src\solve.jl:993
 [11] y_at_1(x0::ForwardDiff.Dual{ForwardDiff.Tag{typeof(y_at_1), Float64}, Float64, 1})
    @ Main C:\Users\herma\.julia\dev\Symboltz\bug.jl:32      
 [12] derivative(f::typeof(y_at_1), x::Float64)
    @ ForwardDiff C:\Users\herma\.julia\packages\ForwardDiff\PcZ48\src\derivative.jl:14
 [13] top-level scope
    @ C:\Users\herma\.julia\dev\Symboltz\bug.jl:36
in expression starting at C:\Users\herma\.julia\dev\Symboltz\bug.jl:36

To me, the Union{Float64, ForwardDiff.Dual{...} looks like a symptom that the type promotion gets confused because x is a dual number and y is just a normal number.

[ChrisRackauckas]

To me, the Union{Float64, ForwardDiff.Dual{...} looks like a symptom that the type promotion gets confused because x is a dual number and y is just a normal number.

Yes, something happened with promotion here. It needs to promote the whole u0 vector to dual, and something went wrong with that.

[AayushSabharwal]

The real fix for this is to pass use_union = false to ODEProblem. See #2973

[hersle]

Thanks! I will try that when it is released.

Why is it a switch? Is it for performance reasons?

[ChrisRackauckas]

The default should be false? We should never default to allowing u0 to be a union? The solvers won't like a union.

[hersle]

Isn't the union handling only for the type promotion during pre-processing, so the solvers will only see pure numbers or duals? Or do I misunderstand?

[AayushSabharwal]

The default should be false? We should never default to allowing u0 to be a union? The solvers won't like a union.

Yeah, default false makes sense but the current default is true and changing it is breaking.

[ChrisRackauckas]

what's a case where it's breaking? Solvers would fail in most cases that you had unions anyways.

[AayushSabharwal]

I should've been clearer: it's breaking according to semver. I haven't tested to see if any tests fail. Will run and check

[ChrisRackauckas]

No DAE solver supports a u0 union as a state. It's not breaking because that case is already broken.