Ok, this thread kind of went on in my head and I felt that one possibility was still missing…
I did some experiments with preventing setting zero entries in jacobian.jl . As a result, we get:
So this gives an improvement of assembly times by a constant factor, but we still see the O(N^2) behavior. This appears to be due to the fact that without further a priori information, ForwardDiff assumes that the jacobian is a full matrix and calculates all these zeros, still leading to O(N^2) complexity.
For a known sparsity pattern (as we have for this case), the remedy suggested is to provide the sparsity pattern along with coloring information such that derivative calculations can be scheduled to handle only the nonzero entry cases.
So instead of
ForwardDiff.jacobian!(jac, ffd, X, Y);
for jac with unknown sparsity pattern we do for jac with known sparsity pattern and nonzero entries:
using SparseDiffTools
...
colors = matrix_colors(jac)
forwarddiff_color_jacobian!(jac, ffd, X, colorvec = colors)
(for unpatched ForwardDiff) resulting in:
As we see, we now get O(N) behavior because now known zero jacobian entries are not calculated to begin with.
So we have two ways to handle the sparsity here: matrix coloring (once the sparsity pattern is known) or atomic assembly (which with ExtendableSparse performs reasonably well without a priori information about the sparsity pattern).