Error when taking the derivative of a form which includes the cell-normals of a manifold.

[JHopeCollins]

When trying to create a LinearVariationalProblem from the derivative of a form, UFL throws an error if:

• The mesh is a manifold
• The coordinate degree is greater than 1
• The form includes the gradient of the cell normals

MFE:

import firedrake as fd

coords_degree = 2  # works when coords_degree=1, breaks when coords_degree=2

mesh = fd.IcosahedralSphereMesh(radius=1,
                                refinement_level=2,
                                degree=coords_degree)
x = fd.SpatialCoordinate(mesh)
mesh.init_cell_orientations(x)
n = fd.CellNormal(mesh)

V = fd.FunctionSpace(mesh, "BDM", 1)

u = fd.Function(V)
v = fd.TestFunction(V)

F = fd.inner(fd.grad(fd.inner(v, n)), u)*fd.dx

A = fd.derivative(F, u)

f = fd.Function(V)
b = fd.inner(v, f)*fd.dx

prob = fd.LinearVariationalProblem(A, b, u)  # error here

Error message:

UFL:ERROR Currently no support for ReferenceGrad in CoefficientDerivative.
Traceback (most recent call last):
File "refgrad_test.py", line 30, in
prob = fd.LinearVariationalProblem(A, b, u) # error here
File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/home/jhc/icl_ra/firedrake/src/firedrake/firedrake/variational_solver.py", line 316, in init
F = ufl_expr.action(J, u) - L
File "PETSc/Log.pyx", line 115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "PETSc/Log.pyx", line 116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/home/jhc/icl_ra/firedrake/src/firedrake/firedrake/ufl_expr.py", line 212, in action
return ufl.action(form, coefficient)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/formoperators.py", line 109, in action
form = expand_derivatives(form)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/ad.py", line 34, in expand_derivatives
form = apply_derivatives(form)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/apply_derivatives.py", line 1136, in apply_derivatives
return map_integrand_dags(rules, expression)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 46, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 27, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 27, in
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 35, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 46, in
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/corealg/map_dag.py", line 36, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/corealg/map_dag.py", line 99, in map_expr_dags
r = handlers[v.ufl_typecode](v, *[vcache[u] for u in v.ufl_operands])
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/apply_derivatives.py", line 1090, in coefficient_derivative
return map_expr_dag(rules, f,
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/corealg/map_dag.py", line 36, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/corealg/map_dag.py", line 97, in map_expr_dags
r = handlersv.ufl_typecode
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/algorithms/apply_derivatives.py", line 889, in reference_grad
error("Currently no support for ReferenceGrad in CoefficientDerivative.")
File "/home/jhc/icl_ra/firedrake/src/ufl/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Currently no support for ReferenceGrad in CoefficientDerivative.

This seems to have only started happening since the UFL update, although I'm not sure of exactly when.

[ReubenHill]

I have recreated the issue. Based on what @tommbendall said about last week's UFL update being the cause, I reset UFL back to 1dddf46eb5d011e7144d8780ac008dc2a87a2ce6 and ran the script without trouble. So it seems that either something that got added to FEniCS UFL since then, or indeed @rckirby 's accidental PR merge to Firedrake UFL is to blame. I think we ought to be able to find the culprit by looking at the diff for firedrakeproject/ufl#30?

[wence-]

This is the culprit:

https://github.com/firedrakeproject/ufl/pull/30/files#diff-911adc50ed9ac0d171fe7ba0b3dc8be35e4f9574d1bd4939be9f956b3df2f138R652-R655

Previously the default was that a geometric quantity was considered to be is_cellwise_constant() => True (aside: this is an unsafe default to apply to an abstract base class, but meh).

That merge brings in a fix for CellNormal (which is only cellwise constant if the element is a linear simplex). But that means that now we have new things appearing in preprocessed forms. So one needs to figure out what the correct derivative lowering is in this case. It would be good to simplify the form. I think something like this should provoke the issue?

from tsfc.ufl_utils import compute_form_data
...

n = CellNormalI(mesh)
u = Function(V)
phi = grad(inner(u, n))[0]*dx
fd = compute_form_data(derivative(phi, u))

?

[michalhabera]

See FEniCS/ufl#119. There was a bug in UFL which incorrectly simplified all CellNormal derivatives to 0.

[michalhabera]

I guess TSFC knows how to compute derivatives of Jacobian, (not needed, below is anyway lowered to reference_grad()) The following premature CellNormal lowering could help you?

import ufl

mesh = ufl.Mesh(ufl.VectorElement("CG", ufl.triangle, 2, dim=3))
cn = ufl.algorithms.apply_geometry_lowering.apply_geometry_lowering(ufl.CellNormal(mesh))

gcn = ufl.grad(cn)

[JHopeCollins]

@wence-, thanks for finding that. Your simpler form does result in the same error too. Full script:

import firedrake as fd                                                                           
from tsfc.ufl_utils import compute_form_data

coords_degree = 2  # works when coords_degree=1, breaks when coords_degree=2

mesh = fd.IcosahedralSphereMesh(radius=1,
                                refinement_level=2,
                                degree=coords_degree)
x = fd.SpatialCoordinate(mesh)
mesh.init_cell_orientations(x)
n = fd.CellNormal(mesh)

V = fd.FunctionSpace(mesh, "BDM", 1)

u = fd.Function(V)

phi = fd.grad(fd.inner(u, n))[0]*fd.dx

form_data = compute_form_data(fd.derivative(phi, u)) 

@michalhabera, thanks for the workaround, I'll see if we can use it for our form for the time being.

[wence-]

(That workaround appears not to help :()

Here's a minimal example that only uses UFL to demonstrate the issue:

from ufl import (Cell, CellNormal, Coefficient, FunctionSpace, Mesh,
                 VectorElement, derivative, dx, grad, inner)
from ufl.algorithms import compute_form_data
from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering

mesh = Mesh(VectorElement("P", Cell("triangle", 3), 2))

try_workaround = False
if try_workaround:
    n = apply_geometry_lowering(CellNormal(mesh))
else:
    n = CellNormal(mesh)

V = FunctionSpace(mesh, VectorElement("P", mesh.ufl_cell(), 2))

u = Coefficient(V)

phi = grad(inner(u, n))[0]*dx

fd = compute_form_data(
    derivative(phi, u),
    do_apply_function_pullbacks=True,
    do_apply_geometry_lowering=True,
)

[michalhabera]

(That workaround appears not to help :()

Here's a minimal example that only uses UFL to demonstrate the issue:

from ufl import (Cell, CellNormal, Coefficient, FunctionSpace, Mesh,
                 VectorElement, derivative, dx, grad, inner)
from ufl.algorithms import compute_form_data
from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering

mesh = Mesh(VectorElement("P", Cell("triangle", 3), 2))

try_workaround = False
if try_workaround:
    n = apply_geometry_lowering(CellNormal(mesh))
else:
    n = CellNormal(mesh)

V = FunctionSpace(mesh, VectorElement("P", mesh.ufl_cell(), 2))

u = Coefficient(V)

phi = grad(inner(u, n))[0]*dx

fd = compute_form_data(
    derivative(phi, u),
    do_apply_function_pullbacks=True,
    do_apply_geometry_lowering=True,
)

Indeed. I now remember in all manifold codes we first interpolate the normal field into a conforming VectorFunction space, this avoids the CoefficientDerivative of ReferenceGrad.

Alternatively, and dangerously, you could make coefficient derivatives of ReferenceGrad = 0 by monkey-patching UFL,

from ufl.algorithms.apply_derivatives import GateauxDerivativeRuleset, GenericDerivativeRuleset

def myrule(self, o):
    print(f"D({o})/D({self._w}) simplified to 0")
    return GenericDerivativeRuleset.independent_terminal(self, o)

GateauxDerivativeRuleset.reference_grad = myrule

[dham]

The comments in apply_derivatives.py appear to indicate that CoefficientDerivative of ReferenceGrad could be implemented for the simple case that causes this error, and could be allowed to still fail in more complex cases. Doing the same thing for ReferenceValue is probably advisable.

[jshipton]

Could I ask for an update on this issue please - I'm afraid I got rather confused trying to follow the discussions referenced above. I'm really hoping to be able to run the vector invariant form of the shallow water equations on the sphere again!

[wence-]

Correct fix: handle derivatives of referencegrad correctly in UFL.

Workaround:

# a quadratic approximation of the cell normal (or whatever you need)
# if you want the previous affine approximation of cell normals you could use VectorFunctionSpace(mesh, "DG", 0) instead.
cn = interpolate(CellNormal(mesh), VectorFunctionSpace(mesh, "DG", 2))

# use cn as normal (hah!)

[blechta]

Note that the same issue is seen for a facet normal of a higher-order cell, thus not depending on manifolds, see #3313.