Sparse Jacobian with dense blocks in differentialEquations.jl or otherwise

ziolai · May 2, 2021, 6:51pm

Dear all,

I am planning to use differentialEquations.jl to solve a stiff system of ordinary differential equations in which the Jacobian is sparse with dense blocks. Such systems are well known in this community as they appear in networks with dynamics (for e.g. biology to chemistry) prescribed on the nodes of the network.

Here is what I am wondering about.

1/ Non-Linear Solver at Each Time Step

1.1/ Can the performance of the non-linear solver at each time step be monitored (equivalent of -snes_monitor in PETSc e.g.)

1.2/ Can the setting of the non-linear solver (tolerances, Jacobians approximations) be changed

2/ Linear Solver at Each Non-Linear Step

Can the linear solver (preconditioned Krylov method) be customized and monitored?

3/ Sparsity structure of the Jacobian

The ODE has the form \dot{x} = A x + f(x) leading to a sparse Jacobian with dense blocks. Can the SplitFunction be used to specify an additive preconditioned that treats the contribution A and Jac(f(x)) separately?

Thanks, Domenico.

ChrisRackauckas · May 2, 2021, 9:24pm

Yes, the heuristics are all exposed:

github.com

SciML/OrdinaryDiffEq.jl/blob/master/src/algorithms.jl#L1754


      
          Second order stabilized Runge-Kutta method.
          Exhibits high stability for real eigenvalues and is smoothened to allow for moderate sized complex eigenvalues.
          
          This method takes optional keyword arguments `min_stages`, `max_stages`, and `eigen_est`.
          The function `eigen_est` should be of the form
          
          `eigen_est = (integrator) -> integrator.eigen_est = upper_bound`,
          
          where `upper_bound` is an estimated upper bound on the spectral radius of the Jacobian matrix. If `eigen_est`
          is not provided, `upper_bound` will be estimated using the power iteration.
          """
          struct ROCK2{E} <: OrdinaryDiffEqAdaptiveAlgorithm
              min_stages::Int
              max_stages::Int
              eigen_est::E
          end
          function ROCK2(; min_stages = 0, max_stages = 200, eigen_est = nothing)
              ROCK2(min_stages, max_stages, eigen_est)
          end
          
          """

https://github.com/SciML/DiffEqBase.jl/blob/master/src/nlsolve/type.jl#L65

though it needs more docs.

Yes, you can even write your own routine.

https://diffeq.sciml.ai/stable/features/linear_nonlinear/

You can just do that in the linear solver.

ziolai · May 3, 2021, 8:26am

Amazing! Thx.