search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
498 stars 157 forks source link

BUG: polynomial interpolation on quadrilateral grids #3208

Open cscmaw opened 10 months ago

cscmaw commented 10 months ago

Describe the bug Certain types of polynomial give an error when interpolated on quadrilateral grids.

Steps to Reproduce The following minimal script shows the error:

from firedrake import *
#mesh = UnitSquareMesh(10,10)
mesh = UnitSquareMesh(10,10,quadrilateral=True)
x,y = SpatialCoordinate( mesh )
V = FunctionSpace( mesh, "CG", 1 )
u = Function(V).interpolate(x*x*y*y)
#u = Function(V).interpolate((x*y)**2)

Expected behavior Any of the other options commented out continue without error and work as expected (triangular grids, more compact expression)

Error message

Traceback (most recent call last):
  File "/usr/not-backed-up/firedrake/poly.py", line 9, in <module>
    u = Function(V).interpolate(x*x*y*y)
  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/firedrake/firedrake/src/firedrake/firedrake/function.py", line 410, in interpolate
    return interpolation.interpolate(
  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/firedrake/firedrake/src/firedrake/firedrake/interpolation.py", line 138, in interpolate
    return Interpolator(
  File "/home/firedrake/firedrake/src/pyadjoint/pyadjoint/tape.py", line 109, in wrapper
    return function(*args, **kwargs)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/interpolation.py", line 636, in __init__
    self.callable, arguments = make_interpolator(expr, V, subset, access, bcs=bcs)
  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/firedrake/firedrake/src/firedrake/firedrake/interpolation.py", line 816, in make_interpolator
    loops.extend(_interpolator(V, tensor, expr, subset, arguments, access, bcs=bcs))
  File "<decorator-gen-31>", line 2, in _interpolator
  File "/home/firedrake/firedrake/src/firedrake/firedrake/utils.py", line 75, in wrapper
    return f(*args, **kwargs)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/interpolation.py", line 892, in _interpolator
    kernel = compile_expression(cell_set.comm, expr, to_element, V.ufl_element(),
  File "/home/firedrake/firedrake/src/PyOP2/pyop2/caching.py", line 300, in wrapper
    v = func(*args, **kwargs)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/interpolation.py", line 1026, in compile_expression
    return compile_expression_dual_evaluation(*args, **kwargs)
  File "/home/firedrake/firedrake/src/tsfc/tsfc/driver.py", line 267, in compile_expression_dual_evaluation
    evaluation, basis_indices = to_element.dual_evaluation(fn)
  File "/home/firedrake/firedrake/src/FInAT/finat/finiteelementbase.py", line 270, in dual_evaluation
    evaluation = gem.optimise.contraction(evaluation, shape_indices)
  File "/home/firedrake/firedrake/src/tsfc/gem/optimise.py", line 584, in contraction
    return rebuild(expression)
  File "/home/firedrake/firedrake/src/tsfc/gem/optimise.py", line 556, in rebuild
    return IndexSum(sum_factorise(to_factor, factors), extra)
  File "/home/firedrake/firedrake/src/tsfc/gem/optimise.py", line 358, in sum_factorise
    raise NotImplementedError("Too many indices for sum factorisation!")
NotImplementedError: Too many indices for sum factorisation!

Environment: This is seen in the Docker image on Windows (via Docker Desktop) and linux (via singularity).

wence- commented 10 months ago

This is due to the large hammer in https://github.com/FInAT/FInAT/blob/f493084b14e6a1b68be3bec0c60af4a85931d896/finat/finiteelementbase.py#L270

Something more targetted should be done, but that needs a bunch of work, machinery.

cscmaw commented 10 months ago

Thanks. Will watch out for this in the future.