feat: use SCCNonlinearProblem for initialization

Author: AayushSabharwal

Project.toml

@@ -44,6 +44,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
+SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
 Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -126,7 +127,8 @@ REPL = "1"
 RecursiveArrayTools = "3.26"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "2.64"
+SCCNonlinearSolve = "1.0.0"
+SciMLBase = "2.65"
 SciMLStructures = "1.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"

src/ModelingToolkit.jl

@@ -50,6 +50,7 @@ using Distributed
 import JuliaFormatter
 using MLStyle
 using NonlinearSolve
+import SCCNonlinearSolve
 using Reexport
 using RecursiveArrayTools
 import Graphs: SimpleDiGraph, add_edge!, incidence_matrix

src/systems/diffeqs/abstractodesystem.jl

@@ -1301,6 +1301,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
         initialization_eqs = [],
         fully_determined = nothing,
         check_units = true,
+        use_scc = true,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
@@ -1318,8 +1319,8 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
     end
 
     ts = get_tearing_state(isys)
-    if warn_initialize_determined &&
-       (unassigned_vars = StructuralTransformations.singular_check(ts); !isempty(unassigned_vars))
+    unassigned_vars = StructuralTransformations.singular_check(ts)
+    if warn_initialize_determined && !isempty(unassigned_vars)
         errmsg = """
         The initialization system is structurally singular. Guess values may \
         significantly affect the initial values of the ODE. The problematic variables \
@@ -1381,9 +1382,11 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem,
                      end
         for (k, v) in u0map)
     end
-    if neqs == nunknown
-        NonlinearProblem(isys, u0map, parammap; kwargs...)
+
+    TProb = if neqs == nunknown && isempty(unassigned_vars)
+        use_scc && neqs > 0 && is_split(isys) ? SCCNonlinearProblem : NonlinearProblem
     else
-        NonlinearLeastSquaresProblem(isys, u0map, parammap; kwargs...)
+        NonlinearLeastSquaresProblem
     end
+    TProb(isys, u0map, parammap; kwargs...)
 end

src/systems/nonlinear/nonlinearsystem.jl

@@ -692,7 +692,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map,
         push!(subprobs, prob)
     end
 
-    return SCCNonlinearProblem(subprobs, explicitfuns)
+    return SCCNonlinearProblem(subprobs, explicitfuns, sys, p)
 end
 
 """

src/systems/problem_utils.jl

@@ -541,6 +541,8 @@ Keyword arguments:
   Only applicable if `warn_cyclic_dependency == true`.
 - `substitution_limit`: The number times to substitute initial conditions into each
   other to attempt to arrive at a numeric value.
+- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully
+  determined.
 
 All other keyword arguments are passed as-is to `constructor`.
 """
@@ -554,7 +556,7 @@ function process_SciMLProblem(
         symbolic_u0 = false, warn_cyclic_dependency = false,
         circular_dependency_max_cycle_length = length(all_symbols(sys)),
         circular_dependency_max_cycles = 10,
-        substitution_limit = 100, kwargs...)
+        substitution_limit = 100, use_scc = true, kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys)
     iv = has_iv(sys) ? get_iv(sys) : nothing
@@ -607,7 +609,7 @@ function process_SciMLProblem(
                 sys, t, u0map, pmap; guesses, warn_initialize_determined,
                 initialization_eqs, eval_expression, eval_module, fully_determined,
                 warn_cyclic_dependency, check_units = check_initialization_units,
-                circular_dependency_max_cycle_length, circular_dependency_max_cycles)
+                circular_dependency_max_cycle_length, circular_dependency_max_cycles, use_scc)
             initializeprobmap = getu(initializeprob, unknowns(sys))
 
             punknowns = [p