Skip to content

Commit bf92edb

Browse files
fix: copy initials to u0 if u0 not provided to remake
1 parent b1ccc75 commit bf92edb

File tree

2 files changed

+33
-0
lines changed

2 files changed

+33
-0
lines changed

src/systems/nonlinear/initializesystem.jl

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -593,6 +593,18 @@ function SciMLBase.late_binding_update_u0_p(
593593
prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
594594
supports_initialization(sys) || return newu0, newp
595595
u0 === missing && return newu0, (p === missing ? copy(newp) : newp)
596+
# If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
597+
# and if the system supports initialization (so it has initial parameters),
598+
# and if the initialization solves for `u0`,
599+
# THEN copy the values of `Initial`s to `newu0`.
600+
if u0 === missing && newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
601+
if ArrayInterface.ismutable(newu0)
602+
copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
603+
else
604+
T = StaticArrays.similar_type(newu0)
605+
newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
606+
end
607+
end
596608
# non-symbolic u0 updates initials...
597609
if !(eltype(u0) <: Pair)
598610
# if `p` is not provided or is symbolic

test/initializationsystem.jl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1507,3 +1507,24 @@ end
15071507
@inferred remake(prob; u0 = 2 .* prob.u0, p = prob.p)
15081508
@inferred solve(prob)
15091509
end
1510+
1511+
@testset "Issue#3570: `Initial`s are copied to `u0` if `u0` not provided to `remake`" begin
1512+
@parameters g
1513+
@variables x(t) [state_priority = 10] y(t) λ(t)
1514+
eqs = [D(D(x)) ~ λ * x
1515+
D(D(y)) ~ λ * y - g
1516+
x^2 + y^2 ~ 1]
1517+
@mtkbuild pend = ODESystem(eqs, t)
1518+
1519+
prob = ODEProblem(
1520+
pend, [x => (2 / 2)], (0.0, 1.5), [g => 1], guesses ==> 1, y => 2 / 2])
1521+
sol = solve(prob)
1522+
1523+
setter = setsym_oop(prob, [Initial(x)])
1524+
(u0, p) = setter(prob, [0.8])
1525+
1526+
new_prob = remake(prob; p, initializealg = BrownFullBasicInit())
1527+
@test new_prob[x] 0.8
1528+
new_sol = solve(new_prob)
1529+
@test new_sol[x, 1] 0.8
1530+
end

0 commit comments

Comments
 (0)