1. FEniCS Project Avatar FEniCS Project
  2. Core
  3. FFC
  4. Issues

Crashes when projecting scalar to quadrature space of degree 2.

Issue #182 new
Grant Bruer created an issue

This is related to an issue I raised in the DOLFIN repo, but the bug is actually in the FFC part so I’m making the issue here. (https://bitbucket.org/fenics-project/dolfin/issues/1084/crashes-when-projecting-to-quadrature)

Quadrature representation

The code here tries to project a function onto a quadrature space using the quadrature representation:

from fenics import *

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'P', 1)

u = Function(V)

# Create quadrature space.
quad_params = {
    'quadrature_degree': 2,        # 1 works here, but not 2 or higher.
    'quadrature_rule': 'default',
    'representation': 'quadrature' # Also tried uflacs here as warning suggested.
}
quad_el = TensorElement('Quadrature', mesh.ufl_cell(), shape=(),
    degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])
Q = FunctionSpace(mesh, quad_el)

qdx = dx(metadata=quad_params)

# Set up a system for projecting u onto Q.
v = TestFunction(Q)
q = TrialFunction(Q)

a = q * v * qdx
L = u * v * qdx

# Solve the system.
q_ = Function(Q)
solve(a == L, q_)

It fails with the following error, due to a wrong shape for some table.

/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturerepresentation.py:58: QuadratureRepresentationDeprecationWarning: 
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** representation instead.                               ***
*** ===================================================== ***
  issue_deprecation_warning()
Traceback (most recent call last):
  File "mwe.py", line 29, in <module>
    solve(a == L, q_)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 183, in <listcomp>
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/representation.py", line 460, in _compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturerepresentation.py", line 100, in compute_integral_ir
    ir["optimise_parameters"])
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/optimisedquadraturetransformer.py", line 56, in __init__
    QuadratureTransformerBase.__init__(self, *args)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadraturetransformerbase.py", line 148, in __init__
    self.entity_type)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadratureutils.py", line 49, in create_psi_tables
    name_map, unique_tables = unique_psi_tables(flat_tables, eliminate_zeros)
  File "/home/gbruer/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ffc/quadrature/quadratureutils.py", line 179, in unique_psi_tables
    if abs(vals[r][c]) < format_epsilon:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I narrowed down the bug to the flatten_psi_tables() function in quadratureutils.py , line 118. The function assumes that having a value_shape() == () means that the table is of the form fiat_tables[ip][dof], but for a scalar function defined with a TensorElement, the table will still be in the form fiat_tables[ip][component][dof], with a shape of (#ip, 1, #dof).

                    # Transform fiat_tables to a list of tables on
                    # the form psi_table[dof][ip] for each scalar
                    # component
                    if element.value_shape():
                        # Table is on the form
                        # fiat_tables[ip][component][dof].
                        transposed_table = numpy.transpose(fiat_tables, (1, 2, 0))
                        component_tables = list(enumerate(transposed_table))
                    else:
                        # Scalar element, table is on the form
                        # fiat_tables[ip][dof].  Using () for the
                        # component because generate_psi_name
                        # expects that
                        component_tables = [((),
                                             numpy.transpose(fiat_tables))]

Proposed Fix

A possible fix is to modify this section to take into account that the value_shape does not directly determine the shape of the table.

                        # Transform fiat_tables to a list of tables on
                        # the form psi_table[dof][ip] for each scalar
                        # component
                        if len(numpy.shape(fiat_tables)) == 3:
                            # Table is on the form
                            # fiat_tables[ip][component][dof].
                            transposed_table = numpy.transpose(fiat_tables, (1, 2, 0))
                            if element.value_shape():
                                component_tables = list(enumerate(transposed_table))
                            else:
                                # Need to use () for the component because
                                # generate_psi_name expects that. Note that the
                                # shape of transposed_table is 1 x #dof x #ip.
                                component_tables = [((), transposed_table[0])]
                        else:
                            # Scalar element, table is on the form
                            # fiat_tables[ip][dof].  Using () for the
                            # component because generate_psi_name
                            # expects that
                            component_tables = [((),
                                                 numpy.transpose(fiat_tables))]

Workaround

A workaround for those encountering this issue is to not use a scalar TensorElement.

quad_el = TensorElement('Quadrature', mesh.ufl_cell(), shape=(),
    degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])

Instead, use a scalar FiniteElement.

