Minor fixes

Author: YingboMa

src/systems/diffeqs/modelingtoolkitize.jl

@@ -23,13 +23,17 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     D = Differential(t)
     mm = prob.f.mass_matrix
 
-    lhs = map(mm * vars) do v
-        if iszero(v)
-            0
-        elseif v in var_set
-            D(v)
-        else
-            error("Non-permuation mass matrix is not supported.")
+    if mm === I
+        lhs = map(v->D(v), vars)
+    else
+        lhs = map(mm * vars) do v
+            if iszero(v)
+                0
+            elseif v in var_set
+                D(v)
+            else
+                error("Non-permuation mass matrix is not supported.")
+            end
         end
     end
 
@@ -42,17 +46,17 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
 
     eqs = vcat([lhs[i] ~ rhs[i] for i in eachindex(prob.u0)]...)
 
-    sts = Vector(vec(vars))
+    sts = vec(collect(vars))
     params = if ndims(params) == 0
         [params[1]]
     else
-        Vector(vec(params))
+        vec(collect(params))
     end
 
     de = ODESystem(
         eqs, t, sts, params,
-        default_u0=Dict(sts .=> prob.u0),
-        default_p=Dict(params .=> prob.p)
+        default_u0=Dict(sts .=> vec(collect(prob.u0))),
+        default_p=Dict(params .=> vec(collect(prob.p))),
     )
 
     de

test/modelingtoolkitize.jl

@@ -162,6 +162,7 @@ sol = solve(prob, Vern8())
 modelingtoolkitize(prob)
 
 # Index reduction and mass matrix handling
+using LinearAlgebra
 function pendulum!(du, u, p, t)
     x, dx, y, dy, T = u
     g, L = p
@@ -177,10 +178,10 @@ u0 = [1.0, 0, 0, 0, 0]
 p = [9.8, 1]
 tspan = (0, 10.0)
 pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)
-pendulum_sys  = dae_index_lowering(modelingtoolkitize(pendulum_prob))
+pendulum_sys_org = modelingtoolkitize(pendulum_prob)
+sts = states(pendulum_sys_org)
+pendulum_sys = dae_index_lowering(pendulum_sys_org)
 prob = ODEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Rodas4())
-@parameters t
-@variables x[1:5](t)
-l2 = sol[x[1]].^2 + sol[x[3]].^2
+l2 = sol[sts[1]].^2 + sol[sts[3]].^2
 @test all(l->abs(sqrt(l) - 1) < 0.05, l2)