Thanks Chris - I have tried to do it but I got a bit stuck.
It seems that to initialize the ilu one should know the coefficient matrix A, but how do we get that in our problem? Can it be extracted somehow from the output of the discretize command?
prob = discretize(pdesys, discretization)
using IncompleteLU, LinearAlgebra
LU = ilu(A, τ = 0.1)
sol = solve(prob, Rodas5P(linsolve = IterativeSolversJL_GMRES(), precs=LU));