Error on arrays of static arrays

[antoine-levitt]

Hi,

I'm trying to solve an ODE that operates on Arrays of StaticArrays. I think it implements all the operations DifferentialEquations needs, so should work? MWE:

using OrdinaryDiffEq, DiffEqBase
using StaticArrays

TwoTwo = SMatrix{2,2,Float64}
x0 = zeros(TwoTwo,3,3)
function f(t,x)
    copy(x)
end

tspan = (0.0,400.0)
prob = ODEProblem(f,x0,tspan)
sol = solve(prob,Tsit5(),reltol=1e-4,abstol=1e-4)

gives

ERROR: LoadError: MethodError: no method matching ode_determine_initdt(::Array{StaticArrays.SArray{Tuple{2,2},Complex{Float64},2,4},2}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.ODEProblem{Array{SMatrix{2,2,Complex{Float64}},2},Float64,false,#f,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::Int64, ::OrdinaryDiffEq.Tsit5)
Closest candidates are:
  ode_determine_initdt(::Any, ::tType, ::Any, ::Any, ::Any, ::Any, ::Any, !Matched::DiffEqBase.AbstractODEProblem{uType,tType,true}, ::Any, ::Any) where {tType, uType} at /home/antoine/.julia/v0.6/OrdinaryDiffEq/src/initdt.jl:2
  ode_determine_initdt(::uType, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, !Matched::DiffEqBase.AbstractODEProblem{uType,tType,false}, ::Any, ::Any) where {uType, tType} at /home/antoine/.julia/v0.6/OrdinaryDiffEq/src/initdt.jl:101
Stacktrace:
 [1] auto_dt_reset! at /home/antoine/.julia/v0.6/OrdinaryDiffEq/src/integrators/integrator_interface.jl:211 [inlined]

[ChrisRackauckas]

Yes, Arrays of Arrays is something that currently gives us trouble. But most of the change is almost done. https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/231 has most of it. The last thing that we need to find out is a replacement for:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L98

We need something that gives the element type of u without units, but it needs to be compatible with the element type of u being an array. If you have a solution to that then we'll finish this once and for all 👍

[antoine-levitt]

Assume the input is an abstract array of abstract arrays of... of numbers, and recurse on that using dispatch?

[ChrisRackauckas]

Assume the input is an abstract array of abstract arrays of... of numbers, and recurse on that using dispatch?

No, that's the current strategy:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L92

defined in:

https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/blob/master/src/utils.jl#L95-L96

The issue is that for an abstract array of array of ..., this returns the bottom type, so the cache array is Array{Float64,1}, which if you look at the current error message on JuliaDiffEq/OrdinaryDiffEq.jl#231 it's that it cannot convert the static array to a Float64 because it incorrectly got the eltype of u as its number type. It needs to propagate that up back to the array type somehow, maybe via similar_type.

[ChrisRackauckas]

Related: https://stackoverflow.com/questions/47680441/recursive-unitless-element-type

[ChrisRackauckas]

Fixed by the newest release on a multitude of packages which will go through by the end of tomorrow.

[antoine-levitt]

Thanks! That resulted in a 4x speedup on a ODE in dimension 2x2x30x30 by StaticArraying the 2x2 part - not bad for such a simple change!

[ChrisRackauckas]

Sweet! I'll mark that down as a win!