ILU preconditioned GMRES inside TRBDF2 for ODEs generated by NetworkDynamics.jl

General Usage Performance
#1

We are solving differential equations defined on a graph using NetworkDynamics.jl and DifferentialEquations.jl. In this context we have the two questions outlined below.

Thx. Domenico Lahaye.

1/ Switching the automatic Jacobian off in NetworkDynamics.jl

We would like to convince ourselves that using the Jacobian indeed makes sense by comparing the CPU time of a run with and without the Jacobian. We, however, do not manage to switch off the automatic computation (or use) of the Jacobian. That is, the two runs outlined below give identical CPU time. What are we doing wrong or overseeing?

(ModelingToolkit does not allow to switch off Jacobian computation? Use of Jacobian should be switched off as option for time-integration method? What are solver options for TRBDF2 aside from linsolve= … ?)

RUN 1: prob = ODEProblem(network!, u0, tspan, p, jac=true)

RUN 2: prob = ODEProblem(network!, u0, tspan, p, jac=false)

2/ Specifying ILU preconditioned GMRES as solver within TRBDF2 in NetworkDynamics.jl

We managed to run TRBDF2 with unpreconditioned GMRES as Krylov method. Next we would like to switch on ILU as a preconditioner. The use of ILU requires the detection of the sparsity pattern. The sparsity detection conflicts with the definition by NetworkDynamics of ODEFuntion. Details (including the error message) are given at SIR.jl/preconditioner.jl at 8a7a9a8bb0fdbcc615030a2da755461a39a1d266 · bszhu/SIR.jl · GitHub

(network! = network_dynamics(node, edge, g_directed) # key constructor network_dynamics returns an ODEFunction)

(Should we call the function generate_jacobian(network!,output,input) instead? What does aspreconditioner refer to?)

#2

Any suggestions?

#3

Stiff ODEs use the Jacobian unless you change the jac_prototype to be a matrix-free object, which would then require a compatible linear solver. This is covered in the tutorial on stiff ODEs:

https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/

though admittedly we could (and will) make that a little nicer. So if you’re still using a direct solver, jac=true is just generating the analytical solution to the Jacobian and generating the code for that, vs using AD. In some cases that will have similar performance. The advantage vs AD usually is only clear when sparsity is involved, since even sparse automatic differentiation when done optimally has some limitations.

Since you have MTK working, you can use its sparsejacobian function to extract the sparse Jacobian, or simply add sparse=true in the generation phase and grab the jac_prototype (there it would also have the advantage of specializing the Jacobian computation on the sparsity pattern)

#4

Sincere thx! We will give it a look and get back here.

#5

Spending more time on the documentation proved valuable :wink: Thank you for that suggestion.

The page ODE Solvers · DifferentialEquations.jl list TRBDF2 as an SDIRK method. Is this intentional?

#6

TRBDF2 is an SDIRK method, yes. In fact, it is in the class of eSDIRK methods.

#7

Will have to read up on this. Thx again.