fix: fix several problems raised in multibody-related issues

Project.toml

@@ -78,11 +78,11 @@ Graphs = "1.5.2"
 InteractiveUtils = "1"
 Libdl = "1"
 LinearAlgebra = "1"
-LinearSolve = "3.19.2"
+LinearSolve = "3.56"
 Logging = "1"
 ModelingToolkitBase = "1.4"
 ModelingToolkitStandardLibrary = "2.20"
-ModelingToolkitTearing = "1.1.0"
+ModelingToolkitTearing = "1.2.1"
 Moshi = "0.3.6"
 NonlinearSolve = "4.3"
 OffsetArrays = "1"
@@ -96,8 +96,8 @@ PrecompileTools = "1.2.1"
 REPL = "1"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.12"
-SCCNonlinearSolve = "1.4.0"
-SciMLBase = "2.125.0"
+SCCNonlinearSolve = "1.8.1"
+SciMLBase = "2.132"
 SciMLPublic = "1.0.0"
 Serialization = "1"
 Setfield = "0.7, 0.8, 1"

lib/ModelingToolkitBase/src/problems/linearproblem.jl

@@ -69,7 +69,7 @@ function SciMLBase.LinearProblem{iip}(
         kwargs...
     )
 
-    if any(x -> symbolic_type(x) != NotSymbolic(), u0)
+    if any(x -> symbolic_type(x) != NotSymbolic() || x === nothing, u0)
         u0 = nothing
     end
 
@@ -110,7 +110,15 @@ function get_A_b_from_LinearFunction(
         A = u0_constructor(u0_eltype.(get_A(p)))
         b = u0_constructor(u0_eltype.(get_b(p)))
     else
-        A = u0_constructor(u0_eltype.(interface.update_A!(p)))
+        A = u0_eltype.(interface.update_A!(p))
+        szA = size(A)::NTuple{2, Int}
+        _A = u0_constructor(A)
+        if ArrayInterface.ismutable(_A)
+            A = similar(_A, szA)
+            copyto!(A, _A)
+        else
+            A = StaticArraysCore.similar_type(_A, StaticArraysCore.Size(szA))(_A)
+        end
         b = u0_constructor(u0_eltype.(interface.update_b!(p)))
     end
     if sparse
@@ -125,6 +133,11 @@ function SciMLBase.get_new_A_b(
         sys::AbstractSystem, f::SciMLBase.SymbolicLinearInterface, p, A, b; kw...
     )
     if ArrayInterface.ismutable(A)
+        T = eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1])
+        if eltype(A) !== T
+            A = similar(A, T)
+            b = similar(b, T)
+        end
         f.update_A!(A, p)
         f.update_b!(b, p)
     else

lib/ModelingToolkitBase/src/systems/abstractsystem.jl

@@ -566,7 +566,7 @@ function add_initialization_parameters(sys::AbstractSystem; split = true)
     _sys = unhack_system(sys)
     obs = observed(_sys)
     eqs = equations(_sys)
-    for x in unknowns(sys)
+    for x in unknowns(_sys)
         if !split
             push!(all_initialvars, x)
             continue
@@ -2669,9 +2669,14 @@ end
 
 function default_to_parentscope(v)
     uv = unwrap(v)
-    uv isa SymbolicT || return v
-    return apply_to_variables(v) do sym
-        ParentScope(sym)
+    if uv isa SymbolicT
+        return apply_to_variables(v) do sym
+            ParentScope(sym)
+        end
+    elseif is_array_of_symbolics(uv)
+        return map(default_to_parentscope, v)
+    else
+        return v
     end
 end
 

lib/ModelingToolkitBase/test/initial_values.jl

@@ -204,7 +204,7 @@ end
         pend, [x => 1, g => 1], (0, 1); guesses = [y => λ, λ => y + 1]
     )
     ODEProblem(
-        pend, [x => 1, g => 1], (0, 1);
+        pend, [x => 1, D(y) => 0, g => 1], (0, 1);
         guesses = [y => λ, λ => 0.5], missing_guess_value
     )
 
@@ -346,10 +346,10 @@ end
 
 @testset "`p_constructor` keyword argument" begin
     @parameters g = 1.0
