FEniCS / ffcx

Next generation FEniCS Form Compiler for finite element forms
https://fenicsproject.org
Other
146 stars 39 forks source link

FFCx generates multiple copies of the same kernel #716

Open jpdean opened 2 days ago

jpdean commented 2 days ago

Running FFCx on the following code

from basix.ufl import element
from ufl import (
    FunctionSpace,
    Mesh,
    TestFunction,
    TrialFunction,
    dx,
    grad,
    inner,
)

n = 10

e = element("Lagrange", "triangle", 1)
coord_element = element("Lagrange", "triangle", 1, shape=(2,))
mesh = Mesh(coord_element)

V = FunctionSpace(mesh, e)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(grad(u), grad(v)) * dx(tuple(i for i in range(n)))

generates n copies of the same kernel.

garth-wells commented 2 days ago

This was fixed in https://github.com/FEniCS/ufl/pull/92, but must have crept back in.

jorgensd commented 2 days ago

This is rather strange. Whenever we introduce a derivative in any of these integrals, we do not get the kernel reduction:

from ufl import TrialFunction, TestFunction, dx, Mesh, triangle, FunctionSpace, inner, grad
import ufl.algorithms
import ufl.classes
from ufl.finiteelement import FiniteElement
from ufl.sobolevspace import H1
from ufl.pullback import identity_pullback
cell = triangle
mesh = Mesh(FiniteElement("Lagrange", cell, 1, (2,), identity_pullback, H1))
V = FunctionSpace(mesh, FiniteElement("Lagrange", cell, 1, (), identity_pullback, H1))
u = TrialFunction(V)
v = TestFunction(V)
n = 10
ids = tuple(i for i in range(n))
a = inner(u.dx(0),v)*ufl.dx(ids)
form_data = ufl.algorithms.compute_form_data(
    a,
    do_apply_function_pullbacks=True,
    do_apply_integral_scaling=True,
    do_apply_geometry_lowering=True,
    preserve_geometry_types=(ufl.classes.Jacobian,),
    do_apply_restrictions=True,
    do_append_everywhere_integrals=False)
for i, integral in enumerate(form_data.integral_data):
    print(i, str(integral))

While if we change a to a = inner(u,v)*ufl.dx(ids) we get only one integral

jorgensd commented 2 days ago

Seems to be a fundamental issue in UFL after computing form data. Here is an example

from ufl import TestFunction, Mesh, interval, FunctionSpace
import ufl.algorithms
import ufl.classes
from ufl.finiteelement import FiniteElement
from ufl.sobolevspace import H1
from ufl.pullback import identity_pullback
cell = interval
mesh = Mesh(FiniteElement("Lagrange", cell, 1, (2, ), identity_pullback, H1))
V = FunctionSpace(mesh, FiniteElement("Lagrange", cell, 1, (), identity_pullback, H1))
u = TestFunction(V)
a0 = u.dx(0) * ufl.dx((1))
form_data0 = ufl.algorithms.compute_form_data(
    a0,
    do_apply_function_pullbacks=True,
    do_apply_integral_scaling=True,
    do_apply_geometry_lowering=True,
    preserve_geometry_types=(ufl.classes.Jacobian,),
    do_apply_restrictions=True,
    do_append_everywhere_integrals=False)
form_data1 = ufl.algorithms.compute_form_data(
    a0,
    do_apply_function_pullbacks=True,
    do_apply_integral_scaling=True,
    do_apply_geometry_lowering=True,
    preserve_geometry_types=(ufl.classes.Jacobian,),
    do_apply_restrictions=True,
    do_append_everywhere_integrals=False)
i0 = form_data0.integral_data[0].integrals[0].integrand()
i1 = form_data1.integral_data[0].integrals[0].integrand()
print(i0 == i1)

This yields False. @mscroggs any ideas?