requirement on Jacobian for CVODE_BDF

[dchang24]

When I use CVODE_BDF to solve a coupled stiff ODE (a PDE discretized using method of lines) with around ~1000 variables, is it necessary to provide the exact Jacobian or is it sufficient to just provide the sparsity structure (it is sparse) of the Jacobian? In Matlab using ode15s or in scipy using scipy.integrate.solve_ivp we can just provide the sparsity structure of the Jacobian and the solver will compute the Jacobian numerically using finite difference. This is what I have been doing since I know the sparsity structure (which already gives a pretty decent speedup compared to not giving the sparsity structure) but computing the Jacobian is fairly complicated. What about in Julia DifferentialEquations? It's not clear from the documentation what's actually required. I have something like

prob = ODEProblem(dhdt!, h0, (0.0,tf), p)
sol = solve(prob, alg=CVODE_BDF())

which gives what I expect by comparing with solutions computed with MATLAB. But when I try something like

f = ODEFunction(dhdt!; jac_prototype=jpattern) #jpattern is a sparse matrix which gives the sparsity structure of Jacobian
prob_sparse = ODEProblem(f, h0, (0.0,tf), p)
sol_sparse = solve(prob_sparse, alg=CVODE_BDF(linear_solver=:BCG))

it just gives the wrong answer.

[ChrisRackauckas]

For everything but CVODE we take over and use SparseDiffTools.jl to generate it directly from the sparsity pattern. With CVODE we leave the default Jacobian handling to Sundials, but in 2021 we probably shouldn't because at this point we want to exploit all of these sparsity and AD tools which are missing from Sundials but exist in the SciML ecosystem. So we should really slap our standard DiffEq interface on this portion and take control away from Sundials for the Jacobian building.

@jd-lara might be able to help here.

Moving to Sundials.jl because it's solver-specific.

[jd-lara]

I am not sure what could this be. There are tests in pure Sundials that do work properly, we need to check that the common interface is sending to Sundials the correct thing.

[jd-lara]

@yhchang96 try without specifying the linear solver. I've had issues with the default tolerances of the linear solvers in Sundials before.