quad_el = FiniteElement('Quadrature', mesh.ufl_cell(),
    degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])

UFLACS representation

The UFLACS representation also fails, whether or not a TensorElement or FiniteElement is used.

from fenics import *

mesh = UnitSquareMesh(4, 4)
V = FunctionSpace(mesh, 'P', 1)

u = Function(V)

# Create quadrature space.
quad_params = {
    'quadrature_degree': 2,        # 1 works here, but not 2 or higher.
    'quadrature_rule': 'default',
    'representation': 'uflacs' # Also tried uflacs here as warning suggested.
}
quad_el = FiniteElement('Quadrature', mesh.ufl_cell(),
    degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])
# quad_el = TensorElement('Quadrature', mesh.ufl_cell(), shape=(),
#     degree=quad_params['quadrature_degree'], quad_scheme=quad_params['quadrature_rule'])

Q = FunctionSpace(mesh, quad_el)

qdx = dx(metadata=quad_params)

# Set up a system for projecting u onto Q.
v = TestFunction(Q)
q = TrialFunction(Q)

a = q * v * qdx
L = u * v * qdx

# Solve the system.
q_ = Function(Q)
solve(a == L, q_)

TensorElement

If the scalar TensorElement is used, the following error occurs:

Calling FFC just-in-time (JIT) compiler, this may take some time.
Traceback (most recent call last):
  File "uflacs.py", line 29, in <module>
    solve(a == L, q_)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/jitcompiler.py", line 66, in jit_generate
    prefix=module_name, parameters=parameters, jit=True)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/compiler.py", line 143, in compile_form
    prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/compiler.py", line 190, in compile_ufl_objects
    ir = compute_ir(analysis, prefix, parameters, jit)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/representation.py", line 183, in compute_ir
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/representation.py", line 183, in <listcomp>
    for (form_id, fd) in enumerate(form_datas)]
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/representation.py", line 460, in _compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/uflacsrepresentation.py", line 139, in compute_integral_ir
    parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/build_uflacs_ir.py", line 353, in build_uflacs_ir
    rtol=p["table_rtol"], atol=p["table_atol"])
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 553, in build_optimized_tables
    modified_terminals, rtol=rtol, atol=atol)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 385, in build_element_tables
    name_ignored = add_table(subres)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 362, in add_table
    entitytype, local_derivatives, flat_component)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/uflacs/elementtables.py", line 217, in get_ffc_table_values
    num_dofs, num_points = component_tables[0].shape
ValueError: too many values to unpack (expected 2)

FiniteElement

If a quadrature FiniteElement is used, the following error occurs:

Calling FFC just-in-time (JIT) compiler, this may take some time.
------------------- Start compiler output ------------------------
/tmp/tmp5i2gtxyb/ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58.cpp: In member function 'virtual void ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int) const':
/tmp/tmp5i2gtxyb/ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58.cpp:113:26: error: 'FE5_C0_Q3' was not declared in this scope
         BF0[iq] += fw0 * FE5_C0_Q3[0][iq][iq];
                          ^~~~~~~~~
/tmp/tmp5i2gtxyb/ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58.cpp:113:26: note: suggested alternative: 'FE3_C0_Q3'
         BF0[iq] += fw0 * FE5_C0_Q3[0][iq][iq];
                          ^~~~~~~~~
                          FE3_C0_Q3

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/gbruer/research/crikit-demo/jitfailure-ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58
Traceback (most recent call last):
  File "uflacs.py", line 31, in <module>
    solve(a == L, q_)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/solving.py", line 242, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/problem.py", line 55, in __init__
    L = Form(L, form_compiler_parameters=form_compiler_parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/fem/form.py", line 44, in __init__
    mpi_comm=mesh.mpi_comm())
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dolfin/jit/jit.py", line 97, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/ffc/jitcompiler.py", line 133, in jit_build
    generate=jit_generate)
  File "/home/gbruer/anaconda3/envs/gitlab-ci/lib/python3.7/site-packages/dijitso/jit.py", line 217, in jit
    % err_info['fail_dir'], err_info)
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see '/home/gbruer/research/crikit-demo/jitfailure-ffc_form_ff53cc8b62e5a4496f847e87b12359f54e3f0b58' for details

‌

Comments (1)

  2. Log in to comment
Assignee
–
Type
bug
Priority
major
Status
new
Component
–
Milestone
–
Version
–
Votes
0
Watchers
1
Jira: the preferred issue tracker for Bitbucket. Join the team!