search

firedrakeproject / tsfc

Two-stage form compiler
Other
15 stars 25 forks source link

Certain forms on facets breaking _select_expression in optimise.py on assembly #274

Closed BorisAndrews closed 2 years ago

BorisAndrews commented 2 years ago

Certain forms arising from attempting to specify Lagrange multipliers on facets causes _select_expression in optimise.py to return NotImplementedError: No rule for factorising expressions of this kind..

Example code to return this error:

from firedrake import *
​
degree = 3
mesh_size = 8
​
mesh = UnitSquareMesh(mesh_size, mesh_size, quadrilateral = True)
n = FacetNormal(mesh)
​
cell = mesh.ufl_cell()
UL = FunctionSpace(mesh, MixedElement([
    FiniteElement("RTCE", cell, degree),
    RestrictedElement(FiniteElement("Q", cell, degree), restriction_domain = "facet")
  ]))
​
ul = Function(UL)
u, l = split(ul)
​
energy = 2*avg(inner(u, n))*l('+')*dS
deriv = derivative(energy, ul, TestFunction(UL))
vec = assemble(deriv)

Work around found by @wence-: setting {"mode" : "vanilla"} in the form_compiler_parameters of the specific problem, alongside adding return partial_indexed(ListTensor(expressions), (index,)) to _select_expression before the error is raised, resolved the issue.

wence- commented 2 years ago

Here's a reduced example of the problem:

from finat.point_set import PointSet
from gem import Index, select_expression
from tsfc.finatinterface import create_element
from ufl import FiniteElement, RestrictedElement, quadrilateral

ufl_element = RestrictedElement(FiniteElement("Q", quadrilateral, 2), restriction_domain="facet")
ps = PointSet([[0.5]])
finat_element = create_element(ufl_element)
evaluations = []
for eid in range(4):
    (val,) = finat_element.basis_evaluation(0, ps, (1, eid)).values()
    evaluations.append(val)
expr = select_expression(evaluations, Index())

The idea behind select_expression is that it is meant to push the indexing inside by appropriately factorising a list of expressions. In this case the parts that end up failing are: [Product(...), Zero()] which isn't a handled set of cases. I think that what is happening is that there is an assumption that each of the expressions in the list handed to select_expression have the same structure (in the sense of pattern-matching). In this case however, the tabulation of some parts of the element result in literal Zeros which destroys the structure-equivalence.

One can get around this by not eagerly producing symbolic zeros when building Literal objects and constant folding thereafter. But I don't know if this has any negative consequences down the line.

wence- commented 2 years ago

That merge had bad consequences, so we've reverted it for now.