Allow uEltype to be an immutable vector, e.g. SVector from StaticArrays.jl
#50
Comments
That's fixed by a
Yeah, this should get worked out using |
|
I just noticed something though: you don't want to use this with stiff methods (unless you use the https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/solve.jl#L4 Just to keep in mind. |
|
Great, thanks. Actually, I'm more interested in hyperbolic PDEs where the explicit strong-stability preserving SSPRK methods will be used. |
|
Yeah, I'm getting this done. Actually, I built RecursiveArrayTools.jl at first to make this exact use case work. Finishing this up shouldn't be bad. I have MArrays working. |
|
I got it working. A few things to convert to recursive calls, and the definitions: function RecursiveArrayTools.recursivecopy!{T<:StaticArray,N}(b::AbstractArray{T,N},a::AbstractArray{T,N})
@inbounds for i in eachindex(a)
b[i] = a[i]
end
end
function RecursiveArrayTools.recursivecopy!{T<:MArray,N}(b::AbstractArray{T,N},a::AbstractArray{T,N})
@inbounds for i in eachindex(a)
recursivecopy!(b[i],a[i])
end
endThis means StaticArrays needs to be added to RecursiveArrayTools.jl for this to work. But I'll allow it since it's a Julia-only library. The real solution, which wouldn't require a dependency on each array type like this, would need a Base change: For now, we should be good though. I'll get this all tagged. |
using StaticArrays
using DiffEqBase, OrdinaryDiffEq
using RecursiveArrayTools
u0 = zeros(MVector{2,Float64}, 2) + 1
u0[1] = ones(MVector{2,Float64}) + 1
f = (t,u,du) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
u0 = zeros(SVector{2,Float64}, 2) + 1
u0[1] = ones(SVector{2,Float64}) + 1
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
sol = solve(ode, SSPRK22(), dt=1.e-2)
u0 = zero(MVector{2,Float64}) + 1
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)
u0 = zero(SVector{2,Float64}) + 1
f = (t,u) -> u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1.e-2)Looks good. We may even want to add this to the starting tutorial, since |
|
All works, except adaptive timestepping. That needs a smarter default norm. You can pass a norm in |
|
Great work, thanks a lot. My examples are working now. I'm prototyping some ideas that are planned to become something like |
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}fromStaticArrays.jl. A minimal working example isThis results in
Here, the problem is
uEltype(1). Thus, the problem is essentially the assumption thatuEltypeis a scalar.A first workaround may be to define a suitable constructor as follows
However, this results in
Hence, the assumption that
uEltypeis a scalar appears again.However, if I use an algorithm from
ODE.jl,everything works fine, the result is
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 assumptionuEltype <: Numberseems to be widespread, some guidance from more experienced developers of the DiffEq suite would be very nice.What do you think of something like
uScalarTypeinstead ofuEltypefor tolerances etc.? Additionally, the copying inRecursiveArrayToolswould have to be adapted.The text was updated successfully, but these errors were encountered: