diff --git a/README.md b/README.md index 2a7f13c0..523d1883 100644 --- a/README.md +++ b/README.md @@ -18,7 +18,7 @@ functionality should check out [DifferentialEquations.jl](https://github.com/Jul ## API -BoundaryValueDiffEq.jl is part of the JuliaDiffEq common interface, but can be used independently of DifferentialEquations.jl. The only requirement is that the user passes a BoundaryValueDiffEq.jl algorithm to solve. For example, we can solve the [BVP tutorial from the documentation](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/bvp_example/) using the `MIRK4()` algorithm: +BoundaryValueDiffEq.jl is part of the SciML common interface, but can be used independently of DifferentialEquations.jl. The only requirement is that the user passes a BoundaryValueDiffEq.jl algorithm to solve. For example, we can solve the [BVP tutorial from the documentation](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/bvp_example/) using the `MIRK4()` algorithm: ```julia using BoundaryValueDiffEq diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl index fde4f655..3710ae6a 100644 --- a/src/BoundaryValueDiffEq.jl +++ b/src/BoundaryValueDiffEq.jl @@ -76,10 +76,12 @@ end u0 = [5.0, -3.5] bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) - probs = [BVProblem(f1!, bc1!, u0, tspan; nlls=Val(false)), - BVProblem(f1, bc1, u0, tspan; nlls=Val(false)), - TwoPointBVProblem(f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype, nlls=Val(false)), - TwoPointBVProblem(f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype, nlls=Val(false))] + probs = [BVProblem(f1!, bc1!, u0, tspan; nlls = Val(false)), + BVProblem(f1, bc1, u0, tspan; nlls = Val(false)), + TwoPointBVProblem( + f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype, nlls = Val(false)), + TwoPointBVProblem( + f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype, nlls = Val(false))] algs = [] @@ -130,8 +132,8 @@ end u0, tspan, nlls = Val(true)), BVProblem(BVPFunction(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1), u0, tspan, nlls = Val(true)), - TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, - tspan; bcresid_prototype = bcresid_prototype2, nlls = Val(true)), + TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan; + bcresid_prototype = bcresid_prototype2, nlls = Val(true)), TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan; bcresid_prototype = bcresid_prototype2, nlls = Val(true))] @@ -240,7 +242,8 @@ end if @load_preference("PrecompileShootingNLLS", true) append!(algs, [ - Shooting(Tsit5(); nlsolve = LevenbergMarquardt(; disable_geodesic = Val(true)), + Shooting( + Tsit5(); nlsolve = LevenbergMarquardt(; disable_geodesic = Val(true)), jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))), Shooting(Tsit5(); nlsolve = GaussNewton(), jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))]) diff --git a/src/types.jl b/src/types.jl index c09e2922..d0982fed 100644 --- a/src/types.jl +++ b/src/types.jl @@ -98,8 +98,7 @@ end function concrete_jacobian_algorithm( jac_alg::BVPJacobianAlgorithm, prob_type, prob::BVProblem, alg) - u0 = prob.u0 isa AbstractArray ? prob.u0 : - __initial_guess(prob.u0, prob.p, first(prob.tspan)) + u0 = __extract_u0(prob.u0, prob.p, first(prob.tspan)) diffmode = jac_alg.diffmode === nothing ? __default_sparse_ad(u0) : jac_alg.diffmode bc_diffmode = jac_alg.bc_diffmode === nothing ? (prob_type isa TwoPointBVProblem ? __default_sparse_ad : diff --git a/src/utils.jl b/src/utils.jl index 0576e9e4..bfa2d101 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -132,6 +132,12 @@ function __extract_problem_details(prob, u0::AbstractVector{<:AbstractArray}; kw _u0 = first(u0) return Val(true), eltype(_u0), length(_u0), (length(u0) - 1), _u0 end +function __extract_problem_details( + prob, u0::RecursiveArrayTools.AbstractVectorOfArray; kwargs...) + # Problem has Initial Guess + _u0 = first(u0.u) + return Val(true), eltype(_u0), length(_u0), (length(u0.u) - 1), _u0 +end function __extract_problem_details( prob, u0::AbstractArray; dt = 0.0, check_positive_dt::Bool = false) # Problem does not have Initial Guess diff --git a/test/mirk/mirk_basic_tests.jl b/test/mirk/mirk_basic_tests.jl index 1fcf4cc4..c80e380b 100644 --- a/test/mirk/mirk_basic_tests.jl +++ b/test/mirk/mirk_basic_tests.jl @@ -213,3 +213,40 @@ end @test_nowarn solve(prob, MIRK4(); dt = 0.01) end + +@testitem "Solve using Continuation" begin + using RecursiveArrayTools + + g = 9.81 + L = 1.0 + tspan = (0.0, pi / 2) + function simplependulum!(du, u, p, t) + θ = u[1] + dθ = u[2] + du[1] = dθ + du[2] = -(g / L) * sin(θ) + end + + function bc2a!(resid_a, u_a, p) # u_a is at the beginning of the time span + x0 = p + resid_a[1] = u_a[1] - x0 # the solution at the beginning of the time span should be -pi/2 + end + function bc2b!(resid_b, u_b, p) # u_b is at the ending of the time span + x0 = p + resid_b[1] = u_b[1] - pi / 2 # the solution at the end of the time span should be pi/2 + end + + bvp3 = TwoPointBVProblem( + simplependulum!, (bc2a!, bc2b!), [pi / 2, pi / 2], (pi / 4, pi / 2), + -pi / 2; bcresid_prototype = (zeros(1), zeros(1))) + sol3 = solve(bvp3, MIRK4(), dt = 0.05) + + # Needs a SciMLBase fix + bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2), + pi / 2; bcresid_prototype = (zeros(1), zeros(1))) + @test_broken solve(bvp4, MIRK4(), dt = 0.05) isa SciMLBase.ODESolution + + bvp5 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), DiffEqArray(sol3.u, sol3.t), + (0, pi / 2), pi / 2; bcresid_prototype = (zeros(1), zeros(1))) + @test SciMLBase.successful_retcode(solve(bvp5, MIRK4(), dt = 0.05).retcode) +end diff --git a/test/runtests.jl b/test/runtests.jl index 8ba7978a..68469e0f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,8 @@ using ReTestItems -ReTestItems.runtests(@__DIR__) +ReTestItems.runtests(joinpath(@__DIR__, "mirk/")) +ReTestItems.runtests(joinpath(@__DIR__, "misc/")) +ReTestItems.runtests(joinpath(@__DIR__, "shooting/")) + +# Wrappers like ODEInterface don't support parallel testing +ReTestItems.runtests(joinpath(@__DIR__, "wrappers/"); nworkers = 0) diff --git a/test/wrappers/odeinterface_tests.jl b/test/wrappers/odeinterface_tests.jl index 2e253737..e6c5b680 100644 --- a/test/wrappers/odeinterface_tests.jl +++ b/test/wrappers/odeinterface_tests.jl @@ -25,19 +25,15 @@ u0 = [0.5, 1.0] p = [0.1] tspan = (-π / 2, π / 2) -tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, - p; bcresid_prototype = (zeros(1), zeros(1))) - -# Just generate a solution for bvpsol -sol_ms = solve(tpprob, MultipleShooting(10, DP5(), NewtonRaphson()); - dt = π / 20, abstol = 1e-5, maxiters = 1000, adaptive = false) - -export ex7_f!, ex7_2pbc1!, ex7_2pbc2!, u0, p, tspan, tpprob, sol_ms +export ex7_f!, ex7_2pbc1!, ex7_2pbc2!, u0, p, tspan end @testitem "BVPM2" setup=[ODEInterfaceWrapperTestSetup] begin - using ODEInterface, RecursiveArrayTools + using ODEInterface, RecursiveArrayTools, LinearAlgebra + + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, + p; bcresid_prototype = (zeros(1), zeros(1))) sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) @test SciMLBase.successful_retcode(sol_bvpm2) @@ -51,6 +47,13 @@ end @testitem "BVPSOL" setup=[ODEInterfaceWrapperTestSetup] begin using ODEInterface, RecursiveArrayTools + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, + p; bcresid_prototype = (zeros(1), zeros(1))) + + # Just generate a solution for bvpsol + sol_ms = solve(tpprob, MultipleShooting(10, DP5(), NewtonRaphson()); + dt = π / 20, abstol = 1e-5, maxiters = 1000, adaptive = false) + initial_u0 = [sol_ms(t) .+ rand() for t in tspan[1]:(π / 20):tspan[2]] tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1)))