Merge pull request #1890 from Keno/kf/knotty

state_selection: Correct is_not_present condition

Author: YingboMa

src/structural_transformation/partial_state_selection.jl

@@ -267,6 +267,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     end
     dummy_derivatives_set = BitSet(dummy_derivatives)
 
+    irreducible_set = BitSet()
     if ag !== nothing
         function isreducible(x)
             # `k` is reducible if all lower differentiated variables are.
@@ -277,27 +278,46 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
                 end
                 x = diff_to_var[x]
                 x === nothing && break
+                # We deliberately do not check `isempty(𝑑neighbors(graph, x))`
+                # because when `D(x)` appears in the alias graph, and `x`
+                # doesn't appear in any equations nor in the alias graph, `D(x)`
+                # is not reducible. Consider the system `D(x) ~ 0`.
                 if !haskey(ag, x)
                     isred = false
                 end
             end
             isred
         end
-        irreducible_set = BitSet()
-        for (k, (_, v)) in ag
+        for (k, (c, v)) in ag
             isreducible(k) || push!(irreducible_set, k)
             isreducible(k) || push!(irreducible_set, k)
+            iszero(c) && continue
             push!(irreducible_set, v)
         end
     end
 
-    is_not_present = v -> isempty(𝑑neighbors(graph, v)) &&
-        (ag === nothing || (haskey(ag, v) && !(v in irreducible_set)))
+    is_not_present_non_rec = let graph = graph, irreducible_set = irreducible_set
+        v -> begin
+            not_in_eqs = isempty(𝑑neighbors(graph, v))
+            ag === nothing && return not_in_eqs
+            isirreducible = v in irreducible_set
+            return not_in_eqs && !isirreducible
+        end
+    end
+
+    is_not_present = let var_to_diff = var_to_diff
+        v -> while true
+            # if a higher derivative is present, then it's present
+            is_not_present_non_rec(v) || return false
+            v = var_to_diff[v]
+            v === nothing && return true
+        end
+    end
+
     # Derivatives that are either in the dummy derivatives set or ended up not
     # participating in the system at all are not considered differential
     is_some_diff = let dummy_derivatives_set = dummy_derivatives_set
-        v -> !(v in dummy_derivatives_set) &&
-            !(var_to_diff[v] === nothing && is_not_present(v))
+        v -> !(v in dummy_derivatives_set) && !is_not_present(v)
     end
 
     # We don't want tearing to give us `y_t ~ D(y)`, so we skip equations with
@@ -309,7 +329,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     # We can eliminate variables that are not a selected state (differential
     # variables). Selected states are differentiated variables that are not
     # dummy derivatives.
-    can_eliminate = let var_to_diff = var_to_diff
+    can_eliminate = let var_to_diff = var_to_diff, ag = ag
         v -> begin
             if ag !== nothing
                 haskey(ag, v) && return false

test/structural_transformation/tearing.jl

@@ -218,8 +218,8 @@ u0 = [mass.s => 0.0
       mass.v => 1.0]
 
 sys = structural_simplify(ms_model)
-@test sys.jac[] === ModelingToolkit.EMPTY_JAC
-@test sys.tgrad[] === ModelingToolkit.EMPTY_TGRAD
+@test ModelingToolkit.get_jac(sys)[] === ModelingToolkit.EMPTY_JAC
+@test ModelingToolkit.get_tgrad(sys)[] === ModelingToolkit.EMPTY_TGRAD
 prob_complex = ODAEProblem(sys, u0, (0, 1.0))
 sol = solve(prob_complex, Tsit5())
 @test all(sol[mass.v] .== 1)