apply_geometry_lowering does not always lower Jacobians

[Miklós Homolya]

Latest UFL changes therefore also break Firedrake "bendy" (non-affine support).

from firedrake import *

mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, "RT", 1)

u = TrialFunction(V)
v = TestFunction(V)
L = assemble(div(u)*div(v)*dx)

gives:

Compiling form form

Compiler stage 1: Analyzing form(s)
-----------------------------------

  Geometric dimension:       2
  Number of cell subdomains: 0
  Rank:                      2
  Arguments:                 '(v_0, v_1)'
  Number of coefficients:    1
  Coefficients:              '[Coordinates]'
  Unique elements:           'RT1(?), Vector<2 x CG1(?)>'
  Unique sub elements:       'RT1(?), Vector<2 x CG1(?)>, CG1(?)'

  representation:    quadrature
  quadrature_degree: auto --> 3
  quadrature_rule:   auto --> default

Compiler stage 1 finished in 0.0809739 seconds.

Compiler stage 2: Computing intermediate representation
-------------------------------------------------------
  Computing representation of integrals
  Computing quadrature representation
  Transforming cell integral
This object should be implemented by the child class.
Traceback (most recent call last):
  File "/home/mh1714/Documents/spam.py", line 9, in <module>
    L = assemble(div(u)*div(v)*dx)
  ...
  File "/home/mh1714/git/ffc/ffc/quadrature/optimisedquadraturetransformer.py", line 349, in jacobian
    error("This object should be implemented by the child class.")
  File "<string>", line 1, in <lambda>
  File "/home/mh1714/git/ufl/ufl/log.py", line 151, in error
    raise self._exception_type(self._format_raw(*message))
Exception: This object should be implemented by the child class.

Some explanation: ufl/algorithms/compute_form_data.py:227 introduces Jacobians into the form which are not removed by ufl/algorithms/compute_form_data.py:231 and therefore cause FFC to fail.

[Miklós Homolya]

Adding @amcrae to watch.

[Martin Sandve Alnæs]

I don't see how that can happen. Are you sure you call compute_form_data with do_apply_geometry_lowering=True?

(Also I'm always initially confused when you say "FFC" and mean "the firedrake version of FFC", since those are very different.)

[Miklós Homolya]

I'm sure:

    # Compute form metadata
    form_data = compute_form_data(form,
                                  do_apply_function_pullbacks=True,
                                  do_apply_integral_scaling=True,
                                  do_apply_geometry_lowering=True,
                                  preserve_geometry_types=(),
                                  do_apply_restrictions=True,
                                  )

in ffc/analysis.py from mapdes/ffc/new_fd_bendy https://bitbucket.org/mapdes/ffc/src/ebc1307dd41bbbcbd0b6ec01b00d543055b7bb4f/ffc/analysis.py?at=new_fd_bendy

I was stepping through the above in the debugger.

[Martin Sandve Alnæs]

Here's the jacobian handler. The only possibility I see is that "domain.coordinates() is None" must be true. Can you add a check in there to see? (apply_geometry_lowering.py line 75)

    @memoized_handler
    def jacobian(self, o):
        if self._preserve_types[o._ufl_typecode_]:
            return o

        domain = o.domain()
        if domain.coordinates() is None:
            # Affine case in FEniCS: preserve J if there's no coordinate function
            # (the handling of coordinate functions will soon be refactored)
<add check here>
            return o

        x = self.spatial_coordinate(SpatialCoordinate(domain))
        return ReferenceGrad(x)

[Miklós Homolya]

Indeed, domain.coordinates() is None.

[Miklós Homolya]

Why might this (domain.coordinates() being None) happen though? If I use VectorFunctionSpace CG instead of RT, or change the form to e.g. dot(u, v)*dx, this does not fail.

[Martin Sandve Alnæs]

Try logging in the Domain constructor?

[Miklós Homolya]

Only one domain gets constructed without coordinates before the failure, the coordinates domain, which indeed should not have coordinates.

I suspect that it somehow wrongly leaks into the form during the transformations.

[Miklós Homolya]

I think I have figured out the mechanics of the failure, but I'm not sure how to fix it.

When apply_geometry_lowering is called for the first time, it translates any SpatialCoordinate into a ReferenceValue of the coordinate coefficient whose domain has no coordinates. Then apply_derivates translates the ReferenceValue into something containing the JacobianInverse of this coordinateless domain. When apply_geometry_lowering is called again, jacobian will be called on the domain with no coordinates, hence no lowering will happen.

[Martin Sandve Alnæs]

This type of confusion will be avoided when I redesign Domain a bit and remove Domain.coordinates(). Sounds like that should happen sooner rather than later.

[Miklós Homolya]

I wonder how you plan to solve this. Will you introduce circular references?

[Martin Sandve Alnæs]

On the contrary. I will remove the circular dependency that is in todays design.

[Miklós Homolya]

Is it fixed yet?

[Martin Sandve Alnæs]

I believe it's no longer relevant. There is no longer such a thing as the coordinate Coefficient of a Domain, instead the Mesh has only a coordinate element. To refer to the coordinate field symbolically, SpatialCoordinate(mesh) should be used.

[Miklós Homolya]

Thanks.

[Martin Sandve Alnæs]

Should be resolved. Reopen if still relevant.