Документация Engee

Решение больших жестких уравнений

В этом руководстве рассматриваются дополнительные возможности эффективного решения больших жестких обыкновенных дифференциальных уравнений. Для решения жестких обыкновенных дифференциальных уравнений требуется специализировать линейный решатель в свойствах якобиана, чтобы сократить количество линейных решений и обратных решений . Эти же функции и средства управления распространяются и на жесткие уравнения SDE, DDE, DAE и т. д. Данное руководство предназначено для крупномасштабных моделей, например полученных для полудискретизации дифференциальных уравнений в частных производных (PDE). Например, мы будем использовать жесткое дифференциальное уравнение брюсселятора (BRUSS).

Это руководство предназначено для опытных пользователей, желающих ознакомиться с расширенными возможностями. DifferentialEquations.jl автоматизирует большую часть использования, поэтому пользователям рекомендуется сначала попробовать solve(prob) с автоматическим алгоритмом.

Определение уравнения брюсселятора

Этот раздел можно пропустить: в нем просто сформулирована задача примера.

Уравнение PDE брюсселятора определяется следующим образом:

где

и начальными условиями являются

с периодическим граничным условием

во временном диапазоне ].

Для решения этого уравнения PDE дискретизируем его в систему уравнений ODE методом конечных разностей. Дискретизируем u и v в массивы значений в каждой временной точке: u[i,j] = u(i*dx,j*dy) для некоторого выбора dx/dy, и то же самое для v. Тогда уравнение ODE определяется с помощью U[i,j,k] = [u v]. Оператор второй производной, лапласиан, дискретизируется и становится трехдиагональной матрицей с [1 -2 1] и 1 в правом верхнем и левом нижнем углах. Затем нелинейные функции применяются в каждой точке пространства (они транслируются). Используйте dx=dy=1/32.

Результирующее определение ODEProblem имеет следующий вид:

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

Выбор типов функций Якоби

При использовании неявного или полуявного решателя дифференциальных уравнений якобиан должен быть построен на многих итерациях, и это может быть одним из самых дорогостоящих этапов. Для достижения максимальной эффективности при решении жестких уравнений необходимо оптимизировать две части: шаблон разреженности и построение якобиана. Конструкция — это заполнение матрицы J значениями, а шаблон разреженности — это то, какую матрицу J использовать.

Шаблон разреженности задается матрицей-прототипом jac_prototype, которая будет скопирована для использования в качестве J. По умолчанию J — это Matrix, т. е. плотная матрица. Однако если известна разреженность задачи, можно передать матрицу другого типа. Например, SparseMatrixCSC дает разреженную матрицу. Ниже перечислены другие типы разреженных матриц.

  • Двухдиагональная

  • Трехдиагональная

  • Симметричная трехдиагональная

  • Ленточная (BandedMatrices.jl)

  • Блочно-ленточная (BlockBandedMatrices.jl)

DifferentialEquations.jl будет внутренним образом использовать этот тип матрицы, что позволит ускорить факторизацию за счет применения специализированных форм.

Объявление разреженной матрицы Якоби с помощью автоматического определения разряженности

Разреженность якобиана объявляется аргументом jac_prototype в ODEFunction. Обратите внимание, что это следует делать только в том случае, если разреженность велика, например, 0,1 % матрицы является ненулевым, иначе затраты на разреженные матрицы могут оказаться выше, чем выгода от разреженного дифференцирования.

Одним из полезных инструментов-компаньонов для DifferentialEquations.jl является Symbolics.jl. Это позволяет автоматически объявлять типы разреженности якобиана. Чтобы увидеть это в действии, можно взять du и u и вызвать jacobian_sparsity для функции с примерами аргумента, и будет выдана разреженная матрица с шаблоном, которую можно преобразовать в 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:
⎡⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⠀⎥
⎢⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⎥
⎢⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⎥
⎢⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⎦

Обратите внимание, что Julia прекрасно выводит шаблон разреженности. Это здорово, но создание вручную было бы утомительным. Теперь мы просто передаем его в ODEFunction, как и раньше:

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)

Создадим 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

Теперь посмотрим, как версия с разреженностью отличается от версии без разреженности:

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)

Заметим, что в зависимости от свойств шаблона разреженности можно попробовать альтернативные линейные решатели, такие как TRBDF2(linsolve = KLUFactorization()) или TRBDF2(linsolve = UMFPACKFactorization()).

Использование метода Ньютона — Крылова без якобиана

Совершенно иной способ оптимизации линейных решателей для больших разреженных матриц заключается в использовании метода подпространств Крылова. Для перехода к методу Крылова требуется выбрать линейный решатель. Для смены линейного решателя используем команду linsolve и выбираем линейный решатель GMRES.

