fixing create_array

Author: vyudu

Project.toml

@@ -7,6 +7,9 @@ version = "9.54.0"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
+BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
+BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

src/systems/diffeqs/abstractodesystem.jl

@@ -539,11 +539,19 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = []
         (u, p, t) -> (u[1] - _u0)
     end
 
-    return BVProblem{iip}(f, bc, u0, tspan, p; kwargs1..., kwargs...)
+    return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...)
 end
 
 get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.")
 
+@inline function create_array(::Type{Base.ReinterpretArray}, ::Nothing, ::Val{1}, ::Val{dims}, elems...) where dims
+    [elems...]
+end
+
+@inline function create_array(::Type{Base.ReinterpretArray}, T, ::Val{1}, ::Val{dims}, elems...) where dims
+    T[elems...]
+end
+
 """
 ```julia
 DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys),

test/bvproblem.jl

@@ -2,6 +2,8 @@ using BoundaryValueDiffEq, OrdinaryDiffEq
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
+solvers = [MIRK4, RadauIIa5, LobattoIIIa3]
+
 @parameters α = 7.5 β = 4. γ = 8. δ = 5. 
 @variables x(t) = 1. y(t) = 2. 
 
@@ -13,20 +15,26 @@ parammap = [:α => 7.5, :β => 4, :γ => 8., :δ => 5.]
 tspan = (0., 10.)
 
 @mtkbuild lotkavolterra = ODESystem(eqs, t)
+op = ODEProblem(lotkavolterra, u0map, tspan, parammap)
+osol = solve(op, Vern9())
 
-bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap)
-sol = solve(bvp, MIRK4(), dt = 0.01);
+bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true)
 
-bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap)
-sol2 = solve(bvp, MIRK4(), dt = 0.01);
+for solver in solvers
+    println("$solver")
+    sol = solve(bvp, solver(), dt = 0.01)
+    @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
+    @test sol.u[1] == [1., 2.]
+end
 
-op = ODEProblem(lotkavolterra, u0map, tspan, parammap)
-osol = solve(op, Vern9())
+# Test out of place
+bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(lotkavolterra, u0map, tspan, parammap; eval_expression = true)
 
-@test isapprox(sol.u[end],osol.u[end]; atol = 0.01)
-@test isapprox(sol2.u[end],osol.u[end]; atol = 0.01)
-@test sol.u[1] == [1., 2.]
-@test sol2.u[1] == [1., 2.]
+for solver in solvers
+    sol = solve(bvp2, solver(), dt = 0.01)
+    @test isapprox(sol.u[end],osol.u[end]; atol = 0.01)
+    @test sol.u[1] == [1., 2.]
+end
 
 ### Testing on pendulum
 
@@ -38,19 +46,24 @@ eqs = [D(D(θ)) ~ -(g / L) * sin(θ)]
 @mtkbuild pend = ODESystem(eqs, t)
 
 u0map = [θ => π/2, D(θ) => π/2]
-parammap = [:L => 2., :g => 9.81]
+parammap = [:L => 1., :g => 9.81]
 tspan = (0., 10.)
 
+op = ODEProblem(pend, u0map, tspan, parammap)
+osol = solve(op, Vern9())
+
 bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap)
-sol = solve(bvp, MIRK4(), dt = 0.01);
+for solver in solvers
+    sol = solve(bvp2, solver(), dt = 0.01)
+    @test isapprox(sol.u[end],osol.u[end]; atol = 0.01)
+    @test sol.u[1] == [π/2, π/2]
+end
 
+# Test out-of-place
 bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap)
-sol2 = solve(bvp2, MIRK4(), dt = 0.01);
-
-op = ODEProblem(pend, u0map, tspan, parammap)
-osol = solve(op, Vern9())
 
-@test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
-@test sol.u[1] == [π/2, π/2]
-@test isapprox(sol2.u[end], osol.u[end]; atol = 0.01)
-@test sol2.u[1] == [π/2, π/2]
+for solver in solvers
+    sol = solve(bvp2, solver(), dt = 0.01)
+    @test isapprox(sol.u[end],osol.u[end]; atol = 0.01)
+    @test sol.u[1] == [π/2, π/2]
+end