So to start, there are two other methods I thought of.

One is to build the preconditioner based on an approximation. If you do sparsity detection without the integral, you’ll get a sparse matrix that approximates your Jacobian. You can then use that inside of DifferentialEquations.jl as the `jac_prototype`

and that will generate an approximate Jacobian through color differentiation to that sparse matrix. That approximate Jacobian can be used inside of an ILU/AMG preconditioner. This is fine because preconditioners are only supposed to be an approximation anyways: whether this works is really whether it’s a “good enough” approximation, which depends on your equation and the parameter values. This only takes like 5 lines of code to try (copy-pasted out of the tutorial) so I would try this first.

The second thing you can do is use an IMEX integrator:

https://diffeq.sciml.ai/stable/types/split_ode_types/

https://diffeq.sciml.ai/stable/solvers/split_ode_solve/#OrdinaryDiffEq.jl

I’d recommend KenCarp47, where you split the integral to be part of the explicit part. If the integral has well-behaved eigenvalues, this would let you treat it explicitly, and then your true Jacobian is all of the non-integral stuff, and then once again sparse matrix preconditioners etc. etc. that tutorial applies.

One of those two “should” work. But let’s address the other thing just for full clarity:

Note that those preconditioners are not matrix-free. ILU needs to know the sparse matrix and a `tau`

cutoff value in order to construct the preconditioner. Algebraic Multigrids also need to know the structure of the equations through the sparse matrix. That’s really the issue: preconditioners need to match details of the equation since they need to be an approximate inverse, so the easiest way to construct an algorithm that knows “some structure of the equations” is to have as input a sparse matrix because that naturally will have much of the structure for many PDEs.

That said, as @stevengj discussed, there are preconditioners that do not require you have the sparse matrix. All of those will naturally require that you write a function based solely on your PDE, for example fast multiple methods and geometric multigrids, so it’s hard to ever give someone a software that does this automatically like you can with AMG and ILU (well, there is a way to do it, but that’ll be a research project over the next year ).

I think the easiest to write might be a geometric multigrid, since it just requires knowing how to write your PDE with a different dx, along with expanding and contracting the grid. Gil has some good notes on it:

But other than that, for now it’s kind of “may the force be with you” in terms of automation, which I hope can be solved better in the future.