Solving large rigid equations

This guide discusses additional possibilities for efficiently solving large, rigid ordinary differential equations. To solve rigid ordinary differential equations, it is necessary to specialize the linear solver in the properties of the Jacobian in order to reduce the number of linear solutions. and reverse solutions . The same functions and controls apply to rigid equations SDE, DDE, DAE, etc. This guide is intended for large-scale models, such as those obtained for semi-discretization of partial differential equations (PDEs). For example, we will use the Brussellator rigid differential equation (BRUSS).

This guide is intended for experienced users who want to familiarize themselves with the advanced features. DifferentialEquations.jl automates most of the usage, so users are advised to try solve(prob) with an automated algorithm first.

Definition of the Brusselator equation

You can skip this section: it simply sets out the task of the example.

The PDE equation of the Brusselator is defined as follows:

where

and the initial conditions are

with a periodic boundary condition

in the time range ].

To solve this PDE equation, we discretize it into the ODE system of equations using the finite difference method. We discretize u and v into arrays of values at each time point: u[i,j] = u(i*dx,j*dy) for some choice of dx/dy', and the same for `v'. Then the ODE equation is defined using `U[i,j,k] = [u v]. The second derivative operator, the Laplacian, is discretized and becomes a tridiagonal matrix with [1-2-1] and 1 in the upper-right and lower-left corners. Then nonlinear functions are applied at each point in space (they are translated). Use `dx=dy=1/32'.

The resulting definition of `ODEProblem' looks like this:

using DifferentialEquations, LinearAlgebra, SparseArrays

const N = 32
const xyd_brusselator = range(0, stop = 1, length = N)
brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
    A, B, alpha, dx = p
    alpha = alpha / dx^2
    @inbounds for I in CartesianIndices((N, N))
        i, j = Tuple(I)
        x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
        ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
                             limit(j - 1, N)
        du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
                       4u[i, j, 1]) +
                      B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
                      brusselator_f(x, y, t)
        du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
                       4u[i, j, 2]) +
                      A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
    end
end
p = (3.4, 1.0, 10.0, step(xyd_brusselator))

function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I, 1] = 22 * (y * (1 - y))^(3 / 2)
        u[I, 2] = 27 * (x * (1 - x))^(3 / 2)
    end
    u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p)

ODEProblem with uType Array{Float64, 3} and tType Float64. In-place: true
timespan: (0.0, 11.5)
u0: 32×32×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 ⋮                                  ⋱                      ⋮
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0

[:, :, 2] =
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.148923  0.148923  0.148923  0.148923     0.148923  0.148923  0.148923
 0.400332  0.400332  0.400332  0.400332     0.400332  0.400332  0.400332
 0.697746  0.697746  0.697746  0.697746     0.697746  0.697746  0.697746
 1.01722   1.01722   1.01722   1.01722      1.01722   1.01722   1.01722
 1.34336   1.34336   1.34336   1.34336   …  1.34336   1.34336   1.34336
 1.66501   1.66501   1.66501   1.66501      1.66501   1.66501   1.66501
 1.97352   1.97352   1.97352   1.97352      1.97352   1.97352   1.97352
 2.26207   2.26207   2.26207   2.26207      2.26207   2.26207   2.26207
 2.52509   2.52509   2.52509   2.52509      2.52509   2.52509   2.52509
 ⋮                                       ⋱            ⋮
 2.26207   2.26207   2.26207   2.26207      2.26207   2.26207   2.26207
 1.97352   1.97352   1.97352   1.97352      1.97352   1.97352   1.97352
 1.66501   1.66501   1.66501   1.66501   …  1.66501   1.66501   1.66501
 1.34336   1.34336   1.34336   1.34336      1.34336   1.34336   1.34336
 1.01722   1.01722   1.01722   1.01722      1.01722   1.01722   1.01722
 0.697746  0.697746  0.697746  0.697746     0.697746  0.697746  0.697746
 0.400332  0.400332  0.400332  0.400332     0.400332  0.400332  0.400332
 0.148923  0.148923  0.148923  0.148923  …  0.148923  0.148923  0.148923
 0.0       0.0       0.0       0.0          0.0       0.0       0.0

Choosing the types of Jacobi functions

When using an implicit or semi-implicit solver of differential equations, the Jacobian must be constructed over many iterations, and this can be one of the most expensive steps. To achieve maximum efficiency in solving rigid equations, it is necessary to optimize two parts: the sparsity pattern and the construction of the Jacobian. The construction is filling in the matrix J with values, and the sparsity pattern is which matrix J to use.

The sparsity pattern is set by the prototype matrix jac_prototype, which will be copied for use as J. By default, J is a Matrix, i.e. a dense matrix. However, if the sparsity of the problem is known, you can transfer a different type of matrix. For example, SparseMatrixCSC gives a sparse matrix. Other types of sparse matrices are listed below.

Two-diagonal
Three- diagonal
Symmetrical tridiagonal
Lenta (https://github.com/JuliaLinearAlgebra/BandedMatrices.jl [BandedMatrices.jl])
Block-tape (https://github.com/JuliaLinearAlgebra/BlockBandedMatrices.jl [BlockBandedMatrices.jl])

DifferentialEquations.jl will use this type of matrix internally, which will speed up factorization by using specialized forms.

Declaring a sparse Jacobi matrix using automatic sparsity detection

The sparsity of the Jacobian is declared by the argument jac_prototype in `ODEFunction'. Note that this should only be done if sparsity is high, for example, 0.1% of the matrix is non-zero, otherwise the costs of sparse matrices may be higher than the benefits of sparse differentiation.

