fix: copy initials to u0 if u0 not provided to remake

Author: AayushSabharwal

src/systems/nonlinear/initializesystem.jl

@@ -543,6 +543,7 @@ function SciMLBase.remake_initialization_data(
     if u0 === missing && p === missing
         return odefn.initialization_data
     end
+
     if !(eltype(u0) <: Pair) && !(eltype(p) <: Pair)
         oldinitdata = odefn.initialization_data
         oldinitdata === nothing && return nothing
@@ -658,6 +659,18 @@ function SciMLBase.late_binding_update_u0_p(
         prob, sys::AbstractSystem, u0, p, t0, newu0, newp)
     supports_initialization(sys) || return newu0, newp
     u0 === missing && return newu0, (p === missing ? copy(newp) : newp)
+    # If the user passes `p` to `remake` but not `u0` and `u0` isn't empty,
+    # and if the system supports initialization (so it has initial parameters),
+    # and if the initialization solves for `u0`,
+    # THEN copy the values of `Initial`s to `newu0`.
+    if u0 === missing && newu0 !== nothing && p !== missing && supports_initialization(sys) && prob.f.initialization_data !== nothing && prob.f.initialization_data.initializeprobmap !== nothing
+        if ArrayInterface.ismutable(newu0)
+            copyto!(newu0, getu(sys, Initial.(unknowns(sys)))(newp))
+        else
+            T = StaticArrays.similar_type(newu0)
+            newu0 = T(getu(sys, Initial.(unknowns(sys)))(newp))
+        end
+    end
     # non-symbolic u0 updates initials...
     if !(eltype(u0) <: Pair)
         # if `p` is not provided or is symbolic

test/initializationsystem.jl

@@ -1496,3 +1496,24 @@ end
     @test integ3.u ≈ [2.0, 3.0]
     @test integ3.ps[c1] ≈ 2.0
 end
+
+@testset "Issue#3570: `Initial`s are copied to `u0` if `u0` not provided to `remake`" begin
+    @parameters g
+    @variables x(t) [state_priority = 10] y(t) λ(t)
+    eqs = [D(D(x)) ~ λ * x
+           D(D(y)) ~ λ * y - g
+           x^2 + y^2 ~ 1]
+    @mtkbuild pend = ODESystem(eqs, t)
+
+    prob = ODEProblem(
+        pend, [x => (√2 / 2)], (0.0, 1.5), [g => 1], guesses = [λ => 1, y => √2 / 2])
+    sol = solve(prob)
+
+    setter = setsym_oop(prob, [Initial(x)])
+    (u0, p) = setter(prob, [0.8])
+
+    new_prob = remake(prob; p, initializealg = BrownFullBasicInit())
+    @test new_prob[x] ≈ 0.8
+    new_sol = solve(new_prob)
+    @test new_sol[x, 1] ≈ 0.8
+end