Fix tearing_reassemble for under-determined systems (#2587)

* Fix `tearing_reassemble` for under-determined systems

* Fix minor bug

* Add tests

* Update tests

Author: YingboMa

src/structural_transformation/symbolics_tearing.jl

@@ -217,10 +217,18 @@ function check_diff_graph(var_to_diff, fullvars)
 end
 =#
 
-function tearing_reassemble(state::TearingState, var_eq_matching;
-        simplify = false, mm = nothing)
+function tearing_reassemble(state::TearingState, var_eq_matching,
+        full_var_eq_matching = nothing; simplify = false, mm = nothing)
     @unpack fullvars, sys, structure = state
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
+    extra_vars = Int[]
+    if full_var_eq_matching !== nothing
+        for v in 𝑑vertices(state.structure.graph)
+            eq = full_var_eq_matching[v]
+            eq isa Int && continue
+            push!(extra_vars, v)
+        end
+    end
 
     neweqs = collect(equations(state))
     # Terminology and Definition:
@@ -532,6 +540,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching;
     eq_to_diff = new_eq_to_diff
     diff_to_var = invview(var_to_diff)
 
+    old_fullvars = fullvars
     @set! state.structure.graph = complete(graph)
     @set! state.structure.var_to_diff = var_to_diff
     @set! state.structure.eq_to_diff = eq_to_diff
@@ -543,9 +552,15 @@ function tearing_reassemble(state::TearingState, var_eq_matching;
 
     sys = state.sys
     @set! sys.eqs = neweqs
-    @set! sys.unknowns = Any[v
-                             for (i, v) in enumerate(fullvars)
-                             if diff_to_var[i] === nothing && ispresent(i)]
+    unknowns = Any[v
+                   for (i, v) in enumerate(fullvars)
+                   if diff_to_var[i] === nothing && ispresent(i)]
+    if !isempty(extra_vars)
+        for v in extra_vars
+            push!(unknowns, old_fullvars[v])
+        end
+    end
+    @set! sys.unknowns = unknowns
     @set! sys.substitutions = Substitutions(subeqs, deps)
 
     obs_sub = dummy_sub
@@ -570,19 +585,7 @@ end
 function tearing(state::TearingState; kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
-    @unpack graph = state.structure
-    algvars = BitSet(findall(v -> isalgvar(state.structure, v), 1:ndsts(graph)))
-    aeqs = algeqs(state.structure)
-    var_eq_matching′, = tear_graph_modia(state.structure;
-        varfilter = var -> var in algvars,
-        eqfilter = eq -> eq in aeqs)
-    var_eq_matching = Matching{Union{Unassigned, SelectedState}}(var_eq_matching′)
-    for var in 1:ndsts(graph)
-        if isdiffvar(state.structure, var)
-            var_eq_matching[var] = SelectedState()
-        end
-    end
-    var_eq_matching
+    tearing_with_dummy_derivatives(state.structure, ())
 end
 
 """
@@ -594,8 +597,9 @@ instead, which calls this function internally.
 """
 function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing,
         simplify = false, kwargs...)
-    var_eq_matching = tearing(state)
-    invalidate_cache!(tearing_reassemble(state, var_eq_matching; mm, simplify))
+    var_eq_matching, full_var_eq_matching = tearing(state)
+    invalidate_cache!(tearing_reassemble(
+        state, var_eq_matching, full_var_eq_matching; mm, simplify))
 end
 
 """

test/components.jl

@@ -101,6 +101,10 @@ let
     prob2 = ODEProblem(sys2, [source.p.i => 0.0], (0, 10.0), guesses = u0)
     sol2 = solve(prob2, Rosenbrock23())
     @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ -sol2[capacitor.i]
+
+    prob3 = ODEProblem(sys2, [], (0, 10.0), guesses = u0)
+    sol3 = solve(prob2, Rosenbrock23())
+    @test sol3[unknowns(rc_model2), end] ≈ sol2[unknowns(rc_model2), end]
 end
 
 # Outer/inner connections

test/structural_transformation/tearing.jl

@@ -98,7 +98,7 @@ end
 # e5 [  1           1 |  1 ]
 
 let state = TearingState(sys)
-    torn_matching = tearing(state)
+    torn_matching, = tearing(state)
     S = StructuralTransformations.reordered_matrix(sys, torn_matching)
     @test S == [1 0 0 0 1
                 1 1 0 0 0