Fix Wfact_t for DAEs

Author: YingboMa

src/systems/diffeqs/abstractodesystem.jl

@@ -63,16 +63,17 @@ function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
     isempty(sys.Wfact[]) || return (sys.Wfact[],sys.Wfact_t[])
 
     jac = calculate_jacobian(sys)
+    M = calculate_massmatrix(sys)
     gam = Variable(:__MTKWgamma)()
 
-    W = - LinearAlgebra.I + gam*jac
+    W = - M + gam*jac
     Wfact = lu(W, Val(false), check=false).factors
 
     if simplify
         Wfact = ModelingToolkit.simplify.(Wfact)
     end
 
-    W_t = - LinearAlgebra.I/gam + jac
+    W_t = - M/gam + jac
     Wfact_t = lu(W_t, Val(false), check=false).factors
     if simplify
         Wfact_t = ModelingToolkit.simplify.(Wfact_t)

test/odesystem.jl

@@ -212,9 +212,8 @@ p  = [k₁ => 0.04,
 tspan = (0.0,100000.0)
 prob1 = ODEProblem(sys,u0,tspan,p)
 prob2 = ODEProblem(sys,u0,tspan,p,jac=true)
-# Wfact version is not very stable because of the lack of pivoting
 prob3 = ODEProblem(sys,u0,tspan,p,Wfact=true,Wfact_t=true)
-for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 0.1)]
+for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)]
     sol = solve(prob, Rodas5())
     @test all(x->≈(sum(x), 1.0, atol=atol), sol.u)
 end