Automated sparse analytical example broken

[hersle]

https://docs.sciml.ai/ModelingToolkit/dev/examples/sparse_jacobians/

• Warning: Did not converge after maxiters = 0 substitutions. Either there is a cycle in the rules or maxiters needs to be higher.
• Sparse problem is slower than dense problem.

[ChrisRackauckas]

using BenchmarkTools
@btime solve(prob, FBDF(), save_everystep = false); # 3.248 s (2324 allocations: 65.07 MiB)
@btime solve(sparseprob, FBDF(), save_everystep = false); # 414.998 ms (44005 allocations: 407.14 MiB)
@btime solve(prob, save_everystep = false); # 7.490 s (817537 allocations: 442.43 MiB)
@btime solve(sparseprob, save_everystep = false); # 7.333 s (162922 allocations: 611.78 MiB)
@btime solve(sparseprob, FBDF(linsolve = KrylovJL_GMRES()), save_everystep = false); # 1.414 s (63738 allocations: 121.54 MiB)
@btime solve(prob, FBDF(linsolve = KrylovJL_GMRES()), save_everystep = false); # 990.209 ms (113341 allocations: 123.11 MiB)


solve(prob, save_everystep = false).alg.algs[6]
# FBDF(; max_order = Val{5}(), linsolve = KrylovJL{typeof(Krylov.gmres!), Int64, Nothing, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}(Krylov.gmres!, 0, 0, nothing, (), Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}()), nlsolve = NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}(1//100, 10, 1//5, 1//5, false, true, nothing), precs = DEFAULT_PRECS, κ = nothing, tol = nothing, extrapolant = linear, controller = Standard, step_limiter! = trivial_limiter!, autodiff = AutoFiniteDiff(),)

The key is that the default method is now (it wasn't a long time before) defaulting to using a GMRES linear solver at this size, which is matrix-free and thus doesn't use the sparse matrix much. When the matrix is actually used, you see the sparse direct method is the fastest.

This actually brings up a bigger point that it would be really good if this could automatically construct a preconditioner, using iLU or AMG, since it has the sparse matrix so it has the ability to do so. But there are some API issues since it's not part of the problem specification. @oscardssmith have you thought about this?

[oscardssmith]

no, but now is a good time to. Why does it not being part of the problem make this difficult?

[ChrisRackauckas]

Because MTK builds a problem but doesn't pass an algorithm

[oscardssmith]

Why does this have to be on the MTK side? Can't this be in the solver algorithm?

[ChrisRackauckas]

It requires using the sparse matrix for the definition