Convergence analysis using a reference solution

[fritz-io]

What is the recommended approach in Firedrake for computing convergence rates with H²-conforming elements (e.g., Bell/Argyris) when using a reference solution instead of an analytical solution? I am thinking of a time-dependent biharmonic problem when no analytical solution is known.
When attempting to compare solutions across different mesh resolutions for convergence analysis, I encounter domain ambiguity errors (ValueError: Found multiple domains) during projection/interpolation. This occurs even when using MeshHierarchy and variational projection.

from firedrake import *

base_mesh = UnitSquareMesh(8, 8, quadrilateral=False)
hierarchy = MeshHierarchy(base_mesh, 1) # 2 levels: 8x8, 16x16

fine_mesh = hierarchy[-1]
V_fine = FunctionSpace(fine_mesh, "Bell", 5)
u_ref = Function(V_fine).project(sin(pi*SpatialCoordinate(fine_mesh)[0]))

coarse_mesh = hierarchy[0]
V_coarse = FunctionSpace(coarse_mesh, "Bell", 5)
u_h = Function(V_coarse).project(sin(pi*SpatialCoordinate(coarse_mesh)[0]))

error = errornorm(u_ref, u_h, 'H1') # ValueError: Multiple domains

What I tried:

1. Variational projection:

v_coarse = TestFunction(V_coarse)
u_proj = TrialFunction(V_coarse)
solve(inner(u_proj, v_coarse)*dx == inner(u_ref, v_coarse)*dx, ...)

Fails with Found multiple domains.

2. interpolate with Spatial Coordinates:

x_coarse = SpatialCoordinate(coarse_mesh)
u_ref_on_coarse = interpolate(u_ref(x_coarse[0], x_coarse[1]), V_coarse)

Fails with Indexed.float error.

3. Mesh Hierarchy Tools: With FunctionHierarchy

restrict(u_ref, u_coarse)

AttributeError/Failed imports

What is the canonical way to transfer solutions between non-nested meshes for H² elements? For convergence studies without analytical solutions, what would you recommend?
I use Firedrake v0.13.0 and PETSc v3.16.6 on a Macbook Air M1.