I’m trying to solve a DAE (mass matrix form) for a 2D spectral discretization. The number of DOF is going to be large, and the Jacobian, which I compute explicitly, is not especially sparse.
I’d like to use a preconditioned Krylov solver operating in a matrix-free mode, as applying the Jacobian can be done much faster than the equivalent matrix-vector multiplication. But I can’t figure out how to do that within DifferentialEquations.jl.
Going by the docs, the linsolve option to the native Julia stiff solvers apparently expects to pass matrices around. In SUNDIALS, which does support “ATimes” functions instead of matrices for its iterative solvers, I can’t figure out how to set them up.
Does anyone have a working example of using matrix-free methods in this context?
2 Likes
It’s not ready yet. I have a plan for how to do it:
https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/220
but it needs to get implemented. Until then, Sundials.jl exposes the whole API so you can define the ATimes function just fine if you pass a cfunction for it. So the direct interface should work, just the high level interface isn’t there yet. I wanna get it done this summer though…
1 Like
you can use the python package petsc4py
Or just PETSc.jl. But that doesn’t actually give us something new here. PETSc is actually using Sundials for its BDF implementations, so you might as well just call Sundials directly if that’s what you want since PETSc isn’t much of a higher level interface anyways.
1 Like
Didn’t think to search the repo issues. Thanks.