Merge pull request #2068 from SciML/myb/odeprob

Fix a cycle in dependency chasing

Author: YingboMa

src/structural_transformation/codegen.jl

@@ -114,12 +114,12 @@ function gen_nlsolve!(is_not_prepended_assignment, eqs, vars, u0map::AbstractDic
         tmp = [init_assignments]
         # `deps[init_assignments]` gives the dependency of `init_assignments`
         while true
-            next_assignments = reduce(vcat, deps[init_assignments])
+            next_assignments = unique(reduce(vcat, deps[init_assignments]))
             isempty(next_assignments) && break
             init_assignments = next_assignments
             push!(tmp, init_assignments)
         end
-        needed_assignments_idxs = reduce(vcat, unique(reverse(tmp)))
+        needed_assignments_idxs = unique(reduce(vcat, reverse(tmp)))
         needed_assignments = assignments[needed_assignments_idxs]
     end
 

src/structural_transformation/utils.jl

@@ -165,7 +165,7 @@ function sorted_incidence_matrix(ts::TransformationState, val = true; only_algeq
             push!(J, j)
         end
     end
-    sparse(I, J, val, nsrcs(g), ndsts(g))
+    sparse(I, J, val, nsrcs(graph), ndsts(graph))
 end
 
 ###

test/odaeproblem.jl

@@ -0,0 +1,59 @@
+using ModelingToolkit, ModelingToolkitStandardLibrary, Test
+using OrdinaryDiffEq
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Blocks
+
+function Segment(; name)
+    @named R = Resistor(; R = 1)
+    @named r = Resistor(; R = 1)
+    @named C = Capacitor(; C = 1)
+
+    @named p1 = Pin() # top-left
+    @named p2 = Pin() # top-right
+    @named n = Pin() # bottom
+
+    eqs = [connect(p1, R.p)
+           connect(R.n, p2, r.p)
+           connect(r.n, C.p)
+           connect(C.n, n)]
+
+    return ODESystem(eqs, t, [], [];
+                     name = name,
+                     systems = [r, R, C, n, p1, p2])
+end
+
+function Strip(; name)
+    num_segments = 10
+    # construct `num_segments` segments
+    segments = [Segment(; name = Symbol(:St_, seg))
+                for seg in 1:num_segments]
+
+    @named p1 = Pin() # top-left
+    @named p2 = Pin() # top-right
+    @named n = Pin() # bottom
+
+    eqs = [connect(p1, segments[1].p1)
+           connect(p2, segments[end].p2)
+           [connect(n, seg.n) for seg in segments]...
+           [connect(segments[i].p2, segments[i + 1].p1) for i in 1:(num_segments - 1)]...]
+
+    return ODESystem(eqs, t, [], []; name,
+                     systems = [p1, p2, n, segments...])
+end
+
+@variables t
+@named source = Voltage()
+@named c = Constant(k = 0.01)
+
+@named ground = Ground()
+@named strip = Strip()
+
+rc_eqs = [connect(c.output, source.V)
+          connect(source.p, strip.p1, strip.p2)
+          connect(strip.n, source.n, ground.g)]
+
+@named rc_model = ODESystem(rc_eqs, t, systems = [strip, c, source, ground])
+sys = structural_simplify(rc_model)
+
+prob = ODAEProblem(sys, [], (0, 10))
+@test_nowarn solve(prob, Tsit5())

test/runtests.jl

@@ -23,6 +23,7 @@ using SafeTestsets, Test
 @safetestset "JumpSystem Test" begin include("jumpsystem.jl") end
 @safetestset "Constraints Test" begin include("constraints.jl") end
 @safetestset "Reduction Test" begin include("reduction.jl") end
+@safetestset "ODAEProblem Test" begin include("odaeproblem.jl") end
 @safetestset "Components Test" begin include("components.jl") end
 @safetestset "print_tree" begin include("print_tree.jl") end
 @safetestset "Error Handling" begin include("error_handling.jl") end