diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl index bb8f0a41f..43dd5d334 100644 --- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl +++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl @@ -257,7 +257,16 @@ function init_nested( vecf, vecbc end - prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = u0) : prob + # Initial guess objects (`ODESolution`, `VectorOfArray`, functions, ...) must be + # stripped down to the extracted `u0` vector here: under RecursiveArrayTools v4 + # `AbstractVectorOfArray <: AbstractArray`, so a plain `isa AbstractArray` check + # would embed e.g. an entire previous solution's type in the cache and force + # recompilation of all downstream code against it (issue #500). + prob_ = if !(prob.u0 isa AbstractArray) || prob.u0 isa AbstractVectorOfArray + remake(prob; u0 = u0) + else + prob + end # Somewhat arbitrary initialization of K K0 = __K0_on_u0(prob, stage; tune_parameters = tune_parameters) @@ -424,7 +433,16 @@ function init_expanded( vecf, vecbc end - prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = u0) : prob + # Initial guess objects (`ODESolution`, `VectorOfArray`, functions, ...) must be + # stripped down to the extracted `u0` vector here: under RecursiveArrayTools v4 + # `AbstractVectorOfArray <: AbstractArray`, so a plain `isa AbstractArray` check + # would embed e.g. an entire previous solution's type in the cache and force + # recompilation of all downstream code against it (issue #500). + prob_ = if !(prob.u0 isa AbstractArray) || prob.u0 isa AbstractVectorOfArray + remake(prob; u0 = u0) + else + prob + end return FIRKCacheExpand{iip, T, typeof(diffcache), tune_parameters}( alg_order(alg), stage, M, size(u0), f, bc, prob_, prob.problem_type, prob.p, diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl index 4bcc47f22..7a71fb001 100644 --- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl @@ -566,3 +566,42 @@ end @test bvp2.u0 == u_guess end =# + +# https://github.com/SciML/BoundaryValueDiffEq.jl/issues/500 +# An initial-guess object (previous solution, VectorOfArray, ...) must be +# stripped to a plain vector before being stored in the cache's problem. +# Under RecursiveArrayTools v4 solutions are `AbstractArray`s, so without the +# explicit strip the entire solution type gets embedded in the cache type and +# every downstream method recompiles against it (~20x compile-time blowup). +@testset "Initial-guess object does not leak into cache type (issue 500)" begin + using SciMLBase, RecursiveArrayTools + + tspan = (0.0, 1.0) + function f!(du, u, p, t) + du[1] = -u[2] / p[1] + du[2] = p[2] + du[3] = 0.0 + end + function bca!(res_a, u_a, p) + res_a[1] = u_a[2] + res_a[2] = u_a[1] - 100.0 + end + function bcb!(res_b, u_b, p) + res_b[1] = u_b[3] * (u_b[1] - 20.0) - u_b[2] + end + u_guess = [[100.0 - 0.5 * i, 0.02 * i, 0.006666666666666668] for i in 0:10] + p = [0.002, 0.2] + bvp1 = TwoPointBVProblem( + f!, (bca!, bcb!), u_guess, tspan, p; bcresid_prototype = (zeros(2), zeros(1)) + ) + sol1 = solve(bvp1, RadauIIa5(), dt = 0.1, adaptive = false) + @test SciMLBase.successful_retcode(sol1) + + bvp2 = remake(bvp1, p = [0.0015, 0.03], u0 = sol1) + for nested in (true, false) + cache = init(bvp2, RadauIIa5(; nested_nlsolve = nested), dt = 0.1, adaptive = false) + @test cache.prob.u0 isa Vector{Float64} + end + sol2 = solve(bvp2, RadauIIa5(), dt = 0.1, adaptive = false) + @test SciMLBase.successful_retcode(sol2) +end diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl index 78c401ecd..70ed43105 100644 --- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl +++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl @@ -235,7 +235,16 @@ function SciMLBase.__init( vecf, vecbc end - prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = u0) : prob + # Initial guess objects (`ODESolution`, `VectorOfArray`, functions, ...) must be + # stripped down to the extracted `u0` vector here: under RecursiveArrayTools v4 + # `AbstractVectorOfArray <: AbstractArray`, so a plain `isa AbstractArray` check + # would embed e.g. an entire previous solution's type in the cache and force + # recompilation of all downstream code against it (issue #500). + prob_ = if !(prob.u0 isa AbstractArray) || prob.u0 isa AbstractVectorOfArray + remake(prob; u0 = u0) + else + prob + end return MIRKCache{iip, T, use_both, typeof(diffcache), tune_parameters}( alg_order(alg), stage, N, size(u0), f, bc, prob_, prob.problem_type, prob.p, alg, diff --git a/lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl index 5af87b611..1fc789588 100644 --- a/lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl +++ b/lib/BoundaryValueDiffEqMIRK/test/Core/mirk_basic_tests.jl @@ -382,7 +382,7 @@ end @testset "Compatibility with StaticArrays" begin using StaticArrays - const g = 9.81 + g = 9.81 L = 1.0 tspan = (0.0, pi / 2) function simplependulum!(du, u, p, t) @@ -750,3 +750,46 @@ end @test sol.u[1][1:2] ≈ a1 atol = 1.0e-8 @test sol.u[end][1:2] ≈ a2 atol = 1.0e-8 end + +# https://github.com/SciML/BoundaryValueDiffEq.jl/issues/500 +# An initial-guess object (previous solution, VectorOfArray, ...) must be +# stripped to a plain vector before being stored in the cache's problem. +# Under RecursiveArrayTools v4 solutions are `AbstractArray`s, so without the +# explicit strip the entire solution type gets embedded in the cache type and +# every downstream method recompiles against it (~20x compile-time blowup). +@testset "Initial-guess object does not leak into cache type (issue 500)" begin + using SciMLBase, RecursiveArrayTools + + tspan = (0.0, 1.0) + function f!(du, u, p, t) + du[1] = -u[2] / p[1] + du[2] = p[2] + du[3] = 0.0 + end + function bca!(res_a, u_a, p) + res_a[1] = u_a[2] + res_a[2] = u_a[1] - 100.0 + end + function bcb!(res_b, u_b, p) + res_b[1] = u_b[3] * (u_b[1] - 20.0) - u_b[2] + end + u_guess = [[100.0 - 0.5 * i, 0.02 * i, 0.006666666666666668] for i in 0:10] + p = [0.002, 0.2] + bvp1 = TwoPointBVProblem( + f!, (bca!, bcb!), u_guess, tspan, p; bcresid_prototype = (zeros(2), zeros(1)) + ) + sol1 = solve(bvp1, MIRK5(), dt = 0.1, adaptive = false) + @test SciMLBase.successful_retcode(sol1) + + bvp2 = remake(bvp1, p = [0.0015, 0.03], u0 = sol1) + cache = init(bvp2, MIRK5(), dt = 0.1, adaptive = false) + @test cache.prob.u0 isa Vector{Float64} + sol2 = solve(bvp2, MIRK5(), dt = 0.1, adaptive = false) + @test SciMLBase.successful_retcode(sol2) + + bvp3 = remake(bvp1, u0 = VectorOfArray(u_guess)) + cache = init(bvp3, MIRK5(), dt = 0.1, adaptive = false) + @test cache.prob.u0 isa Vector{Float64} + sol3 = solve(bvp3, MIRK5(), dt = 0.1, adaptive = false) + @test SciMLBase.successful_retcode(sol3) +end diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl index 33be886c4..bd4f684b6 100644 --- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl +++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl @@ -113,7 +113,16 @@ function SciMLBase.__init( vecf, vecbc end - prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = u0) : prob + # Initial guess objects (`ODESolution`, `VectorOfArray`, functions, ...) must be + # stripped down to the extracted `u0` vector here: under RecursiveArrayTools v4 + # `AbstractVectorOfArray <: AbstractArray`, so a plain `isa AbstractArray` check + # would embed e.g. an entire previous solution's type in the cache and force + # recompilation of all downstream code against it (issue #500). + prob_ = if !(prob.u0 isa AbstractArray) || prob.u0 isa AbstractVectorOfArray + remake(prob; u0 = u0) + else + prob + end return MIRKNCache{iip, T}( alg_order(alg), stage, M, size(u0), f, bc, prob_, prob.problem_type,