@btime solve(prob_ode_brusselator_2d, KenCarp47(linsolve = KrylovJL_GMRES()),
             save_everystep = false);
  1.275 s (324459 allocations: 42.38 MiB)

Для такого ускорения не требуется определение шаблона разреженности, и, следовательно, оно может быть более простым способом масштабирования при решении больших задач. Более подробная информация о выборе линейного решателя приведена в документации по линейному решателю. Вариантами linsolve являются любые допустимые решатели LinearSolve.jl.

При переходе на линейный решатель Крылова решатель ODE автоматически переводится в режим без якобиана, что значительно сокращает объем требуемой памяти. Это поведение можно переопределить, добавив concrete_jac=true в алгоритм.

Добавление предобусловливателя

В качестве предобусловливателя в интерфейсе линейного решателя может быть использован любой предобусловливатель, совместимый с LinearSolve.jl. Для определения предобусловливателей в совместимых решателях жестких ODE необходимо определить функцию precs, которая возвращает левый и правый предобусловливатели, матрицы, аппроксимирующие обратные W = I - gamma*J, используемые при решении ODE. Пример этого с использованием IncompleteLU.jl выглядит следующим образом:

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)

Обратите внимание на несколько особенностей этого предобусловливателя. Этот предобусловливатель использует разреженный якобиан, и поэтому мы задаем concrete_jac=true, чтобы указать алгоритму генерировать якобиан (в противном случае по умолчанию в GMRES используется алгоритм без якобиана). Затем newW = true каждый раз, когда вычисляется новая матрица W, и newW=nothing на этапе запуска решателя. Таким образом, мы проверяем newW === nothing || newW и, если выражение верно, то только в этих точках обновляем предобусловливатель, в противном случае мы просто переходим к предыдущей версии. Используя convert(AbstractMatrix,W), мы получаем конкретную матрицу W (совпадающую с jac_prototype, а следовательно, SpraseMatrixCSC), которую можно использовать в определении предобусловливателя. Затем мы используем IncompleteLU.ilu в этой разреженной матрице для генерации предобусловливателя. Мы возвращаем Pl,nothing, что означает, что предобусловливатель является левым предобусловливателем и правого предобусловливания не существует.

Таким образом, в данном методе используется как решатель Крылова, так и разреженный якобиан. Мало того, он быстрее обеих реализаций. IncompleteLU нестабилен в том, что требует хорошо настроенного параметра τ. Другой вариант — использовать пакет AlgebraicMultigrid.jl, который является более автоматическим. Настройка очень похожа на ту, которая описана ранее:

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)

или с помощью сглаживателя Якоби:

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)

Более подробная информация об интерфейсе предобусловливателя приведена в документации по линейному решателю.

Обработка, связанная с SUNDIALS

Несмотря на то, что большая часть настроек позволяет автоматизировать переход к использованию Sundials, существуют некоторые различия между чистыми реализациями Julia и реализациями Sundials, на которые следует обратить внимание. Они подробно описаны в документации по решателю Sundials, но здесь мы выделим основные детали, которые нужно иметь в виду.

Определение разреженной матрицы и якобиана для Sundials происходит так же, как и для любого другого пакета. Основное отличие заключается в выборе линейного решателя. В Sundials выбор линейного решателя осуществляется с помощью символа в linear_solver из списка предустановленных значений. Особо следует отметить выбор :Band для ленточной матрицы и :GMRES для использования GMRES. Если вы используете Sundials, :GMRES не требует определения JacVecOperator и всегда будет использовать метод Ньютона-Крылова без якобиана (с численным дифференцированием). Таким образом, по данной задаче можно сделать следующее.

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)

Обратите внимание, что для использования разреженных матриц в Sundials требуется аналитическая функция якобиана. Для автоматической генерации мы будем использовать modelingtoolkitize пакета 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); # очень медленно компилируется

Использование предобусловливателей с SUNDIALS

Сведения о настройке предобусловливателя с помощью Sundials можно найти на странице решателей Sundials. Алгоритмы Sundials существенно отличаются от стандартных алгоритмов на основе Julia тем, что все операции с матрицей Якобиана выполняет пользователь. Для этого необходимо определить функцию psetup, задающую предобусловливатель, а затем функцию prec, являющуюся действием предобусловливателя c вектором. Для функции psetup необходимо сначала вычислить матрицу W = I - gamma*J, а затем вычислить для нее предобусловливатель. Для приведенного выше примера ILU эти действия выполняются для Sundials следующим образом:

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

Затем действие предобусловливателя сводится к простому использованию ldiv! сгенерированного предобусловливателя:

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

После этого мы просто передаем эти функции решателю Sundials, выбрав prec_side=1, чтобы указать, что он является левым предобусловливателем:

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

И аналогично для алгебраической множественной сетки:

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);