Merge pull request #2203 from Keno/kf/tearingsingle

tearing: Heuristically attempt to assign equations with single solvab…

Author: YingboMa

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -32,17 +32,22 @@ function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int}
             end
         end
     end
-    for eq in es  # iterate only over equations that are not in eSolvedFixed
-        vs = Gsolvable[eq]
-        if check_der
-            # if there're differentiated variables, then only consider them
-            try_assign_eq!(ict, vs, v_active, eq, isder)
-            if has_der[]
-                has_der[] = false
-                continue
+    # Heuristic: As a first pass, try to assign any equations that only have one
+    # solvable variable.
+    for only_single_solvable in (true, false)
+        for eq in es  # iterate only over equations that are not in eSolvedFixed
+            vs = Gsolvable[eq]
+            ((length(vs) == 1) ⊻ only_single_solvable) && continue
+            if check_der
+                # if there're differentiated variables, then only consider them
+                try_assign_eq!(ict, vs, v_active, eq, isder)
+                if has_der[]
+                    has_der[] = false
+                    continue
+                end
             end
+            try_assign_eq!(ict, vs, v_active, eq)
         end
-        try_assign_eq!(ict, vs, v_active, eq)
     end
 
     return ict

test/nonlinearsystem.jl

@@ -113,10 +113,10 @@ D = Differential(t)
 @named subsys = convert_system(ODESystem, lorenz1, t)
 @named sys = ODESystem([D(subsys.x) ~ subsys.x + subsys.x], t, systems = [subsys])
 sys = structural_simplify(sys)
-u0 = [subsys.x => 1, subsys.z => 2.0]
+u0 = [subsys.x => 1, subsys.z => 2.0, subsys.y => 1.0]
 prob = ODEProblem(sys, u0, (0, 1.0), [subsys.σ => 1, subsys.ρ => 2, subsys.β => 3])
-sol = solve(prob, Rodas5())
-@test sol[subsys.x] + sol[subsys.y] - sol[subsys.z] ≈ sol[subsys.u]
+sol = solve(prob, FBDF(), reltol = 1e-7, abstol = 1e-7)
+@test sol[subsys.x] + sol[subsys.y] - sol[subsys.z]≈sol[subsys.u] atol=1e-7
 @test_throws ArgumentError convert_system(ODESystem, sys, t)
 
 @parameters t σ ρ β