Merge pull request #2510 from SciML/dd

Expose a non dummy derivative index lowering

src/structural_transformation/StructuralTransformations.jl

@@ -54,7 +54,7 @@ export tearing, partial_state_selection, dae_index_lowering, check_consistency
 export dummy_derivative
 export build_torn_function, build_observed_function, ODAEProblem
 export sorted_incidence_matrix,
-       pantelides!, tearing_reassemble, find_solvables!,
+       pantelides!, pantelides_reassemble, tearing_reassemble, find_solvables!,
        linear_subsys_adjmat!
 export tearing_assignments, tearing_substitution
 export torn_system_jacobian_sparsity

src/systems/systemstructure.jl

@@ -616,6 +616,7 @@ end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
+        dummy_derivative = true,
         kwargs...)
     check_consistency &= fully_determined
     has_io = io !== nothing
@@ -628,7 +629,13 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     if check_consistency
         ModelingToolkit.check_consistency(state, orig_inputs)
     end
-    if fully_determined
+    if fully_determined && dummy_derivative
+        sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
+    elseif fully_determined
+        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        sys = pantelides_reassemble(state, var_eq_matching)
+        state = TearingState(sys)
+        sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
         sys = ModelingToolkit.dummy_derivative(sys, state; simplify, mm, check_consistency)
     else
         sys = ModelingToolkit.tearing(sys, state; simplify, mm, check_consistency)

test/structural_transformation/index_reduction.jl

@@ -155,3 +155,19 @@ let sys = structural_simplify(pendulum2)
     @test sol.retcode == ReturnCode.Success
     @test norm(sol[x] .^ 2 + sol[y] .^ 2 .- 1) < 1e-2
 end
+
+let
+    @parameters g
+    @variables x(t) [state_priority = 10] y(t) λ(t)
+
+    eqs = [D(D(x)) ~ λ * x
+           D(D(y)) ~ λ * y - g
+           x^2 + y^2 ~ 1]
+    @named pend = ODESystem(eqs, t)
+    sys = complete(structural_simplify(pend; dummy_derivative = false))
+    prob = ODEProblem(
+        sys, [x => 1, y => 0, D(x) => 0.0], (0.0, 10.0), [g => 1], guesses = [λ => 0.0])
+    sol = solve(prob, Rodas5P())
+    @test SciMLBase.successful_retcode(sol)
+    @test sol[x^2 + y^2][end] < 1.1
+end