Initial support for interpolation

[jorgensd]

Feature: Differentiable Function Interpolation Block

AI disclosure statement

This pull-request has been developed by me, with the help of Gemini (august 2026).

Description

This PR introduces the InterpolationBlock and an overloaded dolfinx_adjoint.interpolate wrapper, enabling automatic differentiation through dolfinx.fem.Function.interpolate. It supports both PETSc.Mat and native dolfinx.la.MatrixCSR backends while ensuring robust parallel execution and minimal memory overhead during optimization.

Mathematical Formulation

Because function interpolation is a purely linear mapping between finite element spaces, the forward, adjoint, and Hessian operations are governed entirely by the interpolation matrix and its transpose.

1. Forward and Tangent Linear Model (TLM)
Let $\Pi$ denote the interpolation matrix mapping from a source space $V$ to a target space $W$. For an input coefficient vector $x \in V$, the interpolated vector $y \in W$ is given by:

$$y = \Pi x$$

Because the operator is linear, the Tangent Linear Model (the Jacobian action) applied to a perturbation $\delta x$ is simply the interpolation matrix itself:

$$\delta y = \Pi \delta x$$

2. Adjoint Model
The adjoint of a linear operator propagates downstream sensitivities backwards. Given the sensitivity of a functional with respect to the output, $y^* = \frac{\partial J}{\partial y}$, the adjoint identity $\langle \Pi \delta x, y^* \rangle_W = \langle \delta x, \Pi^T y^* \rangle_V$ dictates that the sensitivity with respect to the input is the transpose of the interpolation matrix:

$$x^* = \Pi^T y^*$$

3. Hessian Vector Product
While the second derivative of the linear interpolation operator is zero, the chain rule dictates how downstream curvature is mapped backwards. Given a downstream Hessian vector product $H_{yy} \delta y$, the incoming Hessian sensitivity is $\delta y^* = H_{yy} \Pi \delta x$. Our block completes the chain rule by applying the adjoint operator to this incoming sensitivity, yielding the full Hessian vector product for the control variable:

$$H_{xx} \delta x = \Pi^T (H_{yy} \Pi \delta x) = \Pi^T \delta y^*$$

Implementation Highlights

• Dual Backend Support: The get_mult closure abstracts away the differences between PETSc and native dolfinx.la.MatrixCSR operations.
• Parallel Safety: When using the native CSR backend, the block explicitly zeroes out ghost degrees of freedom in the input vector prior to transpose multiplication. This physically prevents the double-counting of ghosted row contributions across MPI boundaries.
• Memory Efficiency: The interpolation matrices are globally cached using space ID tuples (_INTERPOLATION_MATRIX_CACHE). Additionally, persistent working C++ arrays (_row_vec and _col_vec) are cached on the matrix object, eliminating dynamic Python allocations during iterative line searches.
• Taylor Tested: Verified to achieve Taylor remainder convergence rates of 1.0, 2.0, and 3.0 across continuous, discontinuous, and Nedelec elements in parallel.