ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

Automated sparse analytical example broken

Open hersle opened this issue 6 months ago • 7 comments

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.

hersle avatar Jun 29 '25 08:06 hersle

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?

ChrisRackauckas avatar Jul 15 '25 23:07 ChrisRackauckas

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

oscardssmith avatar Jul 16 '25 01:07 oscardssmith

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

ChrisRackauckas avatar Jul 16 '25 02:07 ChrisRackauckas

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

oscardssmith avatar Jul 16 '25 02:07 oscardssmith

It requires using the sparse matrix for the definition

ChrisRackauckas avatar Jul 16 '25 02:07 ChrisRackauckas

Is the automatic preconditioner in the works? It would be really nice.

hersle avatar Sep 21 '25 09:09 hersle

Preconditioners can't really be automatic. The closest is probably just iLU which just needs a user-set tolerance.

There's some research we'll be doing soon on automation of preconditioners but that's some really speculative stuff.

ChrisRackauckas avatar Sep 29 '25 06:09 ChrisRackauckas