[pbrubeck]

Summary

TSFC now compiles same-mesh ufl.Interpolate nodes in normal Form kernels, and this PR ports the feature from PyOP2 to pyop3.

• Standalone Interpolate and Form(Interpolate) share one code generation path.
• Same-mesh interpolation across a MixedFunctionSpace supports mat_type="aij" natively.
• Interpolation into a VertexOnlyMesh uses the same TSFC code generation.
• Cross-mesh interpolation stays on BaseFormAssembler.

SameMeshInterpolator no longer builds its own op3.loops. It sends the interpolation through the form assembler:

 SameMeshInterpolator._get_callable
-  _build_interpolation_callables      # hand-built op3.loop per sub-expression
+  get_form_assembler(Interpolate, access=...)
+    tsfc_interface.compile_form(..., access)
+      as_pyop3_local_kernel(ast, nargs, access)
+        op3.Function(..., requires_zeroed_output_arguments=True)

MIN and MAX are reduced into the output as MIN_WRITE and MAX_WRITE. TSFC kernels add their result to the output argument, so pyop3 starts these temporaries from zero (see #5543).

TwoFormAssembler._apply_bc now takes the field label of a BC from the BC's own space. Before, a rectangular form such as an interpolation from W = V * V into V failed.

quadrature_space is now a cached property. The adjoint weight, the reciprocal of the number of cells in the iteration set whose closure contains each node, is cached on the mesh, keyed by the dual space and the subset.

Stacked on #5543, which adds MIN/MAX reductions to pyop3 kernel calls.
Depends on firedrakeproject/fiat#294.

Fixes #5468

Evidence

• Before (merge of main without the fixes): 7 failures in tests/firedrake/regression/test_interpolate.py

  test_subset_interpolation_reduction[min|max]     -> all entries stay at ±finfo.max
  test_interpolate_periodic_coords_reduction[...]  -> [inf]
  test_mixed_space_bcs                             -> MultipleUniqueEntriesException
  

  After:

  serial: interpolation, cross-mesh, vertex-only, MG grid transfer, pyop3 loop and TSFC tests  1006 passed
  serial: BC, mixed, real-space, matrix and assembly regression tests                           740 passed
  mpiexec -n 2  -m "parallel[2]" interpolation and vertex-only tests                            32 passed
  mpiexec -n 3  -m "parallel[3]" interpolation and vertex-only tests                           271 passed
  

test_subset_adjoint_interpolation (new, 1 and 2 ranks) compares the matrix-free adjoint with the transposed aij matrix, with and without a subset. It fails if the subset is left out of the cache key.

Repeated interpolation on BDM2 16×16 through Vq = V.quadrature_space, per call:

	before	after
Function(V).interpolate(q)	4.01 ms	3.87 ms
Cofunction(Vq.dual()).interpolate(r)	15.1 ms	7.8 ms

Merge Danger

Door: two-way

Blast Radius: interpolation

All same-mesh interpolation, and the Form kernels that contain it, now go through TSFC and the form assembler. op3.Function gets a new attribute, so kernel caches from before this change must be cleared with firedrake-clean. V.quadrature_space() is now the property V.quadrature_space.

The CI commits marked DROP BEFORE MERGE install the FIAT branch and must be dropped before merging.

Contribution assisted by Claude Code (Claude Opus 5.5).

🤖 Generated with Claude Code

[pbrubeck]

This PR should reimplement SameMeshInterpolator as a dispatcher to {Zero|One|Two}FormAssembler. We want to unify the TSFC codepaths for Form assembly and same-mesh interpolation.

[sghelichkhani]

Hi Pablo, I ran G-ADOPT on #5566 and repeated interpolations no longer build new loops, forward and adjoint. One case still does. A MIN/MAX interpolation into a new output builds a new loop on every call in _initialize_minmax through _loop_target_nodes. G-ADOPT does this every time step in compute_timestep and it costs 6 to 13 ms per call against 0.8 ms with pyop2. Could that loop be cached, or skipped when covers_target is true?

[pbrubeck]

Hi Pablo, I ran G-ADOPT on #5566 and repeated interpolations no longer build new loops, forward and adjoint. One case still does. A MIN/MAX interpolation into a new output builds a new loop on every call in _initialize_minmax through _loop_target_nodes. G-ADOPT does this every time step in compute_timestep and it costs 6 to 13 ms per call against 0.8 ms with pyop2. Could that loop be cached, or skipped when covers_target is true?

Done.

[sghelichkhani]

Thank you. For the record, in the G-ADOPT regression tests the loop cost about 6 ms per call against 0.7 ms before, which made each time step of forward cases such as 2d_compressible_TALA about 40 % slower. With a local change that writes the identity in the same way, the time per step was back to where it was and the results did not change.

[sghelichkhani]

A small heads-up. With current main merged in, this PR fails test_restricted_function_space.py::test_restrict_action[False] and test_poisson_restricted_mixed_space[Cofunction] with a KeyError on an Interpolate node. Both pass on main. I think it comes from #5260 and #5285 which were merged on 10 November.

[pbrubeck]

Thank you. For the record, in the G-ADOPT regression tests the loop cost about 6 ms per call against 0.7 ms before, which made each time step of forward cases such as 2d_compressible_TALA about 40 % slower. With a local change that writes the identity in the same way, the time per step was back to where it was and the results did not change.

What are the before and after here? What is the local change?

[sghelichkhani]

Before your fix one call of this interpolation took about 6 ms in a small serial benchmark of the expression in the G-ADOPT time step adaptor. With your fix it takes 0.5 ms. The code without this PR took 0.6 ms. It wrote the identity with f.assign(Constant(finfo.min)) and built no loop. My local change was a narrower version of your fix. When there was no subset and the source mesh was the target mesh it wrote the identity value into the whole array with NumPy and skipped the loop. Your fix covers that case and the submesh case so I have dropped it. In the G-ADOPT CI the 2d_compressible_TALA case took 474 s without this PR. It took 778 s with this PR before your fix and 452 s with your fix.

[connorjward]

I think that this is fundamentally there. Please do go over all of the tests and make sure that they actually test useful things and don't do the antipatterns I pointed out.

And in future please please please make smaller PRs. Reviewing this now just took me 40 minutes (and I have reviewed it a bunch of times before).

[connorjward]

raise what?

[pbrubeck]

when you raise inside an except, you raise whatever brought you there.

[connorjward]

This should just be an error no? Otherwise we are silently doing something the user doesn't expect and explicitly didn't ask for.

[pbrubeck]

The problem is the aliases for MIN == MIN_RW.

[connorjward]

I am open to redefining the alias. And whatever happens we're making a breaking change with pyop3. So we could just error if someone does op2.MIN -> op3.MIN and tell them to DTRT.

We could also allow this to accept strings ("min") instead of pyop3 intents. Intents in pyop3 are a more general concept than what we need here.

[pbrubeck]

I would rather open an issue and address it in a separate PR. For now, this does not break the API.

For now this is doing whatever the user expected with MIN before the introduction of MIN_RW / MIN_WRITE.

[connorjward]

Why can this be None? I would prefer something more explicit to disable assembler caching.

[pbrubeck]

The key was None for 0- and 2-forms (we do not cache assemblers for them). I added a _cache_assemblers property (True for 1-form, False otherwise).

[pbrubeck]

And in future please please please make smaller PRs. Reviewing this now just took me 40 minutes (and I have reviewed it a bunch of times before).

Thanks. I always aim to keep them short but something quite pervasive like this can't be split into multiple PRs, as you might understand (pyop3).