Benchmark to investigate for Newton-Krylov performance

[ChrisRackauckas]

using OrdinaryDiffEq, Sundials, DiffEqOperators, BenchmarkTools

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.
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., 10., 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)

Jv = JacVecOperator(brusselator_2d_loop,u0,p,0.0)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     u0,(0.,11.5),p)

f = ODEFunction(brusselator_2d_loop;jac_prototype=Jv)
prob_ode_brusselator_2d_jacfree = ODEProblem(f,u0,(0.,11.5),p)
@btime solve(prob_ode_brusselator_2d_jacfree,TRBDF2(linsolve=LinSolveGMRES()),save_everystep=false)
@btime solve(prob_ode_brusselator_2d,CVODE_BDF(linear_solver=:GMRES),save_everystep=false)

@YingboMa

[ChrisRackauckas]

10s vs 0.5s

[rveltz]

I have another issue pending to be submitted to DiffEqBase which could explain this. Somehow, I got the feeling that the tolerance tol is not well handled in function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false; Pl=nothing, Pr=nothing, tol=nothing, kwargs...). You can see that if you enable the log in IterativeSolvers, hence https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/366

[vpuri3]

Yes, i printed the tolerance for a small case and it is weirdly high:

example:

"""
 OrdinaryDiffEq/test/interface/linear_nonlinear_tests.jl
"""
n = 10
A = Tridiagonal(-ones(n-1),2ones(n),-ones(n-1))
rn = (du, u, p, t) -> begin
    mul!(du, A, u)
end
u0 = rand(n)
prob = ODEProblem(ODEFunction(rn, jac_prototype=JacVecOperator{Float64}(rn, u0; autodiff=false)), u0, (0, 10.))
sol = solve(prob, QNDF(autodiff=false, linsolve=LinSolveGMRES()));

with DiffEqBase/linear_nonlinear.jl modified as follows:

function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false; Pl=nothing, Pr=nothing, reltol=eps(eltype(x)), kwargs...)
  if f.iterable === nothing
    Pl = ComposePreconditioner(f.Pl, Pl, true)
    Pr = ComposePreconditioner(f.Pr, Pr, false)
    
    reltol = checkreltol(reltol)
    f.iterable = f.generate_iterator(x,A,b,f.args...;
                                     initially_zero=true,restart=5,
                                     maxiter=20,abstol=1e-12,reltol=reltol,
                                     Pl=Pl,Pr=Pr,
                                     kwargs...,f.kwargs...)

  end
  x .= false
  iter = f.iterable
  purge_history!(iter, x, b)

  for residual in iter
      println("GMRES step $(iter.k), residual: $(iter.residual.current), tol=$(iter.tol)")
  end
  return nothing
end

julia> include("krylov_test.jl")
GMRES step 1, residual: 0.004468827765402508, tol=5.481711392696291
GMRES step 1, residual: 0.0006976753003995138, tol=5.481711392696291
GMRES step 1, residual: 0.00010892502571461855, tol=5.481711392696291
GMRES step 1, residual: 0.0067410564826430016, tol=5.481711392696291
GMRES step 1, residual: 0.0031643581369913687, tol=5.481711392696291
GMRES step 1, residual: 0.0009081523038929753, tol=5.481711392696291
GMRES step 1, residual: 0.00013766281617396835, tol=5.481711392696291
GMRES step 1, residual: 0.0034366200233135642, tol=5.481711392696291
GMRES step 1, residual: 0.0018549989980046486, tol=5.481711392696291
GMRES step 1, residual: 0.0005839668067025519, tol=5.481711392696291
GMRES step 1, residual: 0.005929950997054515, tol=5.481711392696291
GMRES step 1, residual: 0.16739523946435644, tol=5.481711392696291 
GMRES step 1, residual: 0.09269193011431896, tol=5.481711392696291 
GMRES step 1, residual: 0.0733150942936091, tol=5.481711392696291  
GMRES step 1, residual: 0.07245774505656238, tol=5.481711392696291 
GMRES step 1, residual: 0.3740037683520435, tol=5.481711392696291  
GMRES step 1, residual: 1.571981217219508, tol=5.481711392696291   

[rveltz]

Why iter.tol is 5.4817113??

Otherwise, the behavior is actually normal. The matrix is badly conditioned. See:

using BifurcationKit
n = 10
A = Tridiagonal(-ones(n-1),2ones(n),-ones(n-1)) |> Array
rhs = rand(n)
ls = GMRESIterativeSolvers(verbose = true, log = true, N = n)
sol = ls(A, rhs)

gives

=== gmres ===
rest    iter    resnorm
  1   1 1.65e+00
  1   2 1.51e+00
  1   3 1.21e+00
  1   4 9.02e-01
  1   5 5.33e-01
  1   6 3.43e-01
  1   7 1.66e-01
  1   8 1.11e-01
  1   9 4.21e-02
  1  10 5.86e-15

Better example for iterative methods:

n = 100
A = I + sprand(n,n,1/n)
rhs = rand(n)
ls = GMRESIterativeSolvers(verbose = true, log = true, N = n)
sol = ls(A, rhs)

gives

=== gmres ===
rest    iter    resnorm
  1   1 2.34e+00
  1   2 1.34e+00
  1   3 5.37e-01
  1   4 1.56e-01
  1   5 5.98e-02
  1   6 1.88e-02
  1   7 8.50e-03
  1   8 1.80e-03
  1   9 7.18e-04
  1  10 1.51e-04
  1  11 4.92e-05
  1  12 4.90e-06
  1  13 2.41e-06
  1  14 4.88e-07
  1  15 4.78e-08

I apologize to use BifurcationKit for the linear solver but I am very accustomed to the interface which is unified for several packages. It should be similar to LinSolveIterativeSolvers except I do not use the iterator form of gmres