Adding saveat to DifferentialEquations.jl solve causes conversion error when using custom variable type

General Usage
,
1

I’ve got a simple struct that looks something like this:

mutable struct variable{T,N} <: AbstractArray{T,N}
    values   :: Array{T,N}
    spaceIDs :: Vector{Int}
end

I’ve defined Base.:+, as well as Base.:* and Base:./(with scalar), as well as the basic AbstractArray interface functions. I define an ODEProblem with an in-place function (that does nothing) and pass it to solve. It works great (even when the passed problem actually does something to the variables) if I don’t include the saveat option. However soln = solve(prob,saveat=0.1); throws a conversion error with the start of the stacktrace:

Stacktrace:
  [1] push!(a::Vector{AMPS.vars.variable{Float64, 2}}, item::Matrix{Float64})
    @ Base ./array.jl:928
  [2] copyat_or_push!(a::Vector{AMPS.vars.variable{Float64, 2}}, i::Int64, x::Matrix{Float64}, nc::Type{Val{false}})
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/e9IdM/src/utils.jl:79
  [3] _savevalues!(integrator::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5, true, AMPS.vars.variable{Float64, 2}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{AMPS.vars.variable{Float64, 2}}, SciMLBase.ODESolution{Float64, 3, Vector{AMPS.vars.variable{Float64, 2}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{AMPS.vars.variable{Float64, 2}}}, SciMLBase.ODEProblem{AMPS.vars.variable{Float64, 2}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, Main.finitediff.var"#test2D#2", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.Tsit5, OrdinaryDiffEq.InterpolationData{SciMLBase.ODEFunction{true, Main.finitediff.var"#test2D#2", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{AMPS.vars.variable{Float64, 2}}, Vector{Float64}, Vector{Vector{AMPS.vars.variable{Float64, 2}}}, OrdinaryDiffEq.Tsit5Cache{AMPS.vars.variable{Float64, 2}, AMPS.vars.variable{Float64, 2}, AMPS.vars.variable{Float64, 2}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, SciMLBase.ODEFunction{true, Main.finitediff.var"#test2D#2", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{AMPS.vars.variable{Float64, 2}, AMPS.vars.variable{Float64, 2}, AMPS.vars.variable{Float64, 2}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, DiffEqBase.CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Float64, Tuple{}}, AMPS.vars.variable{Float64, 2}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, force_save::Bool, reduce_size::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AgIKX/src/integrators/integrator_utils.jl:64

Looks like it tries to push! a matrix instead of my defined struct. Anyone have an idea why this might be happening? I guess that the solver at some point interprets my array type as a standard matrix, but why would this happen only when I specify saveat?

2

Adding MWE:

using Lazy, DifferentialEquations

mutable struct variable{T,N} <: AbstractArray{T,N}
    values   :: Array{T,N}
    spaceIDs :: Vector{Int}
end

@forward variable.values (Base.size, Base.getindex, Base.setindex!, Base.length)

Base.similar(v::variable) = variable(similar(v.values),v.spaceIDs)
Base.similar(v::variable,::Type{T}) where{T} = variable(similar(v.values,T),v.spaceIDs)

function Base.:+(x::variable{T,N},y::variable{T,N}) where {T,N}
    all(x.spaceIDs == y.spaceIDs) ? variable(x.values+y.values,x.spaceIDs) : error("Attempted to add variables defined on different spaces")
end

Base.:*(x::Number,y::variable) = variable(x*y.values,y.spaceIDs)
Base.:*(x::variable,y::Number) = variable(y*x.values,x.spaceIDs)
Base.:/(x::variable,y::Number) = variable(x.values/y,x.spaceIDs)

u0mat = zeros(100,100)
u0 = variable(u0mat,[1,2])
du = variable(u0mat,[1,2])

function test(du,u,p,t) end 

prob = ODEProblem(test, u0, (0.,5.));
soln = solve(prob,saveat=0.1);
3

Broadcast definitions are probably missing. https://github.com/SciML/DEDataArrays.jl probably would help.

4

It was the broadcast definitions, thanks! Added the first two from here https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting and it no longer breaks when using saveat.