InterpolateBlock doesn't stash interpolators

[ReubenHill]

The adjoint InterpolateBlock creates a new interpolator every time it is evaluated, e.g. for calculating the hessian.

Extra interpolators should probably be stashed on the Interpolator as necessary - though this is perhaps not the only solution. I believe @nbouziani has worked on a similar fix elsewhere - comments welcome.

[nbouziani]

Sorry for the delay! I have commented what I think is the same issue, cf. #1981.

[ReubenHill]

(Closed the duplicate issue)

[nbouziani]

I have not dealt with similar issues in #1674 since the interpolate is part of the external operator evaluation process (happens inside the external operator), that is there is no additional InterpolateBlock induced.

Just to be sure I understand correctly: You want to stash Interpolator(expr, space) right ? Is it not just when you interpolate via Interpolator(expr, space).interpolate(...) that there is a computational cost ? If yes, then I don't understand the problem with defining a new Interpolator every time something needs interpolating.

[ReubenHill]

If expr and space are the same multiple times (as I think typically they will be when doing hessian evaluations as part of a firedrake_adjoint.minimize(Ĵ, method='Newton-CG')) you really ought to stash and reuse the Interpolator to save compile time costs

[nbouziani]

Ok I see. If this is just for adjoint purposes, then the information about which interpolator has already been made could also be stashed directly to the InterpolateBlock class. Although I am not sure which option is the most judicious.

[dham]

That doesn't quite work, because you get a new InterpolateBlock for every interpolate operation.

[nbouziani]

I meant having a kind of cache mechanism via shared attributes as it is done for forms with TSFCKernel.

[dham]

Global caches for interpolation matrices are going to get very expensive very fast.

[nbouziani]

Ok I see.

[wence-]

But kernel caching should be implemented for interpolation as we do for form assembly

[dham]

We should do that. However we can also avoid reassembling the matrices associated with a single interpolator object.

[wence-]

If we tie the lifetime of these things to the lifetime of the object in the forward model then that should DTRT.

[wence-]

Actually, that doesn't work

[ReubenHill]

Why doesn't it work?

[wence-]

Imagine you have the following setup:

def forward(mesh):
    # set up forward model
    interpolator = Interpolator(...)
    while ...:
        interpolator.interpolate()
    return something_that_is_not_the_interpolator

J = forward(mesh)

compute_gradient(J)

So the interpolator doesn't live past the end of the forward model run, but we need it in the adjoint if that's where we're stashing the other thing.

[ReubenHill]

Is this an edge case where stashing on the Interpolator would at least help?

Another idea, how about a global cache but only for the Interpolators in the InterpolateBlock to avoid the expense of caching all interpolators all the time?

[wence-]

Is this an edge case where stashing on the Interpolator would at least help?

No? The interpolator is the thing we're trying to keep alive.

Another idea, how about a global cache but only for the Interpolators in the InterpolateBlock to avoid the expense of caching all interpolators all the time?

How do you know this has happened?

[ReubenHill]

How do you know this has happened?

By doing the caching in the InterpolateBlock, then when a new InterpolateBlock is made, it looks in the cache of existing Interpolate objects for InterpolateBlocks to see if there's one it can reuse

[wence-]

How do those get collected?

[ReubenHill]

Well I don't know the details how caching is done at the moment, but I'm imagining some runtime-only global dictionary that's filled via an annotation of an Interpolate

[wence-]

what I mean is, how do you decide to throw those interpolate objects away from the interpolateblock?

[ReubenHill]

I don't think I understand your question? Do you mean getting rid of the cached objects? Presumably if they lived in a dictionary in some suitably global scope, they would only exist as long as a script took to run?

[wence-]

That sounds like too long. Each Interpolator, if it is reusable, carries a sparse matrix around with it.

[ReubenHill]

Solution: Cache all pyop2 Interpolation kernels created in firedrake/interpolation.py based on a hash of expression and function space (@wence- please correct if wrong)

[wence-]

Yes, I think that is right. ufl.core.Expr needs to gain a .signature() method that is similar to that for ufl.Form (which is stable wrt coefficient numbering). Then we can do something similar to the TSFCKernelCache stuff.

[ReubenHill]

Do function spaces have a suitable hash?

[wence-]

No, but their UFL element does (or should).