# Error when using sparse Jacobians/Krylov solvers in RadauIIA methods

**URL:** <https://discourse.julialang.org/t/error-when-using-sparse-jacobians-krylov-solvers-in-radauiia-methods/94629>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [February 14, 2023, 3:49pm UTC](https://discourse.julialang.org/t/error-when-using-sparse-jacobians-krylov-solvers-in-radauiia-methods/94629 "2023-02-14T15:49:52Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![wpowell96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wpowell96/32/46900_2.png) [@wpowell96](https://discourse.julialang.org/u/wpowell96)\
**Post date:** [February 14, 2023, 3:49pm UTC](https://discourse.julialang.org/t/error-when-using-sparse-jacobians-krylov-solvers-in-radauiia-methods/94629/1 "2023-02-14T15:49:52Z")

</div>

Yesterday I posted an [issue](https://github.com/SciML/DifferentialEquations.jl/issues/934) but realized that I forgot to crosspost this here for visibility.

I was trying to generate some work-precision diagrams for MOL PDEs based on this [old SciML benchmark](https://github.com/SciML/SciMLBenchmarks.jl/blob/master/benchmarks/MOLPDE/ks_fdm_wpd.jmd). The syntax for declaring linear solvers has changed, but I believe this is unrelated. The same error also occurs on the updated version of this benchmark [here](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/MOLPDE/ks_spectral_wpd/). These benchmarks generate a `SplitODEProblem` for the Kuramoto-Sivashinsky equation and solve it using a variety of methods.

To generate a test solution I tried using `RadauIIA5` with very low tolerances and the following error is thrown:

```julia
sol = solve(prob, RadauIIA5(autodiff=false); abstol=1e-5, reltol=1e-5)
┌ Warning: Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/P7HJO/src/alg_utils.jl:241
ERROR: MethodError: no method matching similar(::SparseDiffTools.JacVec{SciMLBase.UJacobianWrapper{SplitFunction{true, SciMLBase.FullSpecialize, ODEFunction{true, SciMLBase.FullSpecialize, DiffEqArrayOperator{Float64, Matrix{Float64}, typeof(SciMLBase.DEFAULT_UPDATE_FUNC)}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, var"#13#14"{Vector{Float64}, Matrix{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, UniformScaling{Bool}, Vector{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, Vector{Float64}, Vector{Float64}, Vector{Float64}}, ::Type{ComplexF64})

```

Other implicit solvers worked fine. I figured this was an issue with giving it a `SplitODEProblem` with some sparse operator, so I tried this variant of the Brusselator system from the docs with a standard RHS and tried to get it to use a Krylov solver and the same error is thrown. `ImplicitEuler` and `KenCarp4` work fine in this example, and every other implict method I tried works as well except for `RadauIIA3`, which throws the same error. Code is below.

```julia
using DifferentialEquations, LinearAlgebra, SparseArrays, Krylov

const N = 10
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)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,u0,(0.,0.1),p)

sol = solve(prob_ode_brusselator_2d, ImplicitEuler(autodiff=false, linsolve=KrylovJL_GMRES()))

sol = solve(prob_ode_brusselator_2d, KenCarp4(autodiff=false, linsolve=KrylovJL_GMRES()))

sol = solve(prob_ode_brusselator_2d, RadauIIA5(autodiff=false, linsolve=KrylovJL_GMRES()))

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 14, 2023, 4:13pm UTC](https://discourse.julialang.org/t/error-when-using-sparse-jacobians-krylov-solvers-in-radauiia-methods/94629/2 "2023-02-14T16:13:23Z")

</div>

Open an issue please.

---

<div class="post-metadata">

**Author:** ![wpowell96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wpowell96/32/46900_2.png) [@wpowell96](https://discourse.julialang.org/u/wpowell96)\
**Post date:** [February 14, 2023, 4:20pm UTC](https://discourse.julialang.org/t/error-when-using-sparse-jacobians-krylov-solvers-in-radauiia-methods/94629/3 "2023-02-14T16:20:24Z")

</div>

Issue is linked in the post. Here it is again for visibility: [https://github.com/SciML/DifferentialEquations.jl/issues/934](https://github.com/SciML/DifferentialEquations.jl/issues/934)