BUG: CellDiameter on extruded meshes no longer works

[stephankramer]

The following code

from firedrake import *
mesh1d = UnitIntervalMesh(2)
mesh = ExtrudedMesh(mesh1d, 2)
DG0 = FunctionSpace(mesh, "DG", 0)
assemble(CellDiameter(mesh)*TestFunction(DG0)*dx)

produces the following error

/home/skramer/fd/src/ufl/ufl/algorithms/apply_geometry_lowering.py:310: UserWarning: Only know how to compute cell diameter of P1 or Q1 cell.
  warnings.warn("Only know how to compute cell diameter of P1 or Q1 cell.")
Traceback (most recent call last):
  File "/home/skramer/fd/src/PyOP2/pyop2/caching.py", line 198, in __new__
    return cls._cache_lookup(key)
           ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/tsfc_interface.py", line 73, in _cache_lookup
    return cls._cache.get((key, commkey)) or cls._read_from_disk(key, comm)
                                             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/tsfc_interface.py", line 94, in _read_from_disk
    raise KeyError(f"Object with key {key} not found")
KeyError: 'Object with key 43358d716dc6302816532f87f1d29b60 not found'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/skramer/tst.py", line 5, in <module>
    assemble(CellDiameter(mesh)*TestFunction(DG0)*dx)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/skramer/fd/src/firedrake/firedrake/adjoint_utils/assembly.py", line 30, in wrapper
    output = assemble(form, *args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 133, in assemble
    return get_assembler(expr, *args, **kwargs).assemble(tensor=tensor)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 981, in assemble
    self.execute_parloops(tensor)
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 1215, in execute_parloops
    for parloop in self.parloops(tensor):
                   ^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 1011, in parloops
    self._parloops = tuple(parloop_builder.build(tensor) for parloop_builder in self.parloop_builders)
                                                                                ^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/PyOP2/pyop2/utils.py", line 62, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
                                           ^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 1017, in parloop_builders
    for local_kernel, subdomain_id in self.local_kernels:
                                      ^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/PyOP2/pyop2/utils.py", line 62, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
                                           ^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/firedrake/firedrake/assemble.py", line 1052, in local_kernels
    kernels = tsfc_interface.compile_form(
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/skramer/fd/src/firedrake/firedrake/tsfc_interface.py", line 261, in compile_form
    kinfos = TSFCKernel(f, prefix, parameters,
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/PyOP2/pyop2/caching.py", line 200, in __new__
    obj = make_obj()
          ^^^^^^^^^^
  File "/home/skramer/fd/src/PyOP2/pyop2/caching.py", line 190, in make_obj
    obj.__init__(*args, **kwargs)
  File "/home/skramer/fd/src/firedrake/firedrake/tsfc_interface.py", line 146, in __init__
    tree = tsfc_compile_form(form, prefix=name, parameters=parameters,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/driver.py", line 86, in compile_form
    kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, diagonal=diagonal, log=log)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/driver.py", line 152, in compile_integral
    integrand_exprs = builder.compile_integrand(integrand, params, ctx)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/kernel_interface/common.py", line 142, in compile_integrand
    expressions = fem.compile_ufl(integrand,
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/fem.py", line 721, in compile_ufl
    result = map_expr_dags(context.translator, expressions)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/ufl/ufl/corealg/map_dag.py", line 98, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/ufl_utils.py", line 162, in _modified_terminal
    return self.modified_terminal(o)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/fem.py", line 365, in modified_terminal
    return translate(mt.terminal, mt, self.context)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3.11/functools.py", line 909, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/skramer/fd/src/tsfc/tsfc/fem.py", line 387, in translate_geometricquantity
    raise NotImplementedError("Cannot handle geometric quantity type: %s" % type(terminal))
NotImplementedError: Cannot handle geometric quantity type: <class 'ufl.geometry.CellDiameter'>

This does seem to work on an older Firedrake installation (about Jul '23).

[connorjward]

This is very curious. I have looked at the commits for TSFC and UFL in the last year and none of them seem to touch CellDiameter.

[wence-]

My guess is https://github.com/FEniCS/ufl/pull/197

That probably changed the implemented of cell_diameter in apply_geometry_lowering in UFL to return CellDiameter for extruded hex cells, rather than fall-through to the Q1 case.

[connorjward]

Spot on @wence-, thanks! For an extruded mesh we have domain.ufl_coordinate_element().degree() == (1, 1) which was a previously permissible element. Now though we have domain.ufl_coordinate_element().embedded_superdegree == 2 which is no longer allowed.

This might well be a simple fix. Looking at where embedded_superdegree is defined we should probably be using max instead of sum (we even use min for embedded_subdegree).

[wence-]

Nah, sum is correct, since a Q1 element contains some quadratic polynomials, so it is embedded in the complete polynomial space P2.

[wence-]

That implementation is, however, probably wrong because it should probably say sum(sub.embedded_superdegree() for sub in sub_elements()) (which would handle the case of TPE(TPE(interval, interval), interval) (say)

[connorjward]

OK. Do you have any suggestion for what the right thing is to do here? Should CellDiameter for an extruded mesh not be supported?

[wence-]

I think in apply_geometry_lowering the check for Q1 cells needs to be reinstated. So it does something like:

if all(degree == 1 for degree in self.degree()): 
    # can do affine and/or Q1 lowering

I don't know exactly how to check for Q1, but notice that previously the check was (pseudo-code)

if not (P1 or Q1):
    dont_lower
elif P1:
    lower_simplex
else:
    # must be Q1
    lower_hypercube

But now the check is:

if not P1:
    dont_lower
elif P1:
    lower_simplex
else:
    # unreachable
    lower_hypercube

[connorjward]

I believe that degree() is no longer a UFL finite element property and now only embedded_{super,sub}degree exist so this may be tricky to fix. I have added this to the agenda for the next Firedrake meeting.

[dham]

This looks like a mistake in the UFL changes last year. The test for geometry lowering is using the embedded_superdegree but it should be using the embedded_subdegree. The reason subdegree is right is that you care about whether the edges are straight, not whether there are any quadratic functions on the interior.

[connorjward]

Fixed in https://github.com/FEniCS/ufl/pull/295 (though we also need to update our UFL fork).