SciML collocation BVP solvers overwrite the initial guess - bug or feature?

Specific Domains Numerics
1

I’ve encountered a behavior that I didn’t expect when solving BVPs with SciML collocation methods and using the continuation method. When you call “solve” and are using a previous solution object as the initial guess, that initial guess is overwritten by the new solution. Is that intended? Is it documented somewhere?

Here’s a code example:

using BoundaryValueDiffEq, NonlinearSolve
tspan = (0.0, 1.0)
function f!(du, u, p, t)
    cond = p[1]
    vol_heat = p[2]
    du[1] = -u[2]/cond
    du[2] = vol_heat
    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)
    tref = 20.0
    res_b[1] = u_b[3] * (u_b[1] - tref) - u_b[2]
 end
u_guess = [ [100.0, 0.0, 0.006666666666666668],
    [99.5, 0.020000000000000004, 0.006666666666666668],
    [98.0, 0.04000000000000001, 0.006666666666666668],
    [95.5, 0.060000000000000005, 0.006666666666666668],
    [92.0, 0.08000000000000002, 0.006666666666666668],
    [87.5, 0.1, 0.006666666666666668],
    [82.0, 0.12000000000000001, 0.006666666666666668],
    [75.5, 0.14, 0.006666666666666668],
    [68.0, 0.16000000000000003, 0.006666666666666668],
    [59.49999999999999, 0.18000000000000002, 0.006666666666666668],
    [50.0, 0.2, 0.006666666666666668] ]
p1 = [0.002, 0.2]
bvp1 = TwoPointBVProblem(f!, (bca!, bcb!), u_guess, tspan, p1; bcresid_prototype = (zeros(2), zeros(1)))
sol1 = solve(bvp1, MIRK5(nlsolve = NewtonRaphson()), dt = 0.1, adaptive = false; nlsolve_kwargs = (; maxiters = 10, show_trace = Val(false)))

println("sol1.u[1][1]=     ",sol1.u[1][1], "  sol1.u[1][2]=     ", sol1.u[1][2], "  sol1.u[1][3]=     ", sol1.u[1][3])

p1[1] = 0.0015
p1[2] = 0.03
sol_guess = deepcopy(sol1)
println("sol_guess.u[1][1]=",sol_guess.u[1][1], "  sol_guess.u[1][2]=", sol_guess.u[1][2], "  sol_guess.u[1][3]=", sol_guess.u[1][3])
println("Solving")
bvp2 = TwoPointBVProblem(f!, (bca!, bcb!), sol_guess, tspan, p1; bcresid_prototype = (zeros(2), zeros(1)))
sol2 = solve(bvp2, MIRK5(nlsolve = NewtonRaphson()), dt = 0.1, adaptive = false; nlsolve_kwargs = (; maxiters = 10, show_trace = Val(false)));
println("sol1.u[1][1]=     ",sol1.u[1][1], "  sol1.u[1][2]=     ", sol1.u[1][2], "  sol1.u[1][3]=     ", sol1.u[1][3])
println("sol_guess.u[1][1]=",sol_guess.u[1][1], "  sol_guess.u[1][2]=", sol_guess.u[1][2], "  sol_guess.u[1][3]=", sol_guess.u[1][3])
println("sol2.u[1][1]=     ",sol2.u[1][1], "  sol2.u[1][2]=     ", sol2.u[1][2], "  sol2.u[1][3]=     ", sol2.u[1][3])

And here is the output:

sol1.u[1][1]=     100.0  sol1.u[1][2]=     0.0  sol1.u[1][3]=     0.006666666666666668
sol_guess.u[1][1]=100.0  sol_guess.u[1][2]=0.0  sol_guess.u[1][3]=0.006666666666666668
Solving
sol1.u[1][1]=     100.0  sol1.u[1][2]=     0.0  sol1.u[1][3]=     0.006666666666666668
sol_guess.u[1][1]=100.0  sol_guess.u[1][2]=3.604278872980289e-17  sol_guess.u[1][3]=0.00042857142857142747
sol2.u[1][1]=     100.0  sol2.u[1][2]=     3.604278872980289e-17  sol2.u[1][3]=     0.00042857142857142747

The values in “sol_guess” have become identical to the values in “sol2” after the call “sol2 = solve(bvp2,…)”.

In my continuation scheme it may be that I try to solve using one method and it fails to converge. I then want to try to solve using a different method (one that might be more robust) and the same initial guess as was used for the failed solve. Unless I’ve saved my previous initial guess as a separate data item (such as “sol1” in the example I’ve given) that initial guess has been destroyed. This issue can be easily fixed, but if you don’t expect the initial guess to be overwritten it can cause a lot of debugging in order to find the issue.

2

This likely isn’t intended. Sometime in the last ~year or so we’ve added an actual aliasing system to most of the sciml solvers, but it is very possible that the BVP solvers haven’t gotten it yet.

Our default generally is to not alias user provided arrays unless the user specifically tells us to.

3

Should I file a bug report on GitHub or not?

4

Yes please.