Merge pull request #286 from ErikQQY/qqy/fix_dt

Author: ErikQQY

lib/BoundaryValueDiffEqCore/src/utils.jl

@@ -207,7 +207,6 @@ end
 function __extract_problem_details(
         prob, u0::SciMLBase.ODESolution; dt = 0.0, check_positive_dt::Bool = false)
     # Problem passes in a initial guess function
-    check_positive_dt && dt ≤ 0 && throw(ArgumentError("dt must be positive"))
     _u0 = first(u0.u)
     _t = u0.t
     return Val(true), eltype(_u0), length(_u0), (length(_t) - 1), _u0

lib/BoundaryValueDiffEqFIRK/src/firk.jl

@@ -162,14 +162,7 @@ function init_nested(prob::BVProblem, alg::AbstractFIRK; dt = 0.0,
 
     prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = X) : prob
 
-    if isa(prob.u0, AbstractArray) && eltype(prob.u0) <: AbstractVector
-        u0_mat = hcat(prob.u0...)
-        avg_u0 = vec(sum(u0_mat, dims = 2)) / size(u0_mat, 2)
-    else
-        avg_u0 = prob.u0
-    end
-
-    K0 = repeat(avg_u0, 1, stage) # Somewhat arbitrary initialization of K
+    K0 = __K0_on_u0(prob.u0, stage) # Somewhat arbitrary initialization of K
 
     nestprob_p = zeros(T, M + 2)
     nest_tol = alg.nest_tol

lib/BoundaryValueDiffEqFIRK/src/utils.jl

@@ -39,3 +39,20 @@ function BoundaryValueDiffEqCore.__append_similar!(
     append!(x, VectorOfArray([similar(last(x)) for _ in 1:N]))
     return x
 end
+
+@inline __K0_on_u0(u0::AbstractArray, stage) = repeat(u0, 1, stage)
+@inline function __K0_on_u0(u0::AbstractVector{<:AbstractArray}, stage)
+    u0_mat = hcat(u0...)
+    avg_u0 = vec(sum(u0_mat, dims = 2)) / size(u0_mat, 2)
+    return repeat(avg_u0, 1, stage)
+end
+@inline function __K0_on_u0(u0::AbstractVectorOfArray, stage)
+    u0_mat = hcat(u0.u...)
+    avg_u0 = vec(sum(u0_mat, dims = 2)) / size(u0_mat, 2)
+    return repeat(avg_u0, 1, stage)
+end
+@inline function __K0_on_u0(u0::SciMLBase.ODESolution, stage)
+    u0_mat = hcat(u0.u...)
+    avg_u0 = vec(sum(u0_mat, dims = 2)) / size(u0_mat, 2)
+    return repeat(avg_u0, 1, stage)
+end

lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl

@@ -437,11 +437,9 @@ end =#
 
     bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2),
         pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
-    @test_broken SciMLBase.successful_retcode(solve(
-        bvp4, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
+    SciMLBase.successful_retcode(solve(bvp4, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
 
     bvp5 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), DiffEqArray(sol3.u, sol3.t),
         (0, pi / 2), pi / 2; bcresid_prototype = (zeros(1), zeros(1)))
-    @test_broken SciMLBase.successful_retcode(solve(
-        bvp5, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
+    SciMLBase.successful_retcode(solve(bvp5, RadauIIa5(; nested_nlsolve = true), dt = 0.05))
 end

test/misc/initial_guess_tests.jl

@@ -0,0 +1,39 @@
+@testitem "Initial guess providing" begin
+    using BoundaryValueDiffEq, RecursiveArrayTools
+    tspan = (0.0, pi / 2)
+    function simplependulum!(du, u, p, t)
+        θ = u[1]
+        dθ = u[2]
+        du[1] = dθ
+        du[2] = -9.81 * sin(θ)
+    end
+    function bc!(residual, u, p, t)
+        residual[1] = u(pi / 4)[1] + pi / 2
+        residual[2] = u(pi / 2)[1] - pi / 2
+    end
+    u0 = [pi / 2, pi / 2]
+    prob = BVProblem(simplependulum!, bc!, u0, tspan)
+    sol1 = solve(prob, MIRK4(), dt = 0.05)
+
+    # Solution
+    prob1 = BVProblem(simplependulum!, bc!, sol1, tspan)
+    sol2 = solve(prob1, MIRK4())
+    @test SciMLBase.successful_retcode(sol2)
+
+    sol3 = solve(prob1, RadauIIa5())
+    @test SciMLBase.successful_retcode(sol3)
+
+    sol4 = solve(prob1, LobattoIIIa4(nested_nlsolve = true))
+    @test SciMLBase.successful_retcode(sol4)
+
+    # VectorOfArray
+    prob2 = BVProblem(simplependulum!, bc!, VectorOfArray(sol1.u), tspan)
+    sol2 = solve(prob2, MIRK4())
+    @test SciMLBase.successful_retcode(sol2)
+
+    sol3 = solve(prob2, RadauIIa5())
+    @test SciMLBase.successful_retcode(sol3)
+
+    sol4 = solve(prob2, LobattoIIIa4(nested_nlsolve = true))
+    @test SciMLBase.successful_retcode(sol4)
+end