diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 16bbaea19..7e8d929a9 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,10 +13,14 @@ concurrency: cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} jobs: test: - runs-on: ubuntu-latest + runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: + os: + - ubuntu-latest + - macos-latest + - windows-latest version: - '1' steps: diff --git a/README.md b/README.md index 5cef548d4..2a7f13c05 100644 --- a/README.md +++ b/README.md @@ -46,11 +46,11 @@ For the list of available solvers, please refer to the [DifferentialEquations.jl Precompilation can be controlled via `Preferences.jl` - `PrecompileMIRK` -- Precompile the MIRK2 - MIRK6 algorithms (default: `true`). - - `PrecompileShooting` -- Precompile the single shooting algorithms (default: `false`). This is triggered when `OrdinaryDiffEq` is loaded. - - `PrecompileMultipleShooting` -- Precompile the multiple shooting algorithms (default: `false`). This is triggered when `OrdinaryDiffEq` is loaded. - - `PrecompileMIRKNLLS` -- Precompile the MIRK2 - MIRK6 algorithms for under-determined and over-determined BVPs (default: `false`). - - `PrecompileShootingNLLS` -- Precompile the single shooting algorithms for under-determined and over-determined BVPs (default: `false`). This is triggered when `OrdinaryDiffEq` is loaded. - - `PrecompileMultipleShootingNLLS` -- Precompile the multiple shooting algorithms for under-determined and over-determined BVPs (default: `false` ). This is triggered when `OrdinaryDiffEq` is loaded. + - `PrecompileShooting` -- Precompile the single shooting algorithms (default: `true`). + - `PrecompileMultipleShooting` -- Precompile the multiple shooting algorithms (default: `true`). + - `PrecompileMIRKNLLS` -- Precompile the MIRK2 - MIRK6 algorithms for under-determined and over-determined BVPs (default: `true`). + - `PrecompileShootingNLLS` -- Precompile the single shooting algorithms for under-determined and over-determined BVPs (default: `true`). + - `PrecompileMultipleShootingNLLS` -- Precompile the multiple shooting algorithms for under-determined and over-determined BVPs (default: `true` ). To set these preferences before loading the package, do the following (replacing `PrecompileShooting` with the preference you want to set, or pass in multiple pairs to set them together): diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqODEInterfaceExt.jl index 06fb822df..a601bef26 100644 --- a/ext/BoundaryValueDiffEqODEInterfaceExt.jl +++ b/ext/BoundaryValueDiffEqODEInterfaceExt.jl @@ -93,7 +93,7 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPM2; dt = 0.0, reltol = 1e-3, evalsol = evalSolution(sol, x_mesh) ivpsol = SciMLBase.build_solution(prob, alg, x_mesh, map(x -> reshape(convert(Vector{eltype(evalsol)}, x), u0_size), eachcol(evalsol)); - retcode, stats = destats) + retcode, stats = destats, original = (sol, retcode, stats)) bvpm2_destroy(obj) bvpm2_destroy(sol) @@ -187,7 +187,8 @@ 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) + retcode = retcode ≥ 0 ? ReturnCode.Success : ReturnCode.Failure, stats, + original = (sol_t, sol_x, retcode, stats)) return SciMLBase.build_solution(prob, ivpsol, nothing) end diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl index 2afa73dcb..fde4f655c 100644 --- a/src/BoundaryValueDiffEq.jl +++ b/src/BoundaryValueDiffEq.jl @@ -53,246 +53,221 @@ function __solve(prob::BVProblem, alg::BoundaryValueDiffEqAlgorithm, args...; kw end @setup_workload begin - # function f1!(du, u, p, t) - # du[1] = u[2] - # du[2] = 0 - # end - # f1(u, p, t) = [u[2], 0] - - # function bc1!(residual, u, p, t) - # residual[1] = u[1][1] - 5 - # residual[2] = u[end][1] - # end - # bc1(u, p, t) = [u[1][1] - 5, u[end][1]] - - # bc1_a!(residual, ua, p) = (residual[1] = ua[1] - 5) - # bc1_b!(residual, ub, p) = (residual[1] = ub[1]) - - # bc1_a(ua, p) = [ua[1] - 5] - # bc1_b(ub, p) = [ub[1]] - - # tspan = (0.0, 5.0) - # u0 = [5.0, -3.5] - # bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) - - # probs = [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)] - - # algs = [] - - # jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) - - # if Preferences.@load_preference("PrecompileMIRK", true) - # append!(algs, - # [MIRK2(; jac_alg), MIRK3(; jac_alg), MIRK4(; jac_alg), - # MIRK5(; jac_alg), MIRK6(; jac_alg)]) - # end - - # @compile_workload begin - # for prob in probs, alg in algs - # solve(prob, alg; dt = 0.2) - # end - # end - - # function f1_nlls!(du, u, p, t) - # du[1] = u[2] - # du[2] = -u[1] - # end - - # f1_nlls(u, p, t) = [u[2], -u[1]] - - # function bc1_nlls!(resid, sol, p, t) - # solₜ₁ = sol[1] - # solₜ₂ = sol[end] - # resid[1] = solₜ₁[1] - # resid[2] = solₜ₂[1] - 1 - # resid[3] = solₜ₂[2] + 1.729109 - # return nothing - # end - # bc1_nlls(sol, p, t) = [sol[1][1], sol[end][1] - 1, sol[end][2] + 1.729109] - - # bc1_nlls_a!(resid, ua, p) = (resid[1] = ua[1]) - # bc1_nlls_b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) - - # bc1_nlls_a(ua, p) = [ua[1]] - # bc1_nlls_b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] - - # tspan = (0.0, 100.0) - # u0 = [0.0, 1.0] - # bcresid_prototype1 = Array{Float64}(undef, 3) - # bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2)) - - # probs = [ - # BVProblem(BVPFunction(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1), - # u0, tspan), - # BVProblem(BVPFunction(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1), - # u0, tspan), - # 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)] - - # jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) - - # nlsolvers = [LevenbergMarquardt(), GaussNewton()] - - # algs = [] - - # if Preferences.@load_preference("PrecompileMIRKNLLS", false) - # for nlsolve in nlsolvers - # append!(algs, - # [MIRK2(; jac_alg, nlsolve), MIRK3(; jac_alg, nlsolve), - # MIRK4(; jac_alg, nlsolve), MIRK5(; jac_alg, nlsolve), - # MIRK6(; jac_alg, nlsolve)]) - # end - # end - - # @compile_workload begin - # for prob in probs, alg in algs - # solve(prob, alg; dt = 0.2) - # end - # end - - # function f1!(du, u, p, t) - # du[1] = u[2] - # du[2] = 0 - # end - # f1(u, p, t) = [u[2], 0] - - # function bc1!(residual, u, p, t) - # residual[1] = u(0.0)[1] - 5 - # residual[2] = u(5.0)[1] - # end - # bc1(u, p, t) = [u(0.0)[1] - 5, u(5.0)[1]] - - # bc1_a!(residual, ua, p) = (residual[1] = ua[1] - 5) - # bc1_b!(residual, ub, p) = (residual[1] = ub[1]) - - # bc1_a(ua, p) = [ua[1] - 5] - # bc1_b(ub, p) = [ub[1]] - - # tspan = (0.0, 5.0) - # u0 = [5.0, -3.5] - # bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) - - # probs = [BVProblem(BVPFunction{true}(f1!, bc1!), u0, tspan; nlls = Val(false)), - # BVProblem(BVPFunction{false}(f1, bc1), u0, tspan; nlls = Val(false)), - # BVProblem( - # BVPFunction{true}( - # f1!, (bc1_a!, bc1_b!); bcresid_prototype, twopoint = Val(true)), - # u0, - # tspan; - # nlls = Val(false)), - # BVProblem( - # BVPFunction{false}(f1, (bc1_a, bc1_b); bcresid_prototype, twopoint = Val(true)), - # u0, tspan; nlls = Val(false))] - - # algs = [] - - # if @load_preference("PrecompileShooting", true) - # push!(algs, - # Shooting(Tsit5(); nlsolve = NewtonRaphson(), - # jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))) - # end - - # if @load_preference("PrecompileMultipleShooting", true) - # push!(algs, - # MultipleShooting(10, - # Tsit5(); - # nlsolve = NewtonRaphson(), - # jac_alg = BVPJacobianAlgorithm(; - # bc_diffmode = AutoForwardDiff(; chunksize = 2), - # nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2)))) - # end - - # @compile_workload begin - # for prob in probs, alg in algs - # solve(prob, alg) - # end - # end - - # function f1_nlls!(du, u, p, t) - # du[1] = u[2] - # du[2] = -u[1] - # end - - # f1_nlls(u, p, t) = [u[2], -u[1]] - - # function bc1_nlls!(resid, sol, p, t) - # solₜ₁ = sol(0.0) - # solₜ₂ = sol(100.0) - # resid[1] = solₜ₁[1] - # resid[2] = solₜ₂[1] - 1 - # resid[3] = solₜ₂[2] + 1.729109 - # return nothing - # end - # bc1_nlls(sol, p, t) = [sol(0.0)[1], sol(100.0)[1] - 1, sol(1.0)[2] + 1.729109] - - # bc1_nlls_a!(resid, ua, p) = (resid[1] = ua[1]) - # bc1_nlls_b!(resid, ub, p) = (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109) - - # bc1_nlls_a(ua, p) = [ua[1]] - # bc1_nlls_b(ub, p) = [ub[1] - 1, ub[2] + 1.729109] - - # tspan = (0.0, 100.0) - # u0 = [0.0, 1.0] - # bcresid_prototype1 = Array{Float64}(undef, 3) - # bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2)) - - # probs = [ - # BVProblem( - # BVPFunction{true}(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1), - # u0, tspan; nlls = Val(true)), - # BVProblem( - # BVPFunction{false}(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1), - # u0, tspan; nlls = Val(true)), - # BVProblem( - # BVPFunction{true}(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!); - # bcresid_prototype = bcresid_prototype2, twopoint = Val(true)), - # u0, - # tspan; - # nlls = Val(true)), - # BVProblem( - # BVPFunction{false}(f1_nlls, (bc1_nlls_a, bc1_nlls_b); - # bcresid_prototype = bcresid_prototype2, twopoint = Val(true)), - # u0, - # tspan; - # nlls = Val(true))] - - # algs = [] - - # if @load_preference("PrecompileShootingNLLS", true) - # append!(algs, - # [ - # Shooting(Tsit5(); nlsolve = TrustRegion(), - # jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))), - # Shooting(Tsit5(); nlsolve = GaussNewton(), - # jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))]) - # end - - # if @load_preference("PrecompileMultipleShootingNLLS", true) - # append!(algs, - # [ - # MultipleShooting(10, - # Tsit5(); - # nlsolve = TrustRegion(), - # jac_alg = BVPJacobianAlgorithm(; - # bc_diffmode = AutoForwardDiff(; chunksize = 2), - # nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2))), - # MultipleShooting(10, - # Tsit5(); - # nlsolve = GaussNewton(), - # jac_alg = BVPJacobianAlgorithm(; - # bc_diffmode = AutoForwardDiff(; chunksize = 2), - # nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2)))]) - # end - - # @compile_workload begin - # for prob in probs, alg in algs - # solve(prob, alg) - # end - # end + f1! = @closure (du, u, p, t) -> begin + du[1] = u[2] + du[2] = 0 + end + f1 = @closure (u, p, t) -> [u[2], 0] + + bc1! = @closure (residual, u, p, t) -> begin + residual[1] = u[1][1] - 5 + residual[2] = u[lastindex(u)][1] + end + + bc1 = @closure (u, p, t) -> [u[1][1] - 5, u[lastindex(u)][1]] + + bc1_a! = @closure (residual, ua, p) -> (residual[1] = ua[1] - 5) + bc1_b! = @closure (residual, ub, p) -> (residual[1] = ub[1]) + + bc1_a = @closure (ua, p) -> [ua[1] - 5] + bc1_b = @closure (ub, p) -> [ub[1]] + + tspan = (0.0, 5.0) + 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))] + + algs = [] + + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) + + if Preferences.@load_preference("PrecompileMIRK", true) + append!(algs, [MIRK2(; jac_alg), MIRK4(; jac_alg), MIRK6(; jac_alg)]) + end + + @compile_workload begin + for prob in probs, alg in algs + solve(prob, alg; dt = 0.2) + end + end + + f1_nlls! = @closure (du, u, p, t) -> begin + du[1] = u[2] + du[2] = -u[1] + end + + f1_nlls = @closure (u, p, t) -> [u[2], -u[1]] + + bc1_nlls! = @closure (resid, sol, p, t) -> begin + solₜ₁ = sol[1] + solₜ₂ = sol[lastindex(sol)] + resid[1] = solₜ₁[1] + resid[2] = solₜ₂[1] - 1 + resid[3] = solₜ₂[2] + 1.729109 + return nothing + end + bc1_nlls = @closure (sol, p, t) -> [ + sol[1][1], sol[lastindex(sol)][1] - 1, sol[lastindex(sol)][2] + 1.729109] + + bc1_nlls_a! = @closure (resid, ua, p) -> (resid[1] = ua[1]) + bc1_nlls_b! = @closure (resid, ub, p) -> (resid[1] = ub[1] - 1; + resid[2] = ub[2] + 1.729109) + + bc1_nlls_a = @closure (ua, p) -> [ua[1]] + bc1_nlls_b = @closure (ub, p) -> [ub[1] - 1, ub[2] + 1.729109] + + tspan = (0.0, 100.0) + u0 = [0.0, 1.0] + bcresid_prototype1 = Array{Float64}(undef, 3) + bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2)) + + probs = [ + BVProblem(BVPFunction(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1), + 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))] + + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) + + nlsolvers = [LevenbergMarquardt(; disable_geodesic = Val(true)), GaussNewton()] + + algs = [] + + if Preferences.@load_preference("PrecompileMIRKNLLS", false) + for nlsolve in nlsolvers + append!(algs, [MIRK2(; jac_alg, nlsolve), MIRK6(; jac_alg, nlsolve)]) + end + end + + @compile_workload begin + for prob in probs, alg in algs + solve(prob, alg; dt = 0.2, abstol = 1e-2) + end + end + + bc1! = @closure (residual, u, p, t) -> begin + residual[1] = u(0.0)[1] - 5 + residual[2] = u(5.0)[1] + end + bc1 = @closure (u, p, t) -> [u(0.0)[1] - 5, u(5.0)[1]] + + tspan = (0.0, 5.0) + u0 = [5.0, -3.5] + bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) + + probs = [BVProblem(BVPFunction{true}(f1!, bc1!), u0, tspan; nlls = Val(false)), + BVProblem(BVPFunction{false}(f1, bc1), u0, tspan; nlls = Val(false)), + BVProblem( + BVPFunction{true}( + f1!, (bc1_a!, bc1_b!); bcresid_prototype, twopoint = Val(true)), + u0, + tspan; + nlls = Val(false)), + BVProblem( + BVPFunction{false}(f1, (bc1_a, bc1_b); bcresid_prototype, twopoint = Val(true)), + u0, tspan; nlls = Val(false))] + + algs = [] + + if @load_preference("PrecompileShooting", true) + push!(algs, + Shooting(Tsit5(); nlsolve = NewtonRaphson(), + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))) + end + + if @load_preference("PrecompileMultipleShooting", true) + push!(algs, + MultipleShooting(10, + Tsit5(); + nlsolve = NewtonRaphson(), + jac_alg = BVPJacobianAlgorithm(; + bc_diffmode = AutoForwardDiff(; chunksize = 2), + nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2)))) + end + + @compile_workload begin + for prob in probs, alg in algs + solve(prob, alg) + end + end + + bc1_nlls! = @closure (resid, sol, p, t) -> begin + solₜ₁ = sol(0.0) + solₜ₂ = sol(100.0) + resid[1] = solₜ₁[1] + resid[2] = solₜ₂[1] - 1 + resid[3] = solₜ₂[2] + 1.729109 + return nothing + end + bc1_nlls = @closure (sol, p, t) -> [ + sol(0.0)[1], sol(100.0)[1] - 1, sol(1.0)[2] + 1.729109] + + tspan = (0.0, 100.0) + u0 = [0.0, 1.0] + bcresid_prototype1 = Array{Float64}(undef, 3) + bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2)) + + probs = [ + BVProblem( + BVPFunction{true}(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1), + u0, tspan; nlls = Val(true)), + BVProblem( + BVPFunction{false}(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1), + u0, tspan; nlls = Val(true)), + BVProblem( + BVPFunction{true}(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!); + bcresid_prototype = bcresid_prototype2, twopoint = Val(true)), + u0, + tspan; + nlls = Val(true)), + BVProblem( + BVPFunction{false}(f1_nlls, (bc1_nlls_a, bc1_nlls_b); + bcresid_prototype = bcresid_prototype2, twopoint = Val(true)), + u0, + tspan; + nlls = Val(true))] + + algs = [] + + if @load_preference("PrecompileShootingNLLS", true) + append!(algs, + [ + Shooting(Tsit5(); nlsolve = LevenbergMarquardt(; disable_geodesic = Val(true)), + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))), + Shooting(Tsit5(); nlsolve = GaussNewton(), + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)))]) + end + + if @load_preference("PrecompileMultipleShootingNLLS", true) + append!(algs, + [ + MultipleShooting(10, + Tsit5(); + nlsolve = LevenbergMarquardt(; disable_geodesic = Val(true)), + jac_alg = BVPJacobianAlgorithm(; + bc_diffmode = AutoForwardDiff(; chunksize = 2), + nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2))), + MultipleShooting(10, + Tsit5(); + nlsolve = GaussNewton(), + jac_alg = BVPJacobianAlgorithm(; + bc_diffmode = AutoForwardDiff(; chunksize = 2), + nonbc_diffmode = AutoSparseForwardDiff(; chunksize = 2)))]) + end + + @compile_workload begin + for prob in probs, alg in algs + solve(prob, alg; nlsolve_kwargs = (; abstol = 1e-2)) + end + end end export Shooting, MultipleShooting diff --git a/src/solve/mirk.jl b/src/solve/mirk.jl index 7736f207b..8d82633e9 100644 --- a/src/solve/mirk.jl +++ b/src/solve/mirk.jl @@ -160,10 +160,12 @@ function SciMLBase.solve!(cache::MIRKCache) return __build_solution(cache.prob, odesol, sol_nlprob) end -function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive; kwargs...) +function __perform_mirk_iteration( + cache::MIRKCache, abstol, adaptive; nlsolve_kwargs = (;), kwargs...) nlprob = __construct_nlproblem(cache, recursive_flatten(cache.y₀)) nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve) - sol_nlprob = __solve(nlprob, nlsolve_alg; abstol, kwargs..., alias_u0 = true) + sol_nlprob = __solve( + nlprob, nlsolve_alg; abstol, kwargs..., nlsolve_kwargs..., alias_u0 = true) recursive_unflatten!(cache.y₀, sol_nlprob.u) defect_norm = 2 * abstol diff --git a/test/misc/type_stability.jl b/test/misc/type_stability.jl deleted file mode 100644 index ef20e6325..000000000 --- a/test/misc/type_stability.jl +++ /dev/null @@ -1,62 +0,0 @@ -using BoundaryValueDiffEq, OrdinaryDiffEq, LinearAlgebra, Test - -f(u, p, t) = [p[1] * u[1] - p[2] * u[1] * u[2], p[3] * u[1] * u[2] - p[4] * u[2]] -function f!(du, u, p, t) - du[1] = p[1] * u[1] - p[2] * u[1] * u[2] - du[2] = p[3] * u[1] * u[2] - p[4] * u[2] -end - -bc(sol, p, t) = [sol[1][1] - 1, sol[end][2] - 2] -function bc!(res, sol, p, t) - res[1] = sol[1][1] - 1 - res[2] = sol[end][2] - 2 -end -twobc_a(ua, p) = [ua[1] - 1] -twobc_b(ub, p) = [ub[2] - 2] -twobc_a!(resa, ua, p) = (resa[1] = ua[1] - 1) -twobc_b!(resb, ub, p) = (resb[1] = ub[2] - 2) - -u0 = Float64[0, 0] -tspan = (0.0, 1.0) -p = [1.0, 1.0, 1.0, 1.0] -bcresid_prototype = (zeros(1), zeros(1)) - -# Multi-Point BVP -@testset "Multi-Point BVP" begin - mpbvp_iip = BVProblem(f!, bc!, u0, tspan, p) - mpbvp_oop = BVProblem(f, bc, u0, tspan, p) - - @testset "Shooting Methods" begin - @inferred solve(mpbvp_iip, Shooting(Tsit5())) - @inferred solve(mpbvp_oop, Shooting(Tsit5())) - @inferred solve(mpbvp_iip, MultipleShooting(5, Tsit5())) - @inferred solve(mpbvp_oop, MultipleShooting(5, Tsit5())) - end - - @testset "MIRK Methods" begin - for solver in (MIRK2(), MIRK3(), MIRK4(), MIRK5(), MIRK6()) - @inferred solve(mpbvp_iip, solver; dt = 0.2) - @inferred solve(mpbvp_oop, solver; dt = 0.2) - end - end -end - -# Two-Point BVP -@testset "Two-Point BVP" begin - tpbvp_iip = TwoPointBVProblem(f!, (twobc_a!, twobc_b!), u0, tspan, p; bcresid_prototype) - tpbvp_oop = TwoPointBVProblem(f, (twobc_a, twobc_b), u0, tspan, p) - - @testset "Shooting Methods" begin - @inferred solve(tpbvp_iip, Shooting(Tsit5())) - @inferred solve(tpbvp_oop, Shooting(Tsit5())) - @inferred solve(tpbvp_iip, MultipleShooting(5, Tsit5())) - @inferred solve(tpbvp_oop, MultipleShooting(5, Tsit5())) - end - - @testset "MIRK Methods" begin - for solver in (MIRK2(), MIRK3(), MIRK4(), MIRK5(), MIRK6()) - @inferred solve(tpbvp_iip, solver; dt = 0.2) - @inferred solve(tpbvp_oop, solver; dt = 0.2) - end - end -end diff --git a/test/misc/type_stability_tests.jl b/test/misc/type_stability_tests.jl new file mode 100644 index 000000000..f05b10637 --- /dev/null +++ b/test/misc/type_stability_tests.jl @@ -0,0 +1,70 @@ +@testitem "Type Stability" begin + using LinearAlgebra + + f(u, p, t) = [p[1] * u[1] - p[2] * u[1] * u[2], p[3] * u[1] * u[2] - p[4] * u[2]] + function f!(du, u, p, t) + du[1] = p[1] * u[1] - p[2] * u[1] * u[2] + du[2] = p[3] * u[1] * u[2] - p[4] * u[2] + end + + bc(sol, p, t) = [sol[1][1] - 1, sol[end][2] - 2] + function bc!(res, sol, p, t) + res[1] = sol[1][1] - 1 + res[2] = sol[end][2] - 2 + end + twobc_a(ua, p) = [ua[1] - 1] + twobc_b(ub, p) = [ub[2] - 2] + twobc_a!(resa, ua, p) = (resa[1] = ua[1] - 1) + twobc_b!(resb, ub, p) = (resb[1] = ub[2] - 2) + + u0 = Float64[0, 0] + tspan = (0.0, 1.0) + p = [1.0, 1.0, 1.0, 1.0] + bcresid_prototype = (zeros(1), zeros(1)) + + jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2)) + + # Multi-Point BVP + @testset "Multi-Point BVP" begin + mpbvp_iip = BVProblem(f!, bc!, u0, tspan, p; nlls = Val(false)) + mpbvp_oop = BVProblem(f, bc, u0, tspan, p; nlls = Val(false)) + + @testset "Shooting Methods" begin + @inferred solve(mpbvp_iip, Shooting(Tsit5(); jac_alg)) + @inferred solve(mpbvp_oop, Shooting(Tsit5(); jac_alg)) + @inferred solve(mpbvp_iip, MultipleShooting(5, Tsit5(); jac_alg)) + @inferred solve(mpbvp_oop, MultipleShooting(5, Tsit5(); jac_alg)) + end + + @testset "MIRK Methods" begin + for solver in (MIRK2(; jac_alg), MIRK3(; jac_alg), MIRK4(; jac_alg), + MIRK5(; jac_alg), MIRK6(; jac_alg)) + @inferred solve(mpbvp_iip, solver; dt = 0.2) + @inferred solve(mpbvp_oop, solver; dt = 0.2) + end + end + end + + # Two-Point BVP + @testset "Two-Point BVP" begin + tpbvp_iip = TwoPointBVProblem( + f!, (twobc_a!, twobc_b!), u0, tspan, p; bcresid_prototype, nlls = Val(false)) + tpbvp_oop = TwoPointBVProblem( + f, (twobc_a, twobc_b), u0, tspan, p; nlls = Val(false)) + + @testset "Shooting Methods" begin + @inferred solve(tpbvp_iip, Shooting(Tsit5(); jac_alg)) + @inferred solve(tpbvp_oop, Shooting(Tsit5(); jac_alg)) + @inferred solve(tpbvp_iip, MultipleShooting(5, Tsit5(); jac_alg)) + @inferred solve(tpbvp_oop, MultipleShooting(5, Tsit5(); jac_alg)) + end + + @testset "MIRK Methods" begin + for solver in (MIRK2(; jac_alg), MIRK3(; jac_alg), MIRK4(; jac_alg), + MIRK5(; jac_alg), MIRK6(; jac_alg)) + @inferred solve(tpbvp_iip, solver; dt = 0.2) + @inferred solve(tpbvp_oop, solver; dt = 0.2) + end + end + end +end