diff --git a/benchmark/simple_pendulum.jl b/benchmark/simple_pendulum.jl index a2d66bb8..8e9df3d0 100644 --- a/benchmark/simple_pendulum.jl +++ b/benchmark/simple_pendulum.jl @@ -43,32 +43,42 @@ function create_simple_pendulum_benchmark() suite["OOP"] = oop_suite if @isdefined(MultipleShooting) - iip_suite["MultipleShooting(100, Tsit5; grid_coarsening = true)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_iip, + iip_suite["MultipleShooting(100, Tsit5; grid_coarsening = true)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_iip, $MultipleShooting(100, Tsit5())) - iip_suite["MultipleShooting(100, Tsit5; grid_coarsening = false)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_iip, + iip_suite["MultipleShooting(100, Tsit5; grid_coarsening = false)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_iip, $MultipleShooting(100, Tsit5(); grid_coarsening = false)) - iip_suite["MultipleShooting(10, Tsit5; grid_coarsening = true)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_iip, + iip_suite["MultipleShooting(10, Tsit5; grid_coarsening = true)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_iip, $MultipleShooting(10, Tsit5())) - iip_suite["MultipleShooting(10, Tsit5; grid_coarsening = false)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_iip, + iip_suite["MultipleShooting(10, Tsit5; grid_coarsening = false)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_iip, $MultipleShooting(10, Tsit5(); grid_coarsening = false)) end if @isdefined(Shooting) - iip_suite["Shooting(Tsit5())"] = @benchmarkable solve($SimplePendulumBenchmark.prob_iip, + iip_suite["Shooting(Tsit5())"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_iip, $Shooting(Tsit5())) end if @isdefined(MultipleShooting) - oop_suite["MultipleShooting(100, Tsit5; grid_coarsening = true)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_oop, + oop_suite["MultipleShooting(100, Tsit5; grid_coarsening = true)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_oop, $MultipleShooting(100, Tsit5())) - oop_suite["MultipleShooting(100, Tsit5; grid_coarsening = false)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_oop, + oop_suite["MultipleShooting(100, Tsit5; grid_coarsening = false)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_oop, $MultipleShooting(100, Tsit5(); grid_coarsening = false)) - oop_suite["MultipleShooting(10, Tsit5; grid_coarsening = true)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_oop, + oop_suite["MultipleShooting(10, Tsit5; grid_coarsening = true)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_oop, $MultipleShooting(10, Tsit5())) - oop_suite["MultipleShooting(10, Tsit5; grid_coarsening = false)"] = @benchmarkable solve($SimplePendulumBenchmark.prob_oop, + oop_suite["MultipleShooting(10, Tsit5; grid_coarsening = false)"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_oop, $MultipleShooting(10, Tsit5(); grid_coarsening = false)) end if @isdefined(Shooting) - oop_suite["Shooting(Tsit5())"] = @benchmarkable solve($SimplePendulumBenchmark.prob_oop, + oop_suite["Shooting(Tsit5())"] = @benchmarkable solve( + $SimplePendulumBenchmark.prob_oop, $Shooting(Tsit5())) end diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqODEInterfaceExt.jl index 960d3b01..562d560f 100644 --- a/ext/BoundaryValueDiffEqODEInterfaceExt.jl +++ b/ext/BoundaryValueDiffEqODEInterfaceExt.jl @@ -3,15 +3,18 @@ module BoundaryValueDiffEqODEInterfaceExt using SciMLBase, BoundaryValueDiffEq, ODEInterface import SciMLBase: __solve import ODEInterface: OptionsODE, OPT_ATOL, OPT_RTOL, OPT_METHODCHOICE, OPT_DIAGNOSTICOUTPUT, - OPT_ERRORCONTROL, OPT_SINGULARTERM, OPT_MAXSTEPS, OPT_BVPCLASS, OPT_SOLMETHOD, - OPT_RHS_CALLMODE, OPT_COLLOCATIONPTS, OPT_MAXSUBINTERVALS, RHS_CALL_INSITU, evalSolution + OPT_ERRORCONTROL, OPT_SINGULARTERM, OPT_MAXSTEPS, OPT_BVPCLASS, + OPT_SOLMETHOD, + OPT_RHS_CALLMODE, OPT_COLLOCATIONPTS, OPT_MAXSUBINTERVALS, + RHS_CALL_INSITU, evalSolution import ODEInterface: Bvpm2, bvpm2_init, bvpm2_solve, bvpm2_destroy, bvpm2_get_x import ODEInterface: bvpsol import ODEInterface: colnew import ForwardDiff -function _test_bvpm2_bvpsol_colnew_problem_criteria(_, ::SciMLBase.StandardBVProblem, alg::Symbol) +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) @@ -54,7 +57,8 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPM2; dt = 0.0, reltol = 1e-3, sol, retcode, stats = bvpm2_solve(initial_guess, bvp2m_f, bvp2m_bc, opt) retcode = retcode ≥ 0 ? ReturnCode.Success : ReturnCode.Failure - destats = SciMLBase.DEStats(stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0) + destats = SciMLBase.DEStats( + stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0) x_mesh = bvpm2_get_x(sol) evalsol = evalSolution(sol, x_mesh) @@ -72,7 +76,7 @@ end #------- function SciMLBase.__solve(prob::BVProblem, alg::BVPSOL; maxiters = 1000, reltol = 1e-3, dt = 0.0, verbose = true, kwargs...) - _test_bvpm2_bvpsol_colnew_problem_criteria(prob, prob.problem_type, :BVPSOL) + _test_bvpm2_bvpsol_colnew_problem_criteria(prob, prob.problem_type, :BVPSOL) @assert isa(prob.p, SciMLBase.NullParameters) "BVPSOL only supports NullParameters!" @assert isa(prob.u0, AbstractVector{<:AbstractArray}) "BVPSOL requires a vector of initial guesses!" n, u0 = (length(prob.u0) - 1), reduce(hcat, prob.u0) @@ -125,7 +129,8 @@ end #------- # COLNEW #------- -function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol=1e-4, dt = 0.0, verbose = true, kwargs...) +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")) @@ -136,18 +141,20 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol 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, + 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) + 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) = if iip - prob.f(du, u, prob.p, t) - else - (du .= prob.f(u, prob.p, t)) - end + rhs(t, u, du) = + if iip + prob.f(du, u, prob.p, t) + else + (du .= prob.f(u, prob.p, t)) + end if prob.f.jac === nothing if iip @@ -169,25 +176,25 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol end Drhs(t, u, df) = jac(df, u, prob.p, t) - #TODO: Fix bc and bcjac for multi-points BVP n_bc_a = length(first(prob.f.bcresid_prototype.x)) n_bc_b = length(last(prob.f.bcresid_prototype.x)) zeta = vcat(fill(first(prob.tspan), n_bc_a), fill(last(prob.tspan), n_bc_b)) bc = function (i, z, resid) - tmpa = copy(z); tmpb = copy(z) + tmpa = copy(z) + tmpb = copy(z) tmp_resid_a = zeros(T, n_bc_a) tmp_resid_b = zeros(T, n_bc_b) prob.f.bc[1](tmp_resid_a, tmpa, prob.p) prob.f.bc[2](tmp_resid_b, tmpb, prob.p) - for j=1:n_bc_a + for j in 1:n_bc_a if i == j resid[1] = tmp_resid_a[j] end end - for j=1:n_bc_b + for j in 1:n_bc_b if i == (j + n_bc_a) resid[1] = tmp_resid_b[j] end @@ -195,12 +202,12 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol end Dbc = function (i, z, dbc) - for j=1:n_bc_a + for j in 1:n_bc_a if i == j dbc[i] = 1.0 end end - for j=1:n_bc_b + for j in 1:n_bc_b if i == (j + n_bc_a) dbc[i] = 1.0 end @@ -222,7 +229,8 @@ function SciMLBase.__solve(prob::BVProblem, alg::COLNEW; maxiters = 1000, reltol end evalsol = evalSolution(sol, mesh) - destats = SciMLBase.DEStats(stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0) + destats = SciMLBase.DEStats( + stats["no_rhs_calls"], 0, 0, 0, stats["no_jac_calls"], 0, 0, 0, 0, 0, 0) return DiffEqBase.build_solution(prob, alg, mesh, collect(Vector{eltype(evalsol)}, eachrow(evalsol)); diff --git a/ext/BoundaryValueDiffEqOrdinaryDiffEqExt.jl b/ext/BoundaryValueDiffEqOrdinaryDiffEqExt.jl index 9adb76ad..b03d47b5 100644 --- a/ext/BoundaryValueDiffEqOrdinaryDiffEqExt.jl +++ b/ext/BoundaryValueDiffEqOrdinaryDiffEqExt.jl @@ -38,7 +38,7 @@ end BVProblem(f1!, bc1!, u0, tspan), BVProblem(f1, bc1, u0, tspan), TwoPointBVProblem(f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype), - TwoPointBVProblem(f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype), + TwoPointBVProblem(f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype) ] algs = [] @@ -101,7 +101,7 @@ end TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan; bcresid_prototype = bcresid_prototype2), TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan; - bcresid_prototype = bcresid_prototype2), + bcresid_prototype = bcresid_prototype2) ] algs = [] @@ -112,7 +112,7 @@ end Shooting(Tsit5(); nlsolve = LevenbergMarquardt(), jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))), Shooting(Tsit5(); nlsolve = GaussNewton(), - jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))), + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))) ]) end @@ -129,7 +129,7 @@ end nlsolve = GaussNewton(; autodiff = AutoForwardDiff(chunksize = 2)), jac_alg = BVPJacobianAlgorithm(; bc_diffmode = AutoForwardDiff(; chunksize = 2), - nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2))), + nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2))) ]) end diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl index 9451db13..9a5cdb91 100644 --- a/src/BoundaryValueDiffEq.jl +++ b/src/BoundaryValueDiffEq.jl @@ -4,15 +4,16 @@ import PrecompileTools: @compile_workload, @setup_workload, @recompile_invalidat @recompile_invalidations begin using ADTypes, Adapt, DiffEqBase, ForwardDiff, LinearAlgebra, NonlinearSolve, - PreallocationTools, Preferences, RecursiveArrayTools, Reexport, SciMLBase, Setfield, - SparseDiffTools, Tricks + PreallocationTools, Preferences, RecursiveArrayTools, Reexport, SciMLBase, + Setfield, + SparseDiffTools, Tricks # Special Matrix Types using BandedMatrices, FastAlmostBandedMatrices, SparseArrays import ADTypes: AbstractADType import ArrayInterface: matrix_colors, - parameterless_type, undefmatrix, fast_scalar_indexing + parameterless_type, undefmatrix, fast_scalar_indexing import ConcreteStructs: @concrete import DiffEqBase: solve import ForwardDiff: pickchunksize @@ -74,7 +75,7 @@ end BVProblem(f1!, bc1!, u0, tspan), BVProblem(f1, bc1, u0, tspan), TwoPointBVProblem(f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype), - TwoPointBVProblem(f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype), + TwoPointBVProblem(f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype) ] algs = [] @@ -129,7 +130,7 @@ end TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan; bcresid_prototype = bcresid_prototype2), TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan; - bcresid_prototype = bcresid_prototype2), + bcresid_prototype = bcresid_prototype2) ] jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) @@ -144,7 +145,7 @@ end [ MIRK2(; jac_alg, nlsolve), MIRK3(; jac_alg, nlsolve), MIRK4(; jac_alg, nlsolve), MIRK5(; jac_alg, nlsolve), - MIRK6(; jac_alg, nlsolve), + MIRK6(; jac_alg, nlsolve) ]) end end diff --git a/src/algorithms.jl b/src/algorithms.jl index 94bd9895..00cb80f8 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -55,9 +55,11 @@ function __propagate_nlsolve_ad_to_jac_alg(nlsolve::N) where {N} ad = hasfield(N, :ad) ? nlsolve.ad : nothing ad === nothing && return BVPJacobianAlgorithm() - Base.depwarn("Setting autodiff to the nonlinear solver in Shooting has been deprecated \ - and will have no effect from the next major release. Update to use \ - `BVPJacobianAlgorithm` directly", :Shooting) + Base.depwarn( + "Setting autodiff to the nonlinear solver in Shooting has been deprecated \ + and will have no effect from the next major release. Update to use \ + `BVPJacobianAlgorithm` directly", + :Shooting) return BVPJacobianAlgorithm(ad) end diff --git a/src/mirk_tableaus.jl b/src/mirk_tableaus.jl index 097bd60e..d9038a7e 100644 --- a/src/mirk_tableaus.jl +++ b/src/mirk_tableaus.jl @@ -31,7 +31,7 @@ function constructMIRK3(::Type{T}) where {T} v = [0, 4 // 9] b = [1 // 4, 3 // 4] x = [0 0 - 2//9 0] + 2//9 0] # Interpolant tableau s_star = 3 @@ -52,8 +52,8 @@ function constructMIRK4(::Type{T}) where {T} v = [0, 1, 1 // 2, 27 // 32] b = [1 // 6, 1 // 6, 2 // 3, 0] x = [0 0 0 0 - 0 0 0 0 - 1//8 -1//8 0 0] + 0 0 0 0 + 1//8 -1//8 0 0] # Interpolant tableau s_star = 4 @@ -74,16 +74,16 @@ function constructMIRK5(::Type{T}) where {T} v = [0, 1, 27 // 32, 837 // 1250] b = [5 // 54, 1 // 14, 32 // 81, 250 // 567] x = [0 0 0 0 - 0 0 0 0 - 3//64 -9//64 0 0 - 21//1000 63//5000 -252//625 0] + 0 0 0 0 + 3//64 -9//64 0 0 + 21//1000 63//5000 -252//625 0] # Interpolant tableau s_star = 6 c_star = [4 // 5, 13 // 23] v_star = [4 // 5, 13 // 23] x_star = [14//1125 -74//875 -128//3375 104//945 0 0 - 1//2 4508233//1958887 48720832//2518569 -27646420//17629983 -11517095//559682 0] + 1//2 4508233//1958887 48720832//2518569 -27646420//17629983 -11517095//559682 0] τ_star = 0.3 TU = MIRKTableau(s, T.(c), T.(v), T.(b), T.(x)) @@ -98,19 +98,19 @@ function constructMIRK6(::Type{T}) where {T} v = [0, 1, 5 // 32, 27 // 32, 1 // 2] b = [7 // 90, 7 // 90, 16 // 45, 16 // 45, 2 // 15, 0, 0, 0, 0] x = [0 0 0 0 0 - 0 0 0 0 0 - 9//64 -3//64 0 0 0 - 3//64 -9//64 0 0 0 - -5//24 5//24 2//3 -2//3 0] + 0 0 0 0 0 + 9//64 -3//64 0 0 0 + 3//64 -9//64 0 0 0 + -5//24 5//24 2//3 -2//3 0] # Interpolant tableau s_star = 9 c_star = [7 // 16, 3 // 8, 9 // 16, 1 // 8] v_star = [7 // 16, 3 // 8, 9 // 16, 1 // 8] x_star = [1547//32768 -1225//32768 749//4096 -287//2048 -861//16384 0 0 0 0 - 83//1536 -13//384 283//1536 -167//1536 -49//512 0 0 0 0 - 1225//32768 -1547//32768 287//2048 -749//4096 861//16384 0 0 0 0 - 233//3456 -19//1152 0 0 0 -5//72 7//72 -17//216 0] + 83//1536 -13//384 283//1536 -167//1536 -49//512 0 0 0 0 + 1225//32768 -1547//32768 287//2048 -749//4096 861//16384 0 0 0 0 + 233//3456 -19//1152 0 0 0 -5//72 7//72 -17//216 0] τ_star = 0.7156 TU = MIRKTableau(s, T.(c), T.(v), T.(b), T.(x)) diff --git a/src/solve/mirk.jl b/src/solve/mirk.jl index ae9eaed5..37721222 100644 --- a/src/solve/mirk.jl +++ b/src/solve/mirk.jl @@ -224,7 +224,8 @@ function __mirk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual, return nothing end -function __mirk_loss!(resid, u, p, y, pt::TwoPointBVProblem, bc!::Tuple{BC1, BC2}, residual, +function __mirk_loss!( + resid, u, p, y, pt::TwoPointBVProblem, bc!::Tuple{BC1, BC2}, residual, mesh, cache) where {BC1, BC2} y_ = recursive_unflatten!(y, u) resids = [get_tmp(r, u) for r in residual] @@ -347,7 +348,8 @@ function __mirk_mpoint_jacobian!(J, _, x, bc_diffmode, nonbc_diffmode, bc_diffca return nothing end -function __mirk_mpoint_jacobian!(J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode, +function __mirk_mpoint_jacobian!( + J::AlmostBandedMatrix, J_c, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache, loss_bc::BC, loss_collocation::C, resid_bc, resid_collocation, L::Int) where {BC, C} J_bc = fillpart(J) diff --git a/src/solve/multiple_shooting.jl b/src/solve/multiple_shooting.jl index f26f9924..565e56fe 100644 --- a/src/solve/multiple_shooting.jl +++ b/src/solve/multiple_shooting.jl @@ -54,12 +54,14 @@ function __solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwargs = (;), end if prob.problem_type isa TwoPointBVProblem - __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, + __solve_nlproblem!( + prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, cur_nshoot, M, N, resida_len, residb_len, solve_internal_odes!, bc[1], bc[2], prob, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...) else - __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, + __solve_nlproblem!( + prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes, cur_nshoot, M, N, prod(resid_size), solve_internal_odes!, bc, prob, f, u0_size, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...) @@ -171,7 +173,8 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_ # Define the functions now ode_fn = (du, u) -> solve_internal_odes!(du, u, prob.p, cur_nshoot, nodes, ode_cache_ode_jac_fn) - bc_fn = (du, u) -> __multiple_shooting_mpoint_loss_bc!(du, u, prob.p, cur_nshoot, nodes, + bc_fn = (du, u) -> __multiple_shooting_mpoint_loss_bc!( + du, u, prob.p, cur_nshoot, nodes, prob, solve_internal_odes!, N, f, bc, u0_size, prob.tspan, alg.ode_alg, u0, ode_cache_bc_jac_fn) @@ -257,7 +260,8 @@ function __multiple_shooting_solve_internal_odes!(resid_nodes, us, cur_nshoots:: Threads.@threads for idx in 1:length(data_partition) cache = odecache[idx] for i in data_partition[idx] - SciMLBase.reinit!(cache, reshape(@view(us[((i - 1) * N + 1):(i * N)]), u0_size); + SciMLBase.reinit!( + cache, reshape(@view(us[((i - 1) * N + 1):(i * N)]), u0_size); t0 = nodes[i], tf = nodes[i + 1]) sol = solve!(cache) us_[i] = deepcopy(sol.u) @@ -355,7 +359,8 @@ end end # Problem has initial guess -@views function __multiple_shooting_initialize!(nodes, prob, alg, ::Val{true}, nshoots::Int, +@views function __multiple_shooting_initialize!( + nodes, prob, alg, ::Val{true}, nshoots::Int, odecache; kwargs...) @unpack u0, tspan = prob diff --git a/src/utils.jl b/src/utils.jl index 5dfa6123..8d872bc2 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -154,7 +154,8 @@ function __initial_guess(f::F, p::P, t::T) where {F, P, T} elseif static_hasmethod(f, Tuple{T}) Base.depwarn("initial guess function must take 2 inputs `(p, t)` instead of just \ `t`. The single argument version has been deprecated and will be \ - removed in the next major release of SciMLBase.", :__initial_guess) + removed in the next major release of SciMLBase.", + :__initial_guess) return f(t) else throw(ArgumentError("`initial_guess` must be a function of the form `f(p, t)`")) diff --git a/test/mirk/mirk_convergence_tests.jl b/test/mirk/mirk_convergence_tests.jl index cdfd879f..c93b11b6 100644 --- a/test/mirk/mirk_convergence_tests.jl +++ b/test/mirk/mirk_convergence_tests.jl @@ -43,12 +43,12 @@ odef1 = ODEFunction(f1, analytic = (u0, p, t) -> [5 - t, -1]) odef2! = ODEFunction(f2!, analytic = (u0, p, t) -> [ 5 * (cos(t) - cot(5) * sin(t)), - 5 * (-cos(t) * cot(5) - sin(t)), + 5 * (-cos(t) * cot(5) - sin(t)) ]) odef2 = ODEFunction(f2, analytic = (u0, p, t) -> [ 5 * (cos(t) - cot(5) * sin(t)), - 5 * (-cos(t) * cot(5) - sin(t)), + 5 * (-cos(t) * cot(5) - sin(t)) ]) bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) @@ -68,7 +68,7 @@ probArr = [ TwoPointBVProblem(odef2!, (boundary_two_point_a!, boundary_two_point_b!), u0, tspan; bcresid_prototype), TwoPointBVProblem(odef2, (boundary_two_point_a, boundary_two_point_b), u0, tspan; - bcresid_prototype), + bcresid_prototype) ]; testTol = 0.2 diff --git a/test/mirk/nonlinear_least_squares.jl b/test/mirk/nonlinear_least_squares.jl index edbbee79..d3a69aa5 100644 --- a/test/mirk/nonlinear_least_squares.jl +++ b/test/mirk/nonlinear_least_squares.jl @@ -1,7 +1,8 @@ using BoundaryValueDiffEq, LinearAlgebra, Test @testset "Overconstrained BVP" begin - SOLVERS = [mirk(; nlsolve) for mirk in (MIRK4, MIRK5, MIRK6), + SOLVERS = [mirk(; nlsolve) + for mirk in (MIRK4, MIRK5, MIRK6), nlsolve in (LevenbergMarquardt(), GaussNewton(), nothing)] # OOP MP-BVP @@ -53,8 +54,11 @@ using BoundaryValueDiffEq, LinearAlgebra, Test bc1a(ua, p) = [ua[1]] bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] - bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + bvp3 = TwoPointBVProblem( + BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), + u0, + tspan) for solver in SOLVERS @time sol = solve(bvp3, solver; verbose = false, dt = 1.0) @@ -65,8 +69,11 @@ using BoundaryValueDiffEq, LinearAlgebra, Test bc1a!(resid, ua, p) = (resid[1] = ua[1]) bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) - bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + bvp4 = TwoPointBVProblem( + BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), + u0, + tspan) for solver in SOLVERS @time sol = solve(bvp3, solver; verbose = false, dt = 1.0, abstol = 1e-3, diff --git a/test/misc/non_vector_inputs.jl b/test/misc/non_vector_inputs.jl index f3a98f9d..101ea232 100644 --- a/test/misc/non_vector_inputs.jl +++ b/test/misc/non_vector_inputs.jl @@ -40,7 +40,7 @@ probs = [ TwoPointBVProblem(f1!, (boundary_a!, boundary_b!), u0, tspan; bcresid_prototype = (Array{Float64}(undef, 1, 1), Array{Float64}(undef, 1, 1))), BVProblem(f1, boundary, u0, tspan), - TwoPointBVProblem(f1, (boundary_a, boundary_b), u0, tspan), + TwoPointBVProblem(f1, (boundary_a, boundary_b), u0, tspan) ]; @testset "Affineness" begin diff --git a/test/misc/odeinterface_wrapper.jl b/test/misc/odeinterface_wrapper.jl index 2dfa85ff..ead00b84 100644 --- a/test/misc/odeinterface_wrapper.jl +++ b/test/misc/odeinterface_wrapper.jl @@ -65,7 +65,8 @@ 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)) +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) diff --git a/test/shooting/nonlinear_least_squares.jl b/test/shooting/nonlinear_least_squares.jl index c45d963b..d1142d46 100644 --- a/test/shooting/nonlinear_least_squares.jl +++ b/test/shooting/nonlinear_least_squares.jl @@ -64,8 +64,11 @@ using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test bc1a(ua, p) = [ua[1]] bc1b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] - bvp3 = TwoPointBVProblem(BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + bvp3 = TwoPointBVProblem( + BVPFunction{false}(f1, (bc1a, bc1b); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), + u0, + tspan) for solver in SOLVERS sol = @time solve(bvp3, solver; verbose = false) @@ -76,8 +79,11 @@ using BoundaryValueDiffEq, LinearAlgebra, LinearSolve, OrdinaryDiffEq, Test bc1a!(resid, ua, p) = (resid[1] = ua[1]) bc1b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) - bvp4 = TwoPointBVProblem(BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), - bcresid_prototype = (zeros(1), zeros(2))), u0, tspan) + bvp4 = TwoPointBVProblem( + BVPFunction{true}(f1!, (bc1a!, bc1b!); twopoint = Val(true), + bcresid_prototype = (zeros(1), zeros(2))), + u0, + tspan) for solver in SOLVERS sol = @time solve(bvp4, solver; verbose = false) diff --git a/test/shooting/orbital.jl b/test/shooting/orbital.jl index aa251dfb..ff4337d9 100644 --- a/test/shooting/orbital.jl +++ b/test/shooting/orbital.jl @@ -7,7 +7,7 @@ y0 = [ -5.3500183933132319E+06, -5528.612564911408, 1216.8442360202787, - 4845.114446429901, + 4845.114446429901 ] init_val = [ -4.7763169762853989E+06, @@ -15,7 +15,7 @@ init_val = [ -5.3500183933132319E+06, 7.0526926403748598E+06, -7.9650476230388973E+05, - -1.1911128863666430E+06, + -1.1911128863666430E+06 ] J2 = 1.08262668E-3 req = 6378137 @@ -66,7 +66,8 @@ resid_f_2p = (Array{Float64, 1}(undef, 3), Array{Float64, 1}(undef, 3)) ### Now use the BVP solver to get closer bvp = BVProblem(orbital!, cur_bc!, y0, tspan) -for autodiff in (AutoForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:central)), +for autodiff in ( + AutoForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:central)), AutoSparseForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:forward)), AutoSparseFiniteDiff()) nlsolve = TrustRegion(; autodiff) @@ -88,7 +89,8 @@ end ### Using the TwoPoint BVP Structure bvp = TwoPointBVProblem(orbital!, (cur_bc_2point_a!, cur_bc_2point_b!), y0, tspan; bcresid_prototype = (Array{Float64}(undef, 3), Array{Float64}(undef, 3))) -for autodiff in (AutoForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:central)), +for autodiff in ( + AutoForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:central)), AutoSparseForwardDiff(; chunksize = 6), AutoFiniteDiff(; fdtype = Val(:forward)), AutoSparseFiniteDiff()) nlsolve = TrustRegion(; autodiff)