BUG: Tape of mesh point dependencies

[GeorgAUT]

Dear all,

We are experiencing an issue with the taping of function dependencies on underlying mesh coordinates, which appear to not be recorded correctly. We were wondering if there is an easy fix for this? I believe the feature of differentiation with respect to mesh coordinates was available in previous versions of Firedrake (based on https://doi.org/10.1007/s00158-019-02281-z). Please find further details below - any help to resolve this would be really greatly appreciated.

Thank you in advance,
Georg

Steps to Reproduce & Expected behavior

The issue appears in the case of interpolating a function onto a mesh. A MWE is included below where we define the $\sin$ function by interpolation on the mesh. The gradient with respect to the mesh point should lead to $\cos$, i.e. the derivative values should be $[\cos(0),\cos(1/2),\cos(1)]\approx [1.0,0.8776,0.5403]$.

from firedrake.adjoint import *
from firedrake.__future__ import interpolate

# Set up the problem
continue_annotation()
mesh = UnitIntervalMesh(2)
spat_coord = SpatialCoordinate(mesh)

# Create function space for mesh coordinates and reassign
V_coord = mesh.coordinates.function_space()
coords = Function(V_coord).interpolate(spat_coord)
mesh.coordinates.assign(coords)

# Interpolate a scalar function (sin of spatial coordinate) into the scalar space
V_scalar = FunctionSpace(mesh, "CG", 1)
u_i = interpolate(sin(spat_coord[0]), V_scalar)
u = assemble(u_i)

# Compute gradient with respect to the mesh points
target = ReducedFunctional(u, Control(coords))
grad_target = target.derivative().dat.data.copy()
print(grad_target)

# Output tape
get_working_tape().visualise(output="tape_interpolation.pdf")

The output we receive here is [0. 0. 0.]. In the visualised tape (tape_interpolation.pdf) $w_3$ is coords, and $w_{15}$ is the function u. The dependency on the mesh points is again not tracked.

Error message
The dependency on mesh coordinates (control variable coords) appears to not be recorded on the tape.

/firedrake/lib/python3.12/site-packages/pytools/persistent_dict.py:63: RecommendedHashNotFoundWarning: Unable to import recommended hash 'siphash24.siphash13', falling back to 'hashlib.sha256'. Run 'python3 -m pip install siphash24' to install the recommended hash.
  warn("Unable to import recommended hash 'siphash24.siphash13', "
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance
WARNING:root:Adjoint value is None, is the functional independent of the control variable?
[0. 0. 0.]

Process finished with exit code 0

Environment:

• OS: MacOS 15.1.1
• Python version: 3.12.6
• Output of firedrake-status:

/Users/georgmaierhofer/Documents/firedrake/bin/firedrake-status:4: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  __import__('pkg_resources').require('firedrake==0.13.0+6193.g5a344e488')
Firedrake Configuration:
    package_manager: True
    minimal_petsc: False
    mpicc: None
    mpicxx: None
    mpif90: None
    mpiexec: None
    disable_ssh: False
    honour_petsc_dir: False
    with_parmetis: False
    slepc: False
    packages: []
    honour_pythonpath: False
    opencascade: False
    torch: False
    jax: False
    petsc_int_type: int32
    cache_dir: /Users/georgmaierhofer/Documents/firedrake/.cache
    complex: False
    remove_build_files: False
    with_blas: download
    netgen: False
Additions:
    None
Environment:
    PYTHONPATH: None
    PETSC_ARCH: None
    PETSC_DIR: None
Status of components:
---------------------------------------------------------------------------
|Package             |Branch                        |Revision  |Modified  |
---------------------------------------------------------------------------
|FInAT               |master                        |e72e533   |False     |
|PyOP2               |master                        |b84b7708  |False     |
|fiat                |master                        |78fc962   |False     |
|firedrake           |master                        |5a344e488 |False     |
|h5py                |firedrake                     |db2ab02f  |False     |
|libsupermesh        |master                        |c07de13   |False     |
|loopy               |main                          |d9876d83  |False     |
|petsc               |firedrake                     |3a8a2250205|False     |
|pyadjoint           |master                        |b9574fb   |False     |
|pytest-mpi          |main                          |f4a9692   |False     |
|tsfc                |master                        |9f42f95   |False     |
|ufl                 |master                        |c1a8afb1  |False     |
---------------------------------------------------------------------------

[GeorgAUT]

Small update: the first point raised in the original post was related to the representation of the derivative and the expected result for the gradient of a mass matrix entry with respect to the mesh coordinates can be found by calling .derivative(options={"riesz_representation": "l2"}). Thus I have removed this from the above. The remaining issue appears to really be that the interpolate block is not recorded on the tape.

[dham]

That is correct. This is related to the change to the new interpolate. Basically, AssembleBlock doesn't see the coordinates as a dependency of interpolate. I'm not yet sure how to fix that.

[APaganini]

I think the old interpolate had the same "issue" (or, rather, feature) of not tracking the dependency on the interpolee in the interpolant coefficients. E.g.,

from firedrake import *

mesh = UnitIntervalMesh(2)
x = SpatialCoordinate(mesh)
f = sin(x[0])

V = FunctionSpace(mesh, "CG", 1)
u = Function(V).interpolate(f)

v = TestFunction(VectorFunctionSpace(mesh, "CG", 1)) 

du = assemble(derivative(u*dx , x, v)) 
print("shape derivative of interpolant", du.vector().array())

df = assemble(derivative(f*dx , x, v)) 
print("shape derivative of interpolee", df.vector().array())

du_analytic = assemble(dot(grad(u), v) * dx + u*div(v)*dx)
print("analytic fomrmula", du_analytic.vector().array())

output:

shape derivative of interpolant [-0.23971277 -0.42073549  0.66044826]
shape derivative of interpolee [-1.15157951e-08 -2.02121220e-08  8.41471023e-01]
analytic fomrmula [0.         0.         0.84147098]