One of the useful companion tools for DifferentialEquations.jl is https://github.com/JuliaSymbolics/Symbolics.jl [Symbolics.jl]. This allows you to automatically declare the sparsity types of the Jacobian. To see this in action, you can take du and u and call jacobian_sparsity for a function with example arguments, and you will get a sparse matrix with a template that can be converted to `jac_prototype'.

using Symbolics
du0 = copy(u0)
jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> brusselator_2d_loop(du, u, p, 0.0),
                                           du0, u0)

2048×2048 SparseArrays.SparseMatrixCSC{Bool, Int64} with 12288 stored entries:
⎡⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⎥
⎢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⎥
⎢⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⎥
⎢⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⎦

Now we just pass the sparsity pattern to the 'ODEFunction', as before.:

f = ODEFunction(brusselator_2d_loop; jac_prototype = float.(jac_sparsity))

(::ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.brusselator_2d_loop), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}) (generic function with 1 method)

Creating an ODEProblem:

prob_ode_brusselator_2d_sparse = ODEProblem(f, u0, (0.0, 11.5), p)

ODEProblem with uType Array{Float64, 3} and tType Float64. In-place: true
timespan: (0.0, 11.5)
u0: 32×32×2 Array{Float64, 3}:
[:, :, 1] =
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 ⋮                                  ⋱                      ⋮
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534  …  0.568534  0.326197  0.121344  0.0
 0.0  0.121344  0.326197  0.568534     0.568534  0.326197  0.121344  0.0

[:, :, 2] =
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.148923  0.148923  0.148923  0.148923     0.148923  0.148923  0.148923
 0.400332  0.400332  0.400332  0.400332     0.400332  0.400332  0.400332
 0.697746  0.697746  0.697746  0.697746     0.697746  0.697746  0.697746
 1.01722   1.01722   1.01722   1.01722      1.01722   1.01722   1.01722
 1.34336   1.34336   1.34336   1.34336   …  1.34336   1.34336   1.34336
 1.66501   1.66501   1.66501   1.66501      1.66501   1.66501   1.66501
 1.97352   1.97352   1.97352   1.97352      1.97352   1.97352   1.97352
 2.26207   2.26207   2.26207   2.26207      2.26207   2.26207   2.26207
 2.52509   2.52509   2.52509   2.52509      2.52509   2.52509   2.52509
 ⋮                                       ⋱            ⋮
 2.26207   2.26207   2.26207   2.26207      2.26207   2.26207   2.26207
 1.97352   1.97352   1.97352   1.97352      1.97352   1.97352   1.97352
 1.66501   1.66501   1.66501   1.66501   …  1.66501   1.66501   1.66501
 1.34336   1.34336   1.34336   1.34336      1.34336   1.34336   1.34336
 1.01722   1.01722   1.01722   1.01722      1.01722   1.01722   1.01722
 0.697746  0.697746  0.697746  0.697746     0.697746  0.697746  0.697746
 0.400332  0.400332  0.400332  0.400332     0.400332  0.400332  0.400332
 0.148923  0.148923  0.148923  0.148923  …  0.148923  0.148923  0.148923
 0.0       0.0       0.0       0.0          0.0       0.0       0.0

Now let’s see how the sparse version differs from the non-sparse version.:

using BenchmarkTools # для @btime
@btime solve(prob_ode_brusselator_2d, TRBDF2(), save_everystep = false);
@btime solve(prob_ode_brusselator_2d_sparse, TRBDF2(), save_everystep = false);
@btime solve(prob_ode_brusselator_2d_sparse, KenCarp47(linsolve = KLUFactorization()),
             save_everystep = false);

  26.847 s (4836 allocations: 97.99 MiB)
  826.241 ms (36023 allocations: 242.45 MiB)
  755.176 ms (35525 allocations: 14.67 MiB)

