DifferentialEquations.jl and ReverseDiff.jl: MethodError: Cannot `convert` an object of type Vector{ReverseDiff.TrackedReal{

Specific Domains Numerics
1

I am trying to get gradients in a heat equation solved with DifferentialEquations.jl:

using DiffEqOperators, OrdinaryDiffEq, Plots, Zygote, DiffEqSensitivity

nknots = 10
h = 1.0 / (nknots + 1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2

const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots - 2)

τ = 0.1
ξ = 0.5
t0 = 0.0
tlimit = 1.0 / ξ
t1 = 2.0 * tlimit
vlimit = 1.0
p = [1.0, 2.0]

function heat!(res, du, u, p, t)
    res[1] = u[1] - t
    res[2:end - 1] = p[1] * p[2] * du[2:end - 1] - 1.0 * Δ *  u
    res[end] = u[end] - u[end - 1]
end

u0 = zeros(nknots)
du0 = zeros(nknots)
prob = DAEProblem{true}(heat!, du0, u0, (t0, t1), p) # , differential_vars=diff_var)
sol = solve(prob, DABDF2(), abstol=1e-3, reltol=1e-3)

function sum_of_solution(p)
    _prob = remake(prob, p=p)
    sum(solve(_prob, DABDF2(), rtol=1e-6, atol=1e-6, saveat=0.01, sensealg=ReverseDiffAdjoint()))
end
println(sum_of_solution(p))
dp1 = Zygote.gradient(sum_of_solution,  p)

but I am getting the error

ERROR: LoadError: MethodError: Cannot `convert` an object of type Vector{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}} to an object of type ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}
Closest candidates are:
  convert(::Type{T}, ::LinearAlgebra.Factorization) where T<:AbstractArray at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/factorization.jl:58
  convert(::Type{T}, ::T) where T<:ReverseDiff.TrackedArray at /Users/salazardetro1/.julia/packages/ReverseDiff/E4Tzn/src/tracked.jl:270
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:14

Is this because ReverseDiff does not support Vectors? What are my options to get the gradients? Thanks.

2

With two parameters forward-mode AD is likely faster anyways, so sensealg=ForwardDiffSensitivty() is likely a good idea. We need to do more on fully implicit DAEs. The tools are much more optimized right now for Index-1 DAEs in mass matrix form, so I would recommend you use that formulation if you can (and that will generally be more efficient even after the fully implicit form is more optimized BTW).

That said, I’d still like to look into this so if you can open an issue on DiffEqSensitivity with this example I might look over it next week. ReverseDiffAdjoint should work here, though since it’s mutating it won’t be a fast way to take derivatives since it will need to do an Array of Structs form.