I’ve found an odd issue while using a previous solution as the initial guess in solving a boundary value problem using the SciML MIRK solvers and TwoPointBVProblem (part of the BoundaryValueDiffEq.jl package. According to the documentation I should be able to simply pass in the solution object as the initial guess, but doing so causes the behavior I describe below.
The first solution uses an initial guess function, initialGuess(p, t). For the subsequent solutions my first technique evaluates the previously calculated solution as a vector of vectors, evaluated at uniform grid points in t (actually “xt” for my problem) as “sol(xt).u” and uses that for the initial guess in the next iteration. Here is a section of my code:
function initialGuess(p, t)
x = t * p[3]
θe = (234.26988 * log(-4.16207e4 * x + 1.25761e3))/p[1]
n̂ = ((-5.9645e13 * x)/(1 - 4.5192e2 * x^2) + 2.8991e12)*sqrt((xk * p[1])/(2 * π * me))/(p[6]/ec)
ψ = (-6.4426e-2 * log(-33.8580 * x + 1.0))/(xke * p[1])
Γe = 0.05/p[6]
ûe = (3.1446 * x^2 - 0.1693 * x + 0.02556)/(xk2 * p[1] * (p[6]/ec))
return [θe, n̂, ψ, Γe, ûe]
end
function smstpc(te, tr, tc, jtot, din, ϕ₀, xphcl)
d, phie, js, pcs, phic, jc, jie = setup(te, tc, tr, din, ϕ₀, xphcl)
p = [te, tc/te, d, pcs, jtot/js, js, jie/js, jc/js]
Δx = 0.01
xt = (0.0 : Δx : 1.0)
prob = TwoPointBVProblem(plasma_transport_eqns!, (emitter_bcs!, collector_bcs!), initialGuess, (0.0, 1.0), p; bcresid_prototype = (zeros(3), zeros(2)))
@time sol1 = solve(prob, MIRK4(nlsolve = RobustMultiNewton()), dt = Δx; adaptive = true, maxiters = 1000)
v̂_e = emitter_sheath(p[5], p[7], sol1.u[1][2], sol1.u[1][1])
v̂_c = collector_sheath(p[5], p[8], sol1.u[end][2], sol1.u[end][1], p[2])
volts = phie - phic + xke*p[1]*(v̂_c - v̂_e + sol1.u[end][3])
Δj = 0.002
jparam = jtot - Δj
j_store = [jtot]
v_store = [volts]
while jparam > 0.02
p[5] = jparam/js
prob = TwoPointBVProblem(plasma_transport_eqns!, (emitter_bcs!, collector_bcs!), sol1(xt).u, (0.0, 1.0), p; bcresid_prototype = (zeros(3), zeros(2)))
@time sol = solve(prob, MIRK4(nlsolve = RobustMultiNewton()), dt = Δx; adaptive = true, maxiters = 1000)
v̂_e = emitter_sheath(p[5], p[7], sol.u[1][2], sol.u[1][1])
v̂_c = collector_sheath(p[5], p[8], sol.u[end][2], sol.u[end][1], p[2])
volts = phie - phic + xke*p[1]*(v̂_c - v̂_e + sol.u[end][3])
if abs(volts) < 3.0 && SciMLBase.successful_retcode(sol.retcode)
push!(j_store, jparam)
push!(v_store, volts)
end
jparam -= Δj
sol1 = sol
end
return (v_store, j_store)
end
and here is the timing output from @time:
72.004986 seconds (7.77 M allocations: 382.433 MiB, 0.19% gc time, 99.91% compilation time: 99% of which was recompilation)
74.480700 seconds (6.31 M allocations: 311.488 MiB, 0.08% gc time, 99.94% compilation time: 100% of which was recompilation)
0.222142 seconds (3.68 M allocations: 162.491 MiB, 12.21% gc time)
0.048849 seconds (801.63 k allocations: 39.341 MiB, 10.25% gc time)
0.031551 seconds (458.75 k allocations: 23.289 MiB, 18.92% gc time)
0.047677 seconds (778.69 k allocations: 38.460 MiB, 10.09% gc time)
0.030927 seconds (457.73 k allocations: 23.236 MiB, 19.58% gc time)
0.031191 seconds (455.68 k allocations: 23.135 MiB, 15.47% gc time)
0.025282 seconds (452.62 k allocations: 22.981 MiB)
0.033302 seconds (506.38 k allocations: 25.389 MiB, 16.48% gc time)
0.037231 seconds (566.10 k allocations: 28.082 MiB, 13.05% gc time)
0.037937 seconds (600.36 k allocations: 29.484 MiB, 15.06% gc time)
0.041759 seconds (671.51 k allocations: 32.716 MiB, 14.02% gc time)
0.050617 seconds (810.85 k allocations: 38.744 MiB, 10.71% gc time)
0.061268 seconds (1.02 M allocations: 47.769 MiB, 9.42% gc time)
So after the recompilation is over the solutions are obtained quickly. However, according to the documentation this is not the recommended method. The solution should be passed in directly, like this:
function smstpc(te, tr, tc, jtot, din, ϕ₀, xphcl)
d, phie, js, pcs, phic, jc, jie = setup(te, tc, tr, din, ϕ₀, xphcl)
p = [te, tc/te, d, pcs, jtot/js, js, jie/js, jc/js]
Δx = 0.01
# xt = (0.0 : Δx : 1.0)
prob = TwoPointBVProblem(plasma_transport_eqns!, (emitter_bcs!, collector_bcs!), initialGuess, (0.0, 1.0), p; bcresid_prototype = (zeros(3), zeros(2)))
@time sol = solve(prob, MIRK4(nlsolve = RobustMultiNewton()), dt = Δx; adaptive = true, maxiters = 1000)
v̂_e = emitter_sheath(p[5], p[7], sol.u[1][2], sol.u[1][1])
v̂_c = collector_sheath(p[5], p[8], sol.u[end][2], sol.u[end][1], p[2])
volts = phie - phic + xke*p[1]*(v̂_c - v̂_e + sol.u[end][3])
Δj = 0.002
jparam = jtot - Δj
j_store = [jtot]
v_store = [volts]
while jparam > 0.02
p[5] = jparam/js
prob = TwoPointBVProblem(plasma_transport_eqns!, (emitter_bcs!, collector_bcs!), sol, (0.0, 1.0), p; bcresid_prototype = (zeros(3), zeros(2)))
@time sol = solve(prob, MIRK4(nlsolve = RobustMultiNewton()), dt = Δx; adaptive = true, maxiters = 1000)
v̂_e = emitter_sheath(p[5], p[7], sol.u[1][2], sol.u[1][1])
v̂_c = collector_sheath(p[5], p[8], sol.u[end][2], sol.u[end][1], p[2])
volts = phie - phic + xke*p[1]*(v̂_c - v̂_e + sol.u[end][3])
if abs(volts) < 3.0 && SciMLBase.successful_retcode(sol.retcode)
push!(j_store, jparam)
push!(v_store, volts)
end
jparam -= Δj
end
return (v_store, j_store)
end
And here is the @time output for that:
89.509736 seconds (56.56 M allocations: 2.838 GiB, 2.08% gc time, 99.95% compilation time)
154.698938 seconds (17.89 M allocations: 949.767 MiB, 0.55% gc time, 99.99% compilation time)
233.326477 seconds (17.89 M allocations: 949.830 MiB, 0.15% gc time, 99.99% compilation time)
4419.765366 seconds (17.89 M allocations: 949.537 MiB, 0.01% gc time, 100.00% compilation time)
The execution time increases dramatically with each iteration, and it is nearly all compilation time (not re-compilation time anymore). Is this the expected behavior, or have I got something wrong?