Merge pull request #3028 from hersle/fix_ics_unknowns

Override defaults when setting initial conditions with mapped unknowns(sys)

src/systems/nonlinear/initializesystem.jl

@@ -44,24 +44,24 @@ function generate_initializesystem(sys::ODESystem;
                 y = get(schedule.dummy_sub, x[1], x[1])
                 y = ModelingToolkit.fixpoint_sub(y, full_diffmap)
 
-                if y isa Symbolics.Arr
+                if y ∈ set_full_states
+                    # defer initialization until defaults are merged below
+                    push!(filtered_u0, y => x[2])
+                elseif y isa Symbolics.Arr
+                    # scalarize array # TODO: don't scalarize arrays
                     _y = collect(y)
-
-                    # TODO: Don't scalarize arrays
-                    for i in 1:length(_y)
+                    for i in eachindex(_y)
                         push!(filtered_u0, _y[i] => x[2][i])
                     end
-                elseif y isa ModelingToolkit.BasicSymbolic
+                elseif y isa Symbolics.BasicSymbolic
                     # y is a derivative expression expanded
                     # add to the initialization equations
                     push!(eqs_ics, y ~ x[2])
-                elseif y ∈ set_full_states
-                    push!(filtered_u0, y => x[2])
                 else
                     error("Initialization expression $y is currently not supported. If its a higher order derivative expression, then only the dummy derivative expressions are supported.")
                 end
             end
-            filtered_u0 = filtered_u0 isa Pair ? todict([filtered_u0]) : todict(filtered_u0)
+            filtered_u0 = todict(filtered_u0)
         end
     else
         dd_guess = Dict()

test/initializationsystem.jl

@@ -486,3 +486,25 @@ sys = extend(sysx, sysy)
     ssys = structural_simplify(sys)
     @test_throws ArgumentError ODEProblem(ssys, [x => 3], (0, 1), []) # y should have a guess
 end
+
+# https://github.com/SciML/ModelingToolkit.jl/issues/3025
+@testset "Override defaults when setting initial conditions with unknowns(sys) or similar" begin
+    @variables x(t) y(t)
+
+    # system 1 should solve to x = 1
+    ics1 = [x => 1]
+    sys1 = ODESystem([D(x) ~ 0], t; defaults = ics1, name = :sys1) |> structural_simplify
+    prob1 = ODEProblem(sys1, [], (0.0, 1.0), [])
+    sol1 = solve(prob1, Tsit5())
+    @test all(sol1[x] .== 1)
+
+    # system 2 should solve to x = y = 2
+    sys2 = extend(
+        sys1,
+        ODESystem([D(y) ~ 0], t; initialization_eqs = [y ~ 2], name = :sys2)
+    ) |> structural_simplify
+    ics2 = unknowns(sys1) .=> 2 # should be equivalent to "ics2 = [x => 2]"
+    prob2 = ODEProblem(sys2, ics2, (0.0, 1.0), []; fully_determined = true)
+    sol2 = solve(prob2, Tsit5())
+    @test all(sol2[x] .== 2) && all(sol2[y] .== 2)
+end