Merge pull request #2803 from SciML/myb/tearing_cycle

Fix specialize overdetemined system in tearing

Author: ChrisRackauckas

downstream/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -89,45 +89,36 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
 
     ieqs = Int[]
     filtered_vars = BitSet()
-    seen_eqs = falses(nsrcs(graph))
+    free_eqs = free_equations(graph, var_sccs, var_eq_matching, varfilter)
+    is_overdetemined = !isempty(free_eqs)
     for vars in var_sccs
         for var in vars
             if varfilter(var)
                 push!(filtered_vars, var)
                 if var_eq_matching[var] !== unassigned
                     ieq = var_eq_matching[var]
-                    seen_eqs[ieq] = true
                     push!(ieqs, ieq)
                 end
             end
             var_eq_matching[var] = unassigned
         end
-        tear_block!(vargraph, vars,
-            var_eq_matching, ict, solvable_graph,
-            ieqs, filtered_vars, isder)
+        tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, ieqs,
+            filtered_vars, isder)
+        # If the systems is overdetemined, we cannot assume the free equations
+        # will not form algebraic loops with equations in the sccs.
+        if !is_overdetemined
+            vargraph.ne = 0
+            for var in vars
+                vargraph.matching[var] = unassigned
+            end
+        end
+        empty!(ieqs)
+        empty!(filtered_vars)
     end
-    free_eqs = findall(!, seen_eqs)
-    if !isempty(free_eqs)
+    if is_overdetemined
         free_vars = findall(x -> !(x isa Int), var_eq_matching)
-        tear_block!(vargraph, (),
-            var_eq_matching, ict, solvable_graph,
-            free_eqs, BitSet(free_vars), isder)
+        tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, free_eqs,
+            BitSet(free_vars), isder)
     end
     return var_eq_matching, full_var_eq_matching, var_sccs
 end
-
-function tear_block!(vargraph, vars,
-        var_eq_matching, ict, solvable_graph, ieqs,
-        filtered_vars, isder)
-    tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, ieqs,
-        filtered_vars,
-        isder)
-
-    # clear cache
-    vargraph.ne = 0
-    for var in vars
-        vargraph.matching[var] = unassigned
-    end
-    empty!(ieqs)
-    empty!(filtered_vars)
-end

src/structural_transformation/tearing.jl

@@ -68,3 +68,16 @@ function algebraic_variables_scc(state::TearingState)
 
     return var_eq_matching, var_sccs
 end
+
+function free_equations(graph, vars_scc, var_eq_matching, varfilter::F) where {F}
+    ne = nsrcs(graph)
+    seen_eqs = falses(ne)
+    for vars in vars_scc, var in vars
+        varfilter(var) || continue
+        ieq = var_eq_matching[var]
+        if ieq isa Int
+            seen_eqs[ieq] = true
+        end
+    end
+    findall(!, seen_eqs)
+end

test/nonlinearsystem.jl

@@ -266,3 +266,11 @@ alg_eqs = [0 ~ p - d * X]
 sys = @test_nowarn NonlinearSystem(alg_eqs; name = :name)
 @test isequal(only(unknowns(sys)), X)
 @test all(isequal.(parameters(sys), [p, d]))
+
+# Over-determined sys
+@variables u1 u2
+@parameters u3 u4
+eqs = [u3 ~ u1 + u2, u4 ~ 2 * (u1 + u2), u3 + u4 ~ 3 * (u1 + u2)]
+@named ns = NonlinearSystem(eqs, [u1, u2], [u3, u4])
+sys = structural_simplify(ns; fully_determined = false)
+@test length(unknowns(sys)) == 1