fix: better handle allow_algebraic=false

Author: AayushSabharwal

src/structural_transformation/partial_state_selection.jl

@@ -170,13 +170,14 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 end
 
 function dummy_derivative_graph!(state::TransformationState, jac = nothing;
-        state_priority = nothing, log = Val(false), allow_symbolic = false, kwargs...)
+        state_priority = nothing, log = Val(false), allow_symbolic = false, allow_algebraic = true, kwargs...)
     state.structure.solvable_graph === nothing &&
-        find_solvables!(state; allow_symbolic, kwargs...)
+        find_solvables!(state; allow_symbolic, allow_algebraic, kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state; allow_symbolic, kwargs...))
+    var_eq_matching = complete(pantelides!(
+        state; allow_symbolic, allow_algebraic, kwargs...))
     dummy_derivative_graph!(
-        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic)
+        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic, allow_algebraic)
 end
 
 struct DummyDerivativeSummary
@@ -186,7 +187,8 @@ end
 
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
-        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false) where {log}
+        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false,
+        allow_algebraic = true) where {log}
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -345,7 +347,8 @@ function dummy_derivative_graph!(
     end
 
     dummy_derivatives_set = BitSet(dummy_derivatives)
-    make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
+    make_differential_denominators_unsolvable!(
+        structure, dummy_derivatives_set; allow_algebraic)
 
     ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log

src/structural_transformation/symbolics_tearing.jl

@@ -926,7 +926,8 @@ end
 function tearing(state::TearingState; kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
-    make_differential_denominators_unsolvable!(state.structure)
+    make_differential_denominators_unsolvable!(
+        state.structure; allow_algebraic = get(kwargs, :allow_algebraic, true))
     tearing_with_dummy_derivatives(state.structure, ())
 end
 
@@ -939,7 +940,7 @@ instead, which calls this function internally.
 """
 function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing,
         simplify = false, array_hack = true, kwargs...)
-    var_eq_matching, full_var_eq_matching = tearing(state)
+    var_eq_matching, full_var_eq_matching = tearing(state; kwargs...)
     invalidate_cache!(tearing_reassemble(
         state, var_eq_matching, full_var_eq_matching; mm, simplify, array_hack))
 end

src/structural_transformation/utils.jl

@@ -284,9 +284,11 @@ differential variables that should be considered as algebraic for the purpose of
 transformation.
 """
 function make_differential_denominators_unsolvable!(
-        structure::SystemStructure, additional_algevars = ())
+        structure::SystemStructure, additional_algevars = (); allow_algebraic)
     for ((eqi, vari), denoms) in structure.denominators
-        all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) && continue
+        if allow_algebraic && all(i -> isalgvar(structure, i) || i in additional_algevars, denoms)
+            continue
+        end
         rem_edge!(structure.solvable_graph, eqi, vari)
     end
 end

src/systems/systemstructure.jl

@@ -722,18 +722,23 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
-        StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
+        var_eq_matching = pantelides!(state; finalize = false, allow_algebraic, kwargs...)
+        StructuralTransformations.make_differential_denominators_unsolvable!(
+            state.structure; allow_algebraic)
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
-        sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
+        sys, mm = ModelingToolkit.alias_elimination!(
+            state; allow_symbolic, allow_algebraic, kwargs...)
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     else
         sys = ModelingToolkit.tearing(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)