Creating BVProblem from ODESystem

Project.toml

@@ -73,14 +73,15 @@ MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
 MTKHomotopyContinuationExt = "HomotopyContinuation"
-MTKLabelledArraysExt = "LabelledArrays"
 MTKInfiniteOptExt = "InfiniteOpt"
+MTKLabelledArraysExt = "LabelledArrays"
 
 [compat]
 AbstractTrees = "0.3, 0.4"
 ArrayInterface = "6, 7"
 BifurcationKit = "0.4"
 BlockArrays = "1.1"
+BoundaryValueDiffEq = "5.12.0"
 ChainRulesCore = "1"
 Combinatorics = "1"
 CommonSolve = "0.2.4"
@@ -148,6 +149,7 @@ julia = "1.9"
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
@@ -178,4 +180,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"]
+test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"]

src/systems/diffeqs/abstractodesystem.jl

@@ -469,6 +469,90 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem,
         initializeprobpmap = initializeprobpmap)
 end
 
+"""
+```julia
+SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan,
+                         parammap = DiffEqBase.NullParameters();
+                         version = nothing, tgrad = false,
+                         jac = true, sparse = true,
+                         simplify = false,
+                         kwargs...) where {iip}
+```
+
+Create a `BVProblem` from the [`ODESystem`](@ref). The arguments `dvs` and
+`ps` are used to set the order of the dependent variable and parameter vectors,
+respectively. `u0map` should be used to specify the initial condition, or be a function returning an initial condition.
+"""
+function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...)
+    BVProblem{true}(sys, args...; kwargs...)
+end
+
+function SciMLBase.BVProblem(sys::AbstractODESystem,
+        u0map::StaticArray,
+        args...;
+        kwargs...)
+    BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...)
+end
+
+function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...)
+    BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...)
+end
+
+function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...)
+    BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...)
+end
+
+function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [],
+        tspan = get_tspan(sys),
+        parammap = DiffEqBase.NullParameters();
+        version = nothing, tgrad = false,
+        callback = nothing,
+        check_length = true,
+        warn_initialize_determined = true,
+        eval_expression = false,
+        eval_module = @__MODULE__,
+        kwargs...) where {iip, specialize}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`")
+    end
+
+    f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
+        t = tspan !== nothing ? tspan[1] : tspan,
+        check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...)
+
+    cbs = process_events(sys; callback, eval_expression, eval_module, kwargs...)
+    kwargs = filter_kwargs(kwargs)
+
+    kwargs1 = (;)
+    if cbs !== nothing
+        kwargs1 = merge(kwargs1, (callback = cbs,))
+    end
+
+    # Construct initial conditions.
+    _u0 = u0 isa Function ? u0(tspan[1]) : u0
+
+    # Define the boundary conditions.
+    bc = if iip
+        (residual, u, p, t) -> (residual .= u[1] .- _u0)
+    else
+        (u, p, t) -> (u[1] - _u0)
+    end
+
+    return BVProblem{iip}(f, bc, _u0, tspan, p; kwargs1..., kwargs...)
+end
+
+get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.")
+
+@inline function SymbolicUtils.Code.create_array(::Type{<:Base.ReinterpretArray}, ::Nothing,
+        ::Val{1}, ::Val{dims}, elems...) where {dims}
+    [elems...]
+end
+
+@inline function SymbolicUtils.Code.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

@@ -0,0 +1,70 @@
+using BoundaryValueDiffEq, OrdinaryDiffEq
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+solvers = [MIRK4, RadauIIa5, LobattoIIIa3]
+
+@parameters α=7.5 β=4.0 γ=8.0 δ=5.0
+@variables x(t)=1.0 y(t)=2.0
+
+eqs = [D(x) ~ α * x - β * x * y,
+    D(y) ~ -γ * y + δ * x * y]
+
+u0map = [:x => 1.0, :y => 2.0]
+parammap = [:α => 7.5, :β => 4, :γ => 8.0, :δ => 5.0]
+tspan = (0.0, 10.0)
+
+@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; eval_expression = true)
+
+for solver in solvers
+    sol = solve(bvp, solver(), dt = 0.01)
+    @test isapprox(sol.u[end], osol.u[end]; atol = 0.01)
+    @test sol.u[1] == [1.0, 2.0]
+end
+
+# Test out of place
+bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}(
+    lotkavolterra, u0map, tspan, parammap; eval_expression = true)
+
+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.0, 2.0]
+end
+
+### Testing on pendulum
+
+@parameters g=9.81 L=1.0
+@variables θ(t) = π / 2
+
+eqs = [D(D(θ)) ~ -(g / L) * sin(θ)]
+
+@mtkbuild pend = ODESystem(eqs, t)
+
+u0map = [θ => π / 2, D(θ) => π / 2]
+parammap = [:L => 1.0, :g => 9.81]
+tspan = (0.0, 6.0)
+
+op = ODEProblem(pend, u0map, tspan, parammap)
+osol = solve(op, Vern9())
+
+bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap)
+for solver in solvers
+    sol = solve(bvp, 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)
+
+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/runtests.jl

@@ -81,6 +81,7 @@ end
             @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl")
             @safetestset "PDE Construction Test" include("pde.jl")
             @safetestset "JumpSystem Test" include("jumpsystem.jl")
+            @safetestset "BVProblem Test" include("bvproblem.jl")
             @safetestset "print_tree" include("print_tree.jl")
             @safetestset "Constraints Test" include("constraints.jl")
         end