Note that, depending on the properties of the sparsity pattern, alternative linear solvers can be tried, such as TRBDF2(linsolve = KLUFactorization()) or `TRBDF2(linsolve = UMFPACKFactorization())'.

Using the Newton-Krylov method without the Jacobian

A completely different way to optimize linear solvers for large sparse matrices is to use the Krylov subspace method. To switch to the Krylov method, it is necessary to choose a linear solver. To change the linear solver, use the linsolve command and select the GMRES linear solver.

@btime solve(prob_ode_brusselator_2d, KenCarp47(linsolve = KrylovJL_GMRES()),
             save_everystep = false);

  1.275 s (324459 allocations: 42.38 MiB)

This acceleration does not require defining a sparsity pattern, and therefore it can be an easier way to scale when solving large tasks. For more information about choosing a linear solver, see documentation on the linear solver. The variants of linsolve are any valid solvers. https://linearsolve.sciml.ai/dev /[LinearSolve.jl].

When switching to the Krylov linear solver, the ODE solver automatically switches to the Jacobian-free mode, which significantly reduces the amount of memory required. This behavior can be redefined by adding concrete_jac=true to the algorithm.

Adding a precondition

Any one can be used as a precondition in the interface of the linear solver. https://docs.sciml.ai/LinearSolve/stable/basics/Preconditioners /[precondition compatible with LinearSolve.jl]. To define preconditioners in compatible solvers of rigid ODES, it is necessary to define the function precs, which returns left and right preconditioners, matrices approximating the inverse W = I - gamma*J used in solving ODES. An example of this using https://github.com/haampie/IncompleteLU.jl [IncompleteLU.jl] looks like this:

using IncompleteLU
function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
    else
        Pl = Plprev
    end
    Pl, nothing
end

# Требуется из-за ошибки в Krylov.jl: https://github.com/JuliaSmoothOptimizers/Krylov.jl/pull/477.
Base.eltype(::IncompleteLU.ILUFactorization{Tv, Ti}) where {Tv, Ti} = Tv

@btime solve(prob_ode_brusselator_2d_sparse,
             KenCarp47(linsolve = KrylovJL_GMRES(), precs = incompletelu,
                       concrete_jac = true), save_everystep = false);

  431.888 ms (72073 allocations: 59.19 MiB)

Pay attention to several features of this precondition. This precondition uses a sparse Jacobian, and therefore we set concrete_jac=true to tell the algorithm to generate a Jacobian (otherwise, GMRES uses a non-Jacobian algorithm by default). Then newW = true each time a new matrix W is calculated, and newW=nothing at the solver startup stage. Thus, we check newW === nothing || newW and, if the expression is true, then we update the precondition only at these points, otherwise we simply switch to the previous version. Using convert(AbstractMatrix,W), we get a specific matrix W (matching jac_prototype, and therefore SpraseMatrixCSC), which can be used in defining a precondition. Then we use IncompleteLU.ilu in this sparse matrix to generate a preconditioner. We return Pl,nothing, which means that the precondition is a left precondition and there is no right precondition.

Thus, this method uses both the Krylov solver and the sparse Jacobian. Moreover, it is faster than both implementations. IncompleteLU is unstable in that it requires a well-configured parameter t. Another option is to use a package https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl [AlgebraicMultigrid.jl], which is more automatic. The setup is very similar to the one described earlier.:

using AlgebraicMultigrid
function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W)))
    else
        Pl = Plprev
    end
    Pl, nothing
end

@btime solve(prob_ode_brusselator_2d_sparse,
             KenCarp47(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid,
                       concrete_jac = true), save_everystep = false);

  684.716 ms (60155 allocations: 124.73 MiB)

or using the Jacobi smoothing:

function algebraicmultigrid2(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        A = convert(AbstractMatrix, W)
        Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(A,
                                                                                presmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
                                                                                                                                  1))),
                                                                                postsmoother = AlgebraicMultigrid.Jacobi(rand(size(A,
                                                                                                                                   1)))))
    else
        Pl = Plprev
    end
    Pl, nothing
end

@btime solve(prob_ode_brusselator_2d_sparse,
             KenCarp47(linsolve = KrylovJL_GMRES(), precs = algebraicmultigrid2,
                       concrete_jac = true), save_everystep = false);

  561.559 ms (67478 allocations: 126.46 MiB)

For more information about the preconditioner interface, see documentation on the linear solver.

Although most of the settings allow you to automate the transition to using Sundials, there are some differences between pure Julia implementations and Sundials implementations that you should pay attention to. They are described in detail in https://docs.sciml.ai/DiffEqDocs/stable/api/sundials /[documentation on the Sundials solver], but here we will highlight the main details that need to be kept in mind.

The definition of a sparse matrix and Jacobian for Sundials is the same as for any other package. The main difference lies in the choice of a linear solver. In Sundials` the linear solver is selected using the symbol in linear_solver' from a list of preset values. Of particular note is the choice of `:Band for the tape matrix and :GMRES for using GMRES. If you use Sundials, :GMRES does not require the definition of a JacVecOperator and will always use the Newton-Krylov method without Jacobian (with numerical differentiation). Thus, the following can be done for this task.

