Handle mass matrix in modelingtoolkitize

Author: YingboMa

src/systems/diffeqs/modelingtoolkitize.jl

@@ -18,19 +18,29 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     vars = ArrayInterface.restructure(prob.u0,[var(:x, i)(ModelingToolkit.value(t)) for i in eachindex(prob.u0)])
     params = p isa DiffEqBase.NullParameters ? [] :
              reshape([Num(Sym{Real}(nameof(Variable(:α, i)))) for i in eachindex(p)],size(p))
+    var_set = Set(vars)
 
     D = Differential(t)
+    mm = prob.f.mass_matrix
 
-    rhs = [D(var) for var in vars]
+    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
 
     if DiffEqBase.isinplace(prob)
-        lhs = similar(vars, Num)
-        prob.f(lhs, vars, params, t)
+        rhs = similar(vars, Num)
+        prob.f(rhs, vars, params, t)
     else
-        lhs = prob.f(vars, params, t)
+        rhs = prob.f(vars, params, t)
     end
 
-    eqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)
+    eqs = vcat([lhs[i] ~ rhs[i] for i in eachindex(prob.u0)]...)
 
     params = if ndims(params) == 0
         [params[1]]