fix: fix simplification of underdetermined systems

Author: AayushSabharwal

lib/ModelingToolkitTearing/src/reassemble.jl

@@ -503,7 +503,7 @@ function generate_system_equations!(state::TearingState, neweqs::Vector{Equation
     offset = 1
     findnextfn = let diff_vars_set = diff_vars_set, solved_vars_set = solved_vars_set,
         diff_to_var = diff_to_var, ispresent = ispresent
-        j -> !(j in diff_vars_set || j in solved_vars_set) && diff_to_var[j] === nothing &&
+        j -> !(j in diff_vars_set || j in solved_vars_set || j in extra_vars) && diff_to_var[j] === nothing &&
             ispresent(j)
     end
     for (i, v) in enumerate(var_ordering)
@@ -980,7 +980,7 @@ Update the system equations, unknowns, and observables after simplification.
 function update_simplified_system!(
         state::TearingState, neweqs::Vector{Equation}, solved_eqs::Vector{Equation},
         dummy_sub::Dict{SymbolicT, SymbolicT}, var_sccs::Vector{Vector{Int}},
-        extra_unknowns::Vector{SymbolicT}, iv::Union{SymbolicT, Nothing},
+        iv::Union{SymbolicT, Nothing},
         D::Union{Differential, Shift, Nothing}; array_hack = true)
     (; fullvars, structure) = state
     (; solvable_graph, var_to_diff, eq_to_diff, graph) = structure
@@ -1027,7 +1027,6 @@ function update_simplified_system!(
     end
     unknown_idxs = filter(filterer, eachindex(state.fullvars))
     unknowns = state.fullvars[unknown_idxs]
-    unknowns = [unknowns; extra_unknowns]
     if StateSelection.is_only_discrete(structure)
         # Algebraic variables are shifted forward by one, so we backshift them.
         _unknowns = SymbolicT[]
@@ -1156,12 +1155,18 @@ function (alg::DefaultReassembleAlgorithm)(state::TearingState,
         inline_linear_sccs = false
     end
 
-    extra_unknowns = state.fullvars[extra_eqs_vars[2]]
     if StateSelection.is_only_discrete(state.structure)
         var_sccs = add_additional_history!(
             state, var_eq_matching, full_var_eq_matching, var_sccs, iv::SymbolicT)
     end
 
+    # Previously, we used to pass `state.fullvars[extra_eqs_vars[2]]` as `extra_unknowns`
+    # to `update_simplified_system!`. These were then appended to the list of unknowns.
+    # However, this is now unnecessary. `generate_system_equations!` uses a `setdiff` to
+    # populate `var_ordering` with all variables that should be in unknowns, including
+    # the extras. This is also required so that `reorder_vars!` subsets `state` correctly
+    # and the `singular_check` in initialization works.
+
     # Structural simplification
     if iv isa SymbolicT # Without iv we don't have derivatives
         D = D::Union{Differential, Shift}
@@ -1185,7 +1190,7 @@ function (alg::DefaultReassembleAlgorithm)(state::TearingState,
         # var_eq_matching and full_var_eq_matching are now invalidated
 
         sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_sccs,
-            extra_unknowns, iv, D; array_hack)
+            iv, D; array_hack)
     else
         D = D::Nothing
         neweqs, solved_eqs,
@@ -1201,7 +1206,7 @@ function (alg::DefaultReassembleAlgorithm)(state::TearingState,
         # var_eq_matching and full_var_eq_matching are now invalidated
 
         sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_sccs,
-            extra_unknowns, iv, D; array_hack)
+            iv, D; array_hack)
     end
 
     @set! state.sys = sys