Sparse jacobians not working #3529
Description
Also raised here: SciML/OrdinaryDiffEq.jl#2654 (comment), but might be an MTK thing:
# Fetch packages
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
# Creates a model.
@parameters p d k1 k2 k3 k4 k5 k6 k7
@variables X(t) Y(t) Z(t) XY(t) XZ2(t) Y3(t) X2(t) Y3Z2(t) XYZ(t) XY3Z4(t) X2Y2Z2(t) V(t) V2(t)
eqs = [
D(X) ~ p - d*X - k1*Y*X - (1//2)*k2*X*(Z^2) - k4*Y*X*Z
D(Y) ~ p - d*Y - k1*Y*X - k4*Y*X*Z
D(Z) ~ p - d*Z - k2*X*(Z^2) - k4*Y*X*Z
D(XY) ~ k1*Y*X
D(XZ2) ~ -k5*Y3Z2*XZ2 + (1//2)*k2*X*(Z^2)
D(Y3) ~ -k3*Y3*X2
D(X2) ~ -k3*Y3*X2
D(Y3Z2) ~ k3*Y3*X2 - k5*Y3Z2*XZ2
D(XYZ) ~ -k6*(XYZ^2) + k4*Y*X*Z
D(XY3Z4) ~ -d*XY3Z4 + k5*Y3Z2*XZ2
D(X2Y2Z2) ~ -d*X2Y2Z2 + (1//2)*k6*(XYZ^2)
D(V) ~ (-k7*(V^2)) / (1 + t) - (Y + X)*V
D(V2) ~ ((1//2)*k7*(V^2)) / (1 + t) - Z*V2
]
@mtkbuild osys = ODESystem(eqs, t)
# Creates sparse/dense Jacobians
u0 = [u => i/10 for (i,u) in enumerate(unknowns(osys))]
ps = [p => i/10 for (i,p) in enumerate(parameters(osys))]
oprob_jac = ODEProblem(osys, u0, 1.0, ps; jac = true, sparse = false)
oprob_sjac = ODEProblem(osys, u0, 1.0, ps; jac = true, sparse = true)
# Checks whether they are the same (they are not).
function get_jac(prob, sparse)
J = sparse ? deepcopy(prob.f.jac_prototype) : zeros(length(prob.u0), length(prob.u0))
ModelingToolkit.is_time_dependent(prob) ? prob.f.jac(J, prob.u0, prob.p, 0.0) : prob.f.jac(J, prob.u0, prob.p)
return J
end
jac = get_jac(oprob_jac, false)
jac_sparse = get_jac(oprob_sjac, true)
Matrix(jac_sparse) == jac # false
maximum(abs2.(Matrix(jac_sparse) .- jac)) # 2.6244000000000005