diff --git a/Manifest.toml b/Manifest.toml index 23df17cd..974c3141 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -808,7 +808,7 @@ version = "0.6.42" [[deps.SciMLBase]] deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLOperators", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "Tables"] -git-tree-sha1 = "0e4fca3dd5de4d4a82c0ffae1e51ab6234af4df0" +git-tree-sha1 = "3fbd6a361ec965a89f1ec64320a8b8a80b7409c3" repo-rev = "ap/nlls_bvp" repo-url = "https://github.com/SciML/SciMLBase.jl.git" uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462" diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqODEInterfaceExt.jl index a601bef2..9e7a76ab 100644 --- a/ext/BoundaryValueDiffEqODEInterfaceExt.jl +++ b/ext/BoundaryValueDiffEqODEInterfaceExt.jl @@ -14,14 +14,6 @@ import ODEInterface: colnew import FastClosures: @closure import ForwardDiff -function _test_bvpm2_bvpsol_colnew_problem_criteria( - _, ::SciMLBase.StandardBVProblem, alg::Symbol) - throw(ArgumentError("$(alg) does not support standard BVProblem. Only TwoPointBVProblem is supported.")) -end -function _test_bvpm2_bvpsol_colnew_problem_criteria(prob, ::TwoPointBVProblem, alg::Symbol) - @assert isinplace(prob) "$(alg) only supports inplace TwoPointBVProblem!" -end - #------ # BVPM2 #------ @@ -54,7 +46,7 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPM2; dt = 0.0, reltol = 1e-3, if prob.u0 isa Function guess_function = @closure (x, y) -> (y .= vec(__initial_guess(prob.u0, prob.p, x))) bvpm2_init(obj, no_odes, no_left_bc, mesh, guess_function, - eltype(u0_)[], alg.max_num_subintervals, prob.u0) + eltype(u0_)[], alg.max_num_subintervals) else u0 = __flatten_initial_guess(prob.u0) bvpm2_init( @@ -98,7 +90,7 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPM2; dt = 0.0, reltol = 1e-3, bvpm2_destroy(obj) bvpm2_destroy(sol) - return SciMLBase.build_solution(prob, ivpsol, nothing) + return ivpsol end #------- @@ -171,9 +163,9 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPSOL; maxiters = 1000, elseif retcode == -5 @warn "Given initial values inconsistent with separable linear bc" elseif retcode == -6 - @warn """Iterative refinement faild to converge for `sol_method=0` - Termination since multiple shooting condition or - condition of Jacobian is too bad for `sol_method=1`""" + @warn "Iterative refinement faild to converge for `sol_method=0` \ + Termination since multiple shooting condition or \ + condition of Jacobian is too bad for `sol_method=1`" elseif retcode == -8 @warn "Condensing algorithm for linear block system fails, try `sol_method=1`" elseif retcode == -9 @@ -187,68 +179,92 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPSOL; maxiters = 1000, ivpsol = SciMLBase.build_solution(prob, alg, sol_t, map(x -> reshape(convert(Vector{eltype(u0_)}, x), u0_size), eachcol(sol_x)); - retcode = retcode ≥ 0 ? ReturnCode.Success : ReturnCode.Failure, stats, - original = (sol_t, sol_x, retcode, stats)) + retcode = retcode ≥ 0 ? ReturnCode.Success : ReturnCode.Failure, + stats, original = (sol_t, sol_x, retcode, stats)) - return SciMLBase.build_solution(prob, ivpsol, nothing) + return ivpsol end #------- # COLNEW #------- -#= TODO: FIX this -function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol=1e-4, dt = 0.0, verbose = true, kwargs...) - _test_bvpm2_bvpsol_colnew_problem_criteria(prob, prob.problem_type, :COLNEW) - has_initial_guess = prob.u0 isa AbstractVector{<:AbstractArray} - dt ≤ 0 && throw(ArgumentError("dt must be positive")) - no_odes, n, u0 = if has_initial_guess - length(first(prob.u0)), (length(prob.u0) - 1), reduce(hcat, prob.u0) - else - length(prob.u0), Int(cld((prob.tspan[2] - prob.tspan[1]), dt)), prob.u0 +function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, + reltol = 1e-3, dt = 0.0, verbose = true, kwargs...) + # FIXME: COLNEW does support MP-BVPs but in a very clunky way + if !(prob.problem_type isa TwoPointBVProblem) + throw(ArgumentError("`COLNEW` only supports `TwoPointBVProblem!`")) end + + dt ≤ 0 && throw(ArgumentError("`dt` must be positive")) + + t₀, t₁ = prob.tspan + u0_ = __extract_u0(prob.u0, prob.p, t₀) + u0_size = size(u0_) + n = __initial_guess_length(prob.u0) + + u0 = __flatten_initial_guess(prob.u0) + mesh = __extract_mesh(prob.u0, t₀, t₁, ifelse(n == -1, dt, n - 1)) + if u0 === nothing + # initial_guess function was provided + u0 = mapreduce(@closure(t->vec(__initial_guess(prob.u0, prob.p, t))), hcat, mesh) + end + + no_odes = length(u0_) + + # has_initial_guess = prob.u0 isa AbstractVector{<:AbstractArray} + # dt ≤ 0 && throw(ArgumentError("dt must be positive")) + # no_odes, n, u0 = if has_initial_guess + # length(first(prob.u0)), (length(prob.u0) - 1), reduce(hcat, prob.u0) + # else + # length(prob.u0), Int(cld((prob.tspan[2] - prob.tspan[1]), dt)), prob.u0 + # end + T = eltype(u0) - mesh = collect(range(prob.tspan[1], stop = prob.tspan[2], length = n + 1)) - opt = OptionsODE( - OPT_BVPCLASS => alg.bvpclass, OPT_COLLOCATIONPTS => alg.collocationpts, + # mesh = collect(range(prob.tspan[1], stop = prob.tspan[2], length = n + 1)) + opt = OptionsODE(OPT_BVPCLASS => alg.bvpclass, OPT_COLLOCATIONPTS => alg.collocationpts, OPT_MAXSTEPS => maxiters, OPT_DIAGNOSTICOUTPUT => alg.diagnostic_output, OPT_MAXSUBINTERVALS => alg.max_num_subintervals, OPT_RTOL => reltol) orders = ones(Int, no_odes) _tspan = [prob.tspan[1], prob.tspan[2]] iip = SciMLBase.isinplace(prob) - rhs(t, u, du) = + rhs = @closure (t, u, du) -> begin if iip prob.f(du, u, prob.p, t) else (du .= prob.f(u, prob.p, t)) end + end if prob.f.jac === nothing if iip - jac = function (df, u, p, t) + jac = (df, u, p, t) -> begin _du = similar(u) prob.f(_du, u, p, t) - _f = (du, u) -> prob.f(du, u, p, t) + _f = @closure (du, u) -> prob.f(du, u, p, t) ForwardDiff.jacobian!(df, _f, _du, u) + return end else - jac = function (df, u, p, t) + jac = (df, u, p, t) -> begin _du = prob.f(u, p, t) - _f = (du, u) -> (du .= prob.f(u, p, t)) + _f = @closure (du, u) -> (du .= prob.f(u, p, t)) ForwardDiff.jacobian!(df, _f, _du, u) + return end end else jac = prob.f.jac end - Drhs(t, u, df) = jac(df, u, prob.p, t) + Drhs = @closure (t, u, df) -> jac(df, u, prob.p, t) - #TODO: Fix bc and bcjac for multi-points BVP + bcresid_prototype, _ = BoundaryValueDiffEq.__get_bcresid_prototype( + prob.problem_type, prob, u0) - n_bc_a = length(first(prob.f.bcresid_prototype.x)) - n_bc_b = length(last(prob.f.bcresid_prototype.x)) + n_bc_a = length(first(bcresid_prototype)) + n_bc_b = length(last(bcresid_prototype)) zeta = vcat(fill(first(prob.tspan), n_bc_a), fill(last(prob.tspan), n_bc_b)) - bc = function (i, z, resid) + bc = @closure (i, z, resid) -> begin tmpa = copy(z) tmpb = copy(z) tmp_resid_a = zeros(T, n_bc_a) @@ -268,7 +284,7 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol end end - Dbc = function (i, z, dbc) + Dbc = @closure (i, z, dbc) -> begin for j in 1:n_bc_a if i == j dbc[i] = 1.0 @@ -287,7 +303,8 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol if retcode == 0 @warn "Collocation matrix is singular" elseif retcode == -1 - @warn "The expected no. of subintervals exceeds storage(try to increase `OPT_MAXSUBINTERVALS`)" + @warn "The expected no. of subintervals exceeds storage(try to increase \ + `OPT_MAXSUBINTERVALS`)" elseif retcode == -2 @warn "The nonlinear iteration has not converged" elseif retcode == -3 @@ -299,11 +316,10 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol destats = SciMLBase.DEStats( stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0, 0, 0) - return DiffEqBase.build_solution(prob, alg, mesh, - collect(Vector{eltype(evalsol)}, eachrow(evalsol)); + return DiffEqBase.build_solution( + prob, alg, mesh, collect(Vector{eltype(evalsol)}, eachrow(evalsol)); retcode = retcode > 0 ? ReturnCode.Success : ReturnCode.Failure, - stats = destats) + stats = destats, original = (sol, retcode, stats)) end -=# end diff --git a/src/algorithms.jl b/src/algorithms.jl index edfdade2..4caf2aaf 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -295,6 +295,10 @@ struct BVPSOL{O} <: BoundaryValueDiffEqAlgorithm end end +function BVPSOL(; bvpclass = 2, sol_method = 0, odesolver = nothing) + return BVPSOL(bvpclass, sol_method, odesolver) +end + """ COLNEW(; bvpclass = 1, collocationpts = 7, diagnostic_output = 1, max_num_subintervals = 3000) diff --git a/test/misc/odeinterface_wrapper.jl b/test/wrappers/odeinterface_tests.jl similarity index 59% rename from test/misc/odeinterface_wrapper.jl rename to test/wrappers/odeinterface_tests.jl index 78bc00c7..2e253737 100644 --- a/test/misc/odeinterface_wrapper.jl +++ b/test/wrappers/odeinterface_tests.jl @@ -1,5 +1,6 @@ -using Test, BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random, OrdinaryDiffEq, - RecursiveArrayTools +@testsetup module ODEInterfaceWrapperTestSetup + +using BoundaryValueDiffEq, LinearAlgebra, ODEInterface, Random, RecursiveArrayTools # Adaptation of https://github.com/luchr/ODEInterface.jl/blob/958b6023d1dabf775033d0b89c5401b33100bca3/examples/BasicExamples/ex7.jl function ex7_f!(du, u, p, t) @@ -27,8 +28,16 @@ tspan = (-π / 2, π / 2) tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1))) -@testset "BVPM2" begin - @info "Testing BVPM2" +# 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 + +end + +@testitem "BVPM2" setup=[ODEInterfaceWrapperTestSetup] begin + using ODEInterface, RecursiveArrayTools sol_bvpm2 = solve(tpprob, BVPM2(); dt = π / 20) @test SciMLBase.successful_retcode(sol_bvpm2) @@ -38,70 +47,61 @@ tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), u0, tspan, @test norm(resid_f, Inf) < 1e-6 end -# Just generate a solution for bvpsol -sol_ms = solve(tpprob, MultipleShooting(10, DP5(), NewtonRaphson()); - dt = π / 20, abstol = 1e-5, maxiters = 1000, - odesolve_kwargs = (; adaptive = false, dt = 0.01, abstol = 1e-6, maxiters = 1000)) - # Just test that it runs. BVPSOL only works with linearly separable BCs. -@testset "BVPSOL" begin - @info "Testing BVPSOL" - - @info "BVPSOL with Vector{<:AbstractArray}" +@testitem "BVPSOL" setup=[ODEInterfaceWrapperTestSetup] begin + using ODEInterface, RecursiveArrayTools 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))) - # Just test that it runs. BVPSOL only works with linearly separable BCs. sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - @info "BVPSOL with VectorOfArray" + @test sol_bvpsol isa SciMLBase.ODESolution initial_u0 = VectorOfArray([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))) - # Just test that it runs. BVPSOL only works with linearly separable BCs. sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - @info "BVPSOL with DiffEqArray" + @test sol_bvpsol isa SciMLBase.ODESolution ts = collect(tspan[1]:(π / 20):tspan[2]) initial_u0 = DiffEqArray([sol_ms(t) .+ rand() for t in ts], ts) tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan, p; bcresid_prototype = (zeros(1), zeros(1))) - sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) - @info "BVPSOL with initial guess function" + @test sol_bvpsol isa SciMLBase.ODESolution initial_u0 = (p, t) -> sol_ms(t) .+ rand() - # FIXME: Upstream fix - # tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, tspan, p; - # bcresid_prototype = (zeros(1), zeros(1))) - # sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) -end - -#= -@info "COLNEW" + tpprob = TwoPointBVProblem(ex7_f!, (ex7_2pbc1!, ex7_2pbc2!), initial_u0, + tspan, p; bcresid_prototype = (zeros(1), zeros(1))) + sol_bvpsol = solve(tpprob, BVPSOL(); dt = π / 20) -function f!(du, u, p, t) - du[1] = u[2] - du[2] = u[1] -end -function bca!(resid_a, u_a, p) - resid_a[1] = u_a[1] - 1 -end -function bcb!(resid_b, u_b, p) - resid_b[1] = u_b[1] + @test sol_bvpsol isa SciMLBase.ODESolution end -fun = BVPFunction( - f!, (bca!, bcb!), bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true)) -tspan = (0.0, 1.0) - -prob = TwoPointBVProblem(fun, [1.0, 0.0], tspan) -sol_colnew = solve(prob, COLNEW(), dt = 0.01) -@test SciMLBase.successful_retcode(sol_colnew) -=# +@testitem "COLNEW" setup=[ODEInterfaceWrapperTestSetup] begin + using ODEInterface, RecursiveArrayTools + + function f!(du, u, p, t) + du[1] = u[2] + du[2] = u[1] + end + function bca!(resid_a, u_a, p) + resid_a[1] = u_a[1] - 1 + end + function bcb!(resid_b, u_b, p) + resid_b[1] = u_b[1] + end + + fun = BVPFunction( + f!, (bca!, bcb!), bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true)) + tspan = (0.0, 1.0) + + prob = TwoPointBVProblem(fun, [1.0, 0.0], tspan) + sol_colnew = solve(prob, COLNEW(), dt = 0.01) + @test SciMLBase.successful_retcode(sol_colnew) +end