AbstractArray support metaissue

[ChrisRackauckas]

I went through to find out what's required for GPUArrays to work well in the *DiffEq packages. Essentially, other than a few GPUArray bugs:

JuliaGPU/GPUArrays.jl#20
JuliaGPU/GPUArrays.jl#21
JuliaGPU/GPUArrays.jl#22

The only issues that come up is that GPUArrays are an AbstractArray type which does not have indexing. Incredibly, in the broadcast branches / masters this mostly works, except for one line. However, that would be fixed by the following tiny new feature in Julia's Base:

JuliaLang/julia#22216

Then the *DiffEq solvers can automatically compile new versions of themselves for the GPU.

(I think this is super cool haha)

[ChrisRackauckas]

DistributedArrays work in the out-of-place form. For the in-place form to work, they just need:

JuliaParallel/DistributedArrays.jl#141

[ChrisRackauckas]

SharedArrays fail in the out-of-place form since u0 .* u0 returns an Array instead of a SharedArray. But there's pretty much no reason to use them out-of-place anyways. They fail in the in-place form since similar(::SharedArray)::Array for some reason, but that's easy to spotfix if necessary.

[ChrisRackauckas]

using GPUArrays, RecursiveArrayTools

u0 = GPUArray(rand(Float32, 32, 32))
v0 = GPUArray(rand(Float32, 32, 32))
uv0 = ArrayPartition(u0,v0)

function f(t,y,dy)
    u,v = y.x
    du,dv = dy.x
    du .= v
    dv .= -u
end
tspan = (Float32(0.0),Float32(1.0))
using OrdinaryDiffEq
prob = ODEProblem(f,uv0,tspan)
sol = solve(prob,BS3())

Almost works in GPUArrays except for JuliaGPU/GPUArrays.jl#41 . But that means that

sol = solve(prob,BS3(),dt=0.01,adaptive=false)

works! Additionally

using OrdinaryDiffEq
u0 = GPUArray(rand(Float32, 32, 32))
v0 = GPUArray(rand(Float32, 32, 32))
f1 = function (t,u,v,du)
  du .= v
end
f2 = function (t,u,v,dv)
  dv .= -u
end
function (::typeof(f2))(::Type{Val{:analytic}}, x, y0)
  u0, v0 = y0
  ArrayPartition(u0*cos(x) + v0*sin(x), -u0*sin(x) + v0*cos(x))
end

prob = DynamicalODEProblem(f1,f2,u0,v0,(Float32(0.0),Float32(5.0)))
sol = solve(prob,McAte4(),dt=1/10)

works.

[ChrisRackauckas]

Tracking not broadcasting: SciML/OrdinaryDiffEq.jl#106

[ChrisRackauckas]

It's all working. No need for the metaissue and we can tackle individual issues.