using Sundials
@btime solve(prob_ode_brusselator_2d, CVODE_BDF(), save_everystep = false);
# Самое простое ускорение: использовать :LapackDense
@btime solve(prob_ode_brusselator_2d, CVODE_BDF(linear_solver = :LapackDense),
             save_everystep = false);
# Версия GMRES: не требует ничего дополнительного!
@btime solve(prob_ode_brusselator_2d, CVODE_BDF(linear_solver = :GMRES),
             save_everystep = false);

  38.983 s (34266 allocations: 1.65 MiB)
  14.505 s (34266 allocations: 1.65 MiB)
  361.438 ms (34346 allocations: 1.60 MiB)

Note that using sparse matrices in Sundials requires a Jacobian analytic function. For automatic generation, we will use the `modelingtoolkitize' package. https://mtk.sciml.ai/dev /[ModelingToolkit.jl]:

using ModelingToolkit
prob_ode_brusselator_2d_mtk = ODEProblem(modelingtoolkitize(prob_ode_brusselator_2d_sparse),
                                         [], (0.0, 11.5), jac = true, sparse = true);
# @btime solve(prob_ode_brusselator_2d_mtk,CVODE_BDF(linear_solver=:KLU),save_everystep=false); # очень медленно компилируется

Using Preconditioners with SUNDIALS

Information about configuring the preconditioner using Sundials can be found on the Sundials solver page. The Sundials algorithms differ significantly from the standard Julia-based algorithms in that all operations with the Jacobian matrix are performed by the user. To do this, define the function psetup, which sets the precondition, and then the function prec, which is the action of the precondition with a vector. For the 'psetup` function, you must first calculate the matrix W = I - gamma*J, and then calculate the precondition for it. For the above ILU example, these actions are performed for Sundials as follows:

using LinearAlgebra
u0 = prob_ode_brusselator_2d_mtk.u0
p = prob_ode_brusselator_2d_mtk.p
const jaccache = prob_ode_brusselator_2d_mtk.f.jac(u0, p, 0.0)
const W = I - 1.0 * jaccache

prectmp = ilu(W, τ = 50.0)
const preccache = Ref(prectmp)

function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
    if jok
        prob_ode_brusselator_2d_mtk.f.jac(jaccache, u, p, t)
        jcurPtr[] = true

        # W = I - gamma*J
        @. W = -gamma * jaccache
        idxs = diagind(W)
        @. @view(W[idxs]) = @view(W[idxs]) + 1

        # Создает предобусловливатель для W
        preccache[] = ilu(W, τ = 5.0)
    end
end

Then the precondition action is reduced to simply using the ldiv! generated precondition.:

function precilu(z, r, p, t, y, fy, gamma, delta, lr)
    ldiv!(z, preccache[], r)
end

After that, we simply pass these functions to the Sundials solver by selecting prec_side=1 to indicate that it is a left-handed preconditioner.:

@btime solve(prob_ode_brusselator_2d_sparse,
             CVODE_BDF(linear_solver = :GMRES, prec = precilu, psetup = psetupilu,
                       prec_side = 1), save_everystep = false);

And similarly for an algebraic multiple grid:

prectmp2 = aspreconditioner(ruge_stuben(W,
                                        presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
                                                                                          1))),
                                        postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
                                                                                           1)))))
const preccache2 = Ref(prectmp2)
function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
    if jok
        prob_ode_brusselator_2d_mtk.f.jac(jaccache, u, p, t)
        jcurPtr[] = true

        # W = I - gamma*J
        @. W = -gamma * jaccache
        idxs = diagind(W)
        @. @view(W[idxs]) = @view(W[idxs]) + 1

        # Создает предобусловливатель для W
        preccache2[] = aspreconditioner(ruge_stuben(W,
                                                    presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
                                                                                                      1))),
                                                    postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
                                                                                                       1)))))
    end
end

function precamg(z, r, p, t, y, fy, gamma, delta, lr)
    ldiv!(z, preccache2[], r)
end

@btime solve(prob_ode_brusselator_2d_sparse,
             CVODE_BDF(linear_solver = :GMRES, prec = precamg, psetup = psetupamg,
                       prec_side = 1), save_everystep = false);