BUG: spatial derivative of facet normal fails ceefficient derivative

[blechta]

Coefficient derivative of a form featuring a spatial derivative of facet normal fails in form preprocessing. This is possibly a UFL and/or TSFC issue, but filing here a Firedrake MWE for reference

MWE:

from firedrake import *

# Create mesh with bendy triangles
mesh = UnitSquareMesh(1, 1)
x, y = SpatialCoordinate(mesh)
coord_space = VectorFunctionSpace(mesh, "P", 2)
coords = interpolate(SpatialCoordinate(mesh), coord_space)
mesh = Mesh(coords)

# Create form featuring spatial derivatives of facet normal
V = FunctionSpace(mesh, "P", 1)
u = Function(V)
v = TestFunction(V)
n = FacetNormal(mesh)
F = (1+n[0].dx(0))*u*v*ds - v*dx

# The culprit: computing coefficient derivative of n[0].dx(0)
J = derivative(F, u)

# This works: no spatial derivative of normal vector in coefficient derivative
#J = derivative(u*v*ds, u)

# Trigger Jacobian preprocessing
solve(F == 0, u, J=J)

Error message:

[...]
  File "/home/firedrake/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 241, in preprocess_form
    form = apply_derivatives(form)
[...]
  File "/home/firedrake/firedrake/src/ufl/ufl/algorithms/apply_derivatives.py", line 979, in reference_grad
    raise NotImplementedError("Currently no support for ReferenceGrad in CoefficientDerivative.")
NotImplementedError: Currently no support for ReferenceGrad in CoefficientDerivative.

Changing the coordinate degree to 1 or supplying the Jacobian without n avoids the problem, thus demonstrating the culprit.

Expected behavior: This error should not happen. Derivative of grad(n) w.r.t. u should simplify to zero.

Environment: This has been reproduced with Firedrake 0.13.0+5906.g7eaa7842d (Docker firedrakeproject/firedrake:2023-11) as well as 0.13.0+5994.gc5e939dde (Docker firedrakeproject/firedrake:latest of January 8, 2024, digest d22b43594971).

[wence-]

I think this is just purely a UFL issue. apply_geometry_lowering will turn n into some function of the jacobian hitting the reference normal, which in turn produces a ReferenceGrad(SpatialCoordinate(...)), so we have CoefficientDerivative(...RGrad(X), u) and this eventually fails.

What's surprising to me is that this doesn't show up in simpler examples.

What about:

x, *_ = SpatialCoordinate(mesh)

F = (1 + x.dx(0))*u*v*dx
J = derivatve(F, u)

?

[michalhabera]

See also https://github.com/firedrakeproject/firedrake/issues/2595, all UFL only related.

[blechta]

What's surprising to me is that this doesn't show up in simpler examples.

What about:

x, *_ = SpatialCoordinate(mesh)

F = (1 + x.dx(0))*u*v*dx
J = derivatve(F, u)

?

I suppose that x.dx(0) gets simplified to 1 here.