-    @variables x(t) y(t) [state_priority = 10, guess = 1.0] λ(t) [guess = 1.0]
+    @variables x(t) y(t) [state_priority = 10] λ(t)
     @named pend = index_reduced_pend()
     pend = complete(pend)
-    guesses(pend)[y] = 1.0
+    guesses(pend)[x] = 1.0
     guesses(pend)[λ] = 1.0
     initial_conditions(pend)[g] = 1.0
 
@@ -365,7 +365,10 @@ end
     @test parameter_values(initdata.initializeprob).tunable isa SVector
 
     pend = complete(pend; split = false)
-    prob = ODEProblem(pend, u0, tspan; u0_constructor, p_constructor, missing_guess_value)
+    prob = ODEProblem(
+        pend, u0, tspan; u0_constructor, p_constructor, missing_guess_value,
+        guesses = [D(x) => 1.0]
+    )
     @test prob.p isa SVector
     initdata = prob.f.initialization_data
     @test state_values(initdata.initializeprob) isa SVector

lib/ModelingToolkitBase/test/variable_scope.jl

@@ -1,6 +1,7 @@
 using ModelingToolkitBase
 using ModelingToolkitBase: SymScope, t_nounits as t, D_nounits as D
-using Symbolics: arguments, value, getname
+using Symbolics: arguments, value, getname, VartypeT
+import SymbolicUtils as SU
 using Test
 
 @variables a b(t) c d e(t)
@@ -157,3 +158,16 @@ end
     sys2 = ModelingToolkitBase.discover_globalscoped(sys)
     @test isequal(only(parameters(sys2)), cond)
 end
+
+@testset "`@named` applies `ParentScope` to arrays of symbolics" begin
+    function Foo(; name, k)
+        tmp = SU.Const{VartypeT}(k)
+        vars = SU.search_variables(tmp)
+        for var in vars
+            @test SU.getmetadata(var, SymScope, LocalScope()) == ParentScope(LocalScope())
+        end
+    end
+    @parameters k
+    @variables x(t)
+    @named foo = Foo(; k = [k + 1, x + 2])
+end

src/problems/sccnonlinearproblem.jl

@@ -119,9 +119,17 @@ function SciMLBase.SCCNonlinearProblem{iip}(
     end
 
     if length(var_sccs) == 1
-        return NonlinearProblem{iip}(
-            sys, op; eval_expression, eval_module, kwargs...
-        )
+        if isaffine(sys)
+            linprob = LinearProblem{iip}(
+                sys, op; eval_expression, eval_module,
+                u0_constructor, cse, kwargs...
+            )
+            return SCCNonlinearProblem((linprob,), (Returns(nothing),), parameter_values(linprob), true; sys)
+        else
+            return NonlinearProblem{iip}(
+                sys, op; eval_expression, eval_module, u0_constructor, cse, kwargs...
+            )
+        end
     end
 
     dvs = unknowns(sys)

test/runtests.jl

@@ -64,6 +64,7 @@ end
         @mtktestset("Guess Propagation", "guess_propagation.jl")
         @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl")
         @mtktestset("InitializationSystem Test", "initializationsystem.jl")
+        @mtktestset("Initial Values Test", "initial_values.jl")
     end
 
     if GROUP == "All" || GROUP == "InterfaceII"

test/structural_transformation/tearing.jl

@@ -219,41 +219,41 @@ prob_complex = ODEProblem(sys, u0, (0, 1.0))
 sol = solve(prob_complex, Tsit5())
 @test all(sol[mass.v] .== 1)
 
-#=
 using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Blocks: Constant
 
-@testset "Inline linear SCCs" begin
-    function RCModel(; name)
-        pars = @parameters begin
-            R = 1.0
-            C = 1.0
-            V = 1.0
-        end
-        systems = @named begin
-            resistor1 = Resistor(R = R)
-            resistor2 = Resistor(R = R)
-            capacitor = Capacitor(C = C, v = 0.0)
-            source = Voltage()
-            constant = Constant(k = V)
-            ground = Ground()
-        end
-        eqs = [
-            connect(constant.output, source.V)
-            connect(source.p, resistor1.p)
-            connect(resistor1.n, resistor2.p)
-            connect(resistor2.n, capacitor.p)
-            connect(capacitor.n, source.n, ground.g)
-        ]
-        return System(eqs, t, [], pars; systems, name)
+function RCModel(; name)
+    pars = @parameters begin
+        R = 1.0
+        C = 1.0
+        V = 1.0
+    end
+    systems = @named begin
+        resistor1 = Resistor(R = R)
+        resistor2 = Resistor(R = R)
+        capacitor = Capacitor(C = C, v = 0.0)
+        source = Voltage()
+        constant = Constant(k = V)
+        ground = Ground()
     end
+    eqs = [
+        connect(constant.output, source.V)
+        connect(source.p, resistor1.p)
+        connect(resistor1.n, resistor2.p)
+        connect(resistor2.n, capacitor.p)
+        connect(capacitor.n, source.n, ground.g)
+    ]
+    return System(eqs, t, [], pars; systems, name)
+end
 
+
+@testset "Inline linear SCCs" begin
     reassemble_alg1 = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true)
     reassemble_alg2 = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 0)
     @mtkcompile sys1 = RCModel()
