I'm experimenting with systems of PDEs. After a semidiscretisation, I would like to store the solution as a vector of some (immutable) vectors, e.g. SVector{2,Float64}
from StaticArrays.jl
. A minimal working example is
using StaticArrays
using DiffEqBase, OrdinaryDiffEq
u0 = zeros(SVector{2,Float64}, 2) + 1
f = (t,u,du) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
This results in
No precise constructor found. Length of input was 1 while length of StaticArrays.SVector{2,Float64} is 2.
in #init#84(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Bool, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##1#2,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/solve.jl:123
in (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##1#2,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
in #solve#83(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##1#2,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/solve.jl:6
in (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##1#2,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
Here, the problem is uEltype(1)
. Thus, the problem is essentially the assumption that uEltype
is a scalar.
A first workaround may be to define a suitable constructor as follows
using StaticArrays
using DiffEqBase, OrdinaryDiffEq
(::Type{SVector{2,T}}){T}(val::Real) = SVector{2,T}(val, val)
u0 = zeros(SVector{2,Float64}, 2) + 1
f = (t,u,du) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
However, this results in
indexed assignment not defined for StaticArrays.SVector{2,Float64}
in copy!(::StaticArrays.SVector{2,Float64}, ::StaticArrays.SVector{2,Float64}) at ./multidimensional.jl:616
in recursivecopy! at /home/hendrik/.julia/v0.5/RecursiveArrayTools/src/utils.jl:2 [inlined]
in recursivecopy!(::Array{StaticArrays.SVector{2,Float64},1}, ::Array{StaticArrays.SVector{2,Float64},1}) at /home/hendrik/.julia/v0.5/RecursiveArrayTools/src/utils.jl:7
in apply_step! at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/integrators/integrator_utils.jl:247 [inlined]
in loopheader! at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/integrators/integrator_utils.jl:10 [inlined]
in solve!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Euler,Array{StaticArrays.SVector{2,Float64},1},Float64,Float64,Float64,Array{Array{StaticArrays.SVector{2,Float64},1},1},DiffEqBase.ODESolution{StaticArrays.SVector{2,Float64},2,Array{Array{StaticArrays.SVector{2,Float64},1},1},Void,Void,Array{Float64,1},Array{Array{Array{StaticArrays.SVector{2,Float64},1},1},1},DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##3#4,Void,UniformScaling{Int64}},OrdinaryDiffEq.Euler,OrdinaryDiffEq.InterpolationData{##3#4,Array{Array{StaticArrays.SVector{2,Float64},1},1},Array{Float64,1},Array{Array{Array{StaticArrays.SVector{2,Float64},1},1},1},OrdinaryDiffEq.EulerCache{Array{StaticArrays.SVector{2,Float64},1},Array{StaticArrays.SVector{2,Float64},1}}}},Array{StaticArrays.SVector{2,Float64},1},##3#4,Void,OrdinaryDiffEq.EulerCache{Array{StaticArrays.SVector{2,Float64},1},Array{StaticArrays.SVector{2,Float64},1}},OrdinaryDiffEq.DEOptions{StaticArrays.SVector{2,Float64},Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void}}) at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/solve.jl:306
in #solve#83(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##3#4,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/hendrik/.julia/v0.5/OrdinaryDiffEq/src/solve.jl:7
in (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Array{StaticArrays.SVector{2,Float64},1},Float64,true,##3#4,Void,UniformScaling{Int64}}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
Hence, the assumption that uEltype
is a scalar appears again.
However, if I use an algorithm from ODE.jl
,
using StaticArrays
using DiffEqBase, OrdinaryDiffEq, ODE
u0 = zeros(SVector{2,Float64}, 2) + 1
f = (t,u,du) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, ode4(), dt=1.e-2)
sol[end]
everything works fine, the result is
2-element Array{StaticArrays.SVector{2,Float64},1}:
[2.70833,2.70833]
[2.70833,2.70833]
However, I would really like to use the algorithms from OrdinaryDiffEq.jl
- of course. I can invest some work and time to get things done. Since the implicit assumption uEltype <: Number
seems to be widespread, some guidance from more experienced developers of the DiffEq suite would be very nice.
What do you think of something like uScalarType
instead of uEltype
for tolerances etc.? Additionally, the copying in RecursiveArrayTools
would have to be adapted.