Yesterday I posted an issue 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. 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. 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:
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.
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()))