DimensionMismatch error in Jacobian for stiff SecondOrderODE systems

[luismunozheinen]

I am trying to solve a stiff system of 200 differential equations of type du=A*u+f(u). sol = solve(prob) runs fine as well as specifying any non-stiff solvers. However specifying alg_hints=[:stiff] or any stiff solver (shown below for Rodas5()) results in

ERROR: DimensionMismatch("parent has 200 elements, which is incompatible with size (200, 200)") Stacktrace: [1] _reshape(::RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}, ::Tuple{Int64,Int64}) at ./reshapedarray.jl:141 [2] reshape_jacobian(::RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}, ::RecursiveArrayTools.ArrayPartition{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2}}}, ::RecursiveArrayTools.ArrayPartition{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2}}}) at ... julia/v0.6/ForwardDiff/src/jacobian.jl:131 [3] chunk_mode_jacobian(::DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void}, ::RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10,RecursiveArrayTools.ArrayPartition{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2}}}}) at ...julia/v0.6/ForwardDiff/src/jacobian.jl:214 [4] jacobian(::Function, ::RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}, ::ForwardDiff.JacobianConfig{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10,RecursiveArrayTools.ArrayPartition{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},Tuple{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},Float64},Float64,10},2}}}}, ::Val{true}) at ...julia/v0.6/ForwardDiff/src/jacobian.jl:19 [5] calc_W!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Float64,Void,Float64,Float64,Float64,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},DiffEqBase.ODESolution{Float64,2,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Void,Void,Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},DiffEqBase.ODEProblem{Tuple{Array{Float64,2},Array{Float64,2}},Float64,false,Void,DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Void,Void,UniformScaling{Int64},DiffEqBase.SecondOrderODEProblem{false}},OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},OrdinaryDiffEq.InterpolationData{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}}}},DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Bool,OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{#condition,#affect!,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Array{Float64,1},Float64,Array{Float64,1}},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}}, ::OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}}, ::Float64, ::Bool, ::Bool) at ...julia/v0.6/OrdinaryDiffEq/src/derivative_utils.jl:99 [6] perform_step!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Float64,Void,Float64,Float64,Float64,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},DiffEqBase.ODESolution{Float64,2,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Void,Void,Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},DiffEqBase.ODEProblem{Tuple{Array{Float64,2},Array{Float64,2}},Float64,false,Void,DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Void,Void,UniformScaling{Int64},DiffEqBase.SecondOrderODEProblem{false}},OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},OrdinaryDiffEq.InterpolationData{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}}}},DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Bool,OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{#condition,#affect!,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Array{Float64,1},Float64,Array{Float64,1}},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}}, ::OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}}, ::Bool) at ...julia/v0.6/OrdinaryDiffEq/src/perform_step/rosenbrock_perform_step.jl:1030 [7] solve!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Float64,Void,Float64,Float64,Float64,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},DiffEqBase.ODESolution{Float64,2,Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Void,Void,Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},DiffEqBase.ODEProblem{Tuple{Array{Float64,2},Array{Float64,2}},Float64,false,Void,DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Void,Void,UniformScaling{Int64},DiffEqBase.SecondOrderODEProblem{false}},OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType},OrdinaryDiffEq.InterpolationData{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},Array{Float64,1},Array{Array{RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},1},1},OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}}}},DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Bool,OrdinaryDiffEq.Rosenbrock5ConstantCache{DiffEqDiffTools.TimeDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}},Void},DiffEqDiffTools.UDerivativeWrapper{DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Float64,Void},OrdinaryDiffEq.Rodas5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{#condition,#affect!,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Float64,Array{Float64,1},Float64,Array{Float64,1}},RecursiveArrayTools.ArrayPartition{Float64,Tuple{Array{Float64,2},Array{Float64,2}}}}) at ...julia/v0.6/OrdinaryDiffEq/src/solve.jl:337 [8] #solve#1239(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Tuple{Array{Float64,2},Array{Float64,2}},Float64,false,Void,DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Void,Void,UniformScaling{Int64},DiffEqBase.SecondOrderODEProblem{false}}, ::OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ...julia/v0.6/OrdinaryDiffEq/src/solve.jl:7 [9] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Tuple{Array{Float64,2},Array{Float64,2}},Float64,false,Void,DiffEqBase.DynamicalODEFunction{false,#f1,DiffEqBase.##193#195},Void,Void,UniformScaling{Int64},DiffEqBase.SecondOrderODEProblem{false}}, ::OrdinaryDiffEq.Rodas5{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)

I cannot reproduce this error with the ODE examples from the DiffEq documentations. I am using Julia 0.6.2 and DifferentialEquations 4.5.0. Thank you for any suggestions.

[ChrisRackauckas]

It looks like you're not using an ODEProblem but a DynamicalODEProblem of sorts? Like a SecondOrderODEProblem?

[luismunozheinen]

Indeed. I did not immediately find out from the error message that it was related to the use SecondOrderODEProblem. Considering the documentation Dynamical, Hamiltonian, and 2nd Order ODE Solvers I presume that the suited way to solve the stiff dynamical ODE is via transformation back to a system of first order ODEs in order to specify stiff solvers

[ChrisRackauckas]

It looks like this is a bug in the stiff solver handling of ArrayPartitions. If you code the 1st order system yourself it should work, but the automated handling likely needs v0.7 to fix this correctly (because of how it changes some AbstractArray type handling).

[ChrisRackauckas]

Fixed by https://github.com/SciML/DiffEqBase.jl/pull/536