apply_geometry_lowering does not always lower Jacobians
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.
Comments (17)
-
reporter -
reporter - edited description
-
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.)
-
reporter 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.
-
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)
-
reporter Indeed, domain.coordinates() is None.
-
reporter 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. -
Try logging in the Domain constructor?
-
reporter 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.
-
reporter 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.
-
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.
-
reporter I wonder how you plan to solve this. Will you introduce circular references?
-
On the contrary. I will remove the circular dependency that is in todays design.
-
reporter Is it fixed yet?
-
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.
-
reporter Thanks.
-
- changed status to resolved
Should be resolved. Reopen if still relevant.
- Log in to comment
Adding @amcrae to watch.