-    @mtkcompile sys2 = RCModel() reassemble_alg=reassemble_alg1
-    @mtkcompile sys3 = RCModel() reassemble_alg=reassemble_alg2
-    
+    @mtkcompile sys2 = RCModel() reassemble_alg = reassemble_alg1
+    @mtkcompile sys3 = RCModel() reassemble_alg = reassemble_alg2
+
     @test length(equations(sys1)) == 2
     @test length(equations(sys2)) == 1
     @test isequal(only(unknowns(sys2)), sys2.capacitor.v)
@@ -263,21 +263,27 @@ using ModelingToolkitStandardLibrary.Blocks: Constant
     idx = findfirst(isequal(sys3.capacitor.i), observables(sys3))
     rhs = observed(sys3)[idx].rhs
     @test operation(rhs) === getindex
-    @test operation(arguments(rhs)[1]) === solve
+    @test operation(arguments(rhs)[1]) === (\)
 
     prob1 = ODEProblem(sys1, [], (0.0, 10.0); guesses = [sys1.resistor1.v => 1.0])
     prob2 = ODEProblem(sys2, [], (0.0, 10.0))
     prob3 = ODEProblem(sys3, [], (0.0, 10.0))
 
-    sol1 = solve(prob1, Rodas5P(); abstol = 1e-8, reltol = 1e-8)
-    sol2 = solve(prob2, Tsit5(), abstol = 1e-8, reltol = 1e-8)
-    sol3 = solve(prob3, Tsit5(), abstol = 1e-8, reltol = 1e-8)
+    sol1 = solve(prob1, Rodas5P(); abstol = 1.0e-8, reltol = 1.0e-8)
+    sol2 = solve(prob2, Tsit5(), abstol = 1.0e-8, reltol = 1.0e-8)
+    sol3 = solve(prob3, Tsit5(), abstol = 1.0e-8, reltol = 1.0e-8)
 
     @test SciMLBase.successful_retcode(sol1)
     @test SciMLBase.successful_retcode(sol2)
     @test SciMLBase.successful_retcode(sol3)
 
-    @test sol2(sol1.t; idxs = unknowns(sys1)).u ≈ sol1.u atol=1e-8
-    @test sol3(sol1.t; idxs = unknowns(sys1)).u ≈ sol1.u atol=1e-8
+    @test sol2(sol1.t; idxs = unknowns(sys1)).u ≈ sol1.u atol = 1.0e-8
+    @test sol3(sol1.t; idxs = unknowns(sys1)).u ≈ sol1.u atol = 1.0e-8
+end
+
+@testset "`Initial` parameters are added for observed variables solved by inline linear SCCS" begin
+    reassemble_alg = StructuralTransformations.DefaultReassembleAlgorithm(; inline_linear_sccs = true, analytical_linear_scc_limit = 1)
+    @mtkcompile sys = RCModel() reassemble_alg = reassemble_alg
+    @test Initial(sys.resistor1.v) in Set(ModelingToolkit.get_ps(sys))
+    @test Initial(sys.resistor2.v) in Set(ModelingToolkit.get_ps(sys))
 end
-=#