[BUG]: Assembling sum of quadrature form on volume and CG form on surface

[SuperFranck]

How to reproduce the bug

Hi everyone,

This code aims to highlights a bug (in my humble opinion) in fenicsx-0.6.0.

We use quadrature and CG elements.
We write an equation mixing these elements.
When trying to form this equation, an error appears. However, using fenicsx-0.5.1 there is no problem.

Note also that the way we write the equation have an impact on whether it works.
We propose three writings that work fine, and one other causing trouble.

Error message is not short, but I hope it may help you

Minimal Example (Python)

import dolfinx as df
import ufl
import numpy as np
import mpi4py

msh = df.mesh.create_rectangle(mpi4py.MPI.COMM_WORLD, [np.array([0,0]), np.array([1, 0.1])], [1,1], cell_type=df.mesh.CellType.triangle) 

CG2_vect = df.fem.FunctionSpace(msh, ("CG", 1))
Qe = ufl.FiniteElement("Quadrature", msh.ufl_cell(), degree=1, quad_scheme='default')
Quad = df.fem.FunctionSpace(msh, Qe)

u = df.fem.Function(Quad)
v = ufl.TrialFunction(CG2_vect)

dx_m = ufl.Measure("dx",domain=msh, metadata={"quadrature_degree": 1, "quadrature_scheme": "default"}) 
ds = ufl.Measure("ds", domain=msh)

residual = u * v * dx_m 
df.fem.form(residual)
residual = v * ds
df.fem.form(residual)
residual = u * v * dx_m - v * dx_m
df.fem.form(residual)
print("ok")


residual = u * v * dx_m - v * ds
df.fem.form(residual)

Output (Python)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[21], line 39
     35 print("ok")
     38 residual = u * v * dx_m - v * ds
---> 39 df.fem.form(residual)

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/dolfinx/fem/forms.py:176, in form(form, dtype, form_compiler_options, jit_options)
    173         return list(map(lambda sub_form: _create_form(sub_form), form))
    174     return form
--> 176 return _create_form(form)

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/dolfinx/fem/forms.py:171, in form.._create_form(form)
    168 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    169 return form argument"""
    170 if isinstance(form, ufl.Form):
--> 171     return _form(form)
    172 elif isinstance(form, collections.abc.Iterable):
    173     return list(map(lambda sub_form: _create_form(sub_form), form))

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/dolfinx/fem/forms.py:145, in form.._form(form)
    142 if mesh is None:
    143     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 145 ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
    146                                        form_compiler_options=form_compiler_options,
    147                                        jit_options=jit_options)
    149 # For each argument in form extract its function space
    150 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/dolfinx/jit.py:56, in mpi_jit_decorator..mpi_jit(comm, *args, **kwargs)
     51 @functools.wraps(local_jit)
     52 def mpi_jit(comm, *args, **kwargs):
     53 
     54     # Just call JIT compiler when running in serial
     55     if comm.size == 1:
---> 56         return local_jit(*args, **kwargs)
     58     # Default status (0 == ok, 1 == fail)
     59     status = 0

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/dolfinx/jit.py:204, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    202 # Switch on type and compile, returning cffi object
    203 if isinstance(ufl_object, ufl.Form):
--> 204     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    205 elif isinstance(ufl_object, ufl.FiniteElementBase):
    206     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:186, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    183     for name in form_names:
    184         decl += form_template.format(name=name)
--> 186     impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
    187                             cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    188 except Exception:
    189     # remove c file so that it will not timeout next time
    190     c_filename = cache_dir.joinpath(module_name + ".c")

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:250, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries)
    246 import ffcx.compiler
    248 # JIT uses module_name as prefix, which is needed to make names of all struct/function
    249 # unique across modules
--> 250 _, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
    252 ffibuilder = cffi.FFI()
    253 ffibuilder.set_source(module_name, code_body, include_dirs=[ffcx.codegeneration.get_include_path()],
    254                       extra_compile_args=cffi_extra_compile_args, libraries=cffi_libraries)

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/compiler.py:102, in compile_ufl_objects(ufl_objects, object_names, prefix, options, visualise)
    100 # Stage 2: intermediate representation
    101 cpu_time = time()
--> 102 ir = compute_ir(analysis, object_names, prefix, options, visualise)
    103 _print_timing(2, time() - cpu_time)
    105 # Stage 3: code generation

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/representation.py:200, in compute_ir(analysis, object_names, prefix, options, visualise)
    190 ir_elements = [
    191     _compute_element_ir(e, analysis.element_numbers, finite_element_names)
    192     for e in analysis.unique_elements
    193 ]
    195 ir_dofmaps = [
    196     _compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
    197     for e in analysis.unique_elements
    198 ]
--> 200 irs = [
    201     _compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
    202                          options, visualise)
    203     for (i, fd) in enumerate(analysis.form_data)
    204 ]
    205 ir_integrals = list(itertools.chain(*irs))
    207 ir_forms = [
    208     _compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers, finite_element_names,
    209                      dofmap_names, object_names)
    210     for (i, fd) in enumerate(analysis.form_data)
    211 ]

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/representation.py:201, in (.0)
    190 ir_elements = [
    191     _compute_element_ir(e, analysis.element_numbers, finite_element_names)
    192     for e in analysis.unique_elements
    193 ]
    195 ir_dofmaps = [
    196     _compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
    197     for e in analysis.unique_elements
    198 ]
    200 irs = [
--> 201     _compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
    202                          options, visualise)
    203     for (i, fd) in enumerate(analysis.form_data)
    204 ]
    205 ir_integrals = list(itertools.chain(*irs))
    207 ir_forms = [
    208     _compute_form_ir(fd, i, prefix, form_names, integral_names, analysis.element_numbers, finite_element_names,
    209                      dofmap_names, object_names)
    210     for (i, fd) in enumerate(analysis.form_data)
    211 ]

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/representation.py:485, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise)
    482 integrands = {rule: integral.integrand() for rule, integral in sorted_integrals.items()}
    484 # Build more specific intermediate representation
--> 485 integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
    486                                   ir["entitytype"], integrands, ir["tensor_shape"],
    487                                   options, visualise)
    489 ir.update(integral_ir)
    491 # Fetch name

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/integral.py:84, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, p, visualise)
     74 # Build terminal_data from V here before factorization. Then we
     75 # can use it to derive table properties for all modified
     76 # terminals, and then use that to rebuild the scalar graph more
     77 # efficiently before argument factorization. We can build
     78 # terminal_data again after factorization if that's necessary.
     80 initial_terminals = {i: analyse_modified_terminal(v['expression'])
     81                      for i, v in S.nodes.items()
     82                      if is_modified_terminal(v['expression'])}
---> 84 mt_table_reference = build_optimized_tables(
     85     quadrature_rule,
     86     cell,
     87     integral_type,
     88     entitytype,
     89     initial_terminals.values(),
     90     ir["unique_tables"],
     91     rtol=p["table_rtol"],
     92     atol=p["table_atol"])
     94 # Fetch unique tables for this quadrature rule
     95 table_types = {v.name: v.ttype for v in mt_table_reference.values()}

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/elementtables.py:361, in build_optimized_tables(quadrature_rule, cell, integral_type, entitytype, modified_terminals, existing_tables, rtol, atol)
    359             t['array'] = numpy.vstack([td['array'] for td in new_table])
    360 else:
--> 361     t = get_ffcx_table_values(quadrature_rule.points, cell,
    362                               integral_type, element, avg, entitytype,
    363                               local_derivatives, flat_component)
    364 # Clean up table
    365 tbl = clamp_table_small_numbers(t['array'], rtol=rtol, atol=atol)

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/elementtables.py:123, in get_ffcx_table_values(points, cell, integral_type, element, avg, entitytype, derivative_counts, flat_component)
    120 component_element, offset, stride = element.get_component_element(flat_component)
    122 for entity in range(num_entities):
--> 123     entity_points = map_integral_points(points, integral_type, cell, entity)
    124     tbl = component_element.tabulate(deriv_order, entity_points)
    125     tbl = tbl[basix_index(derivative_counts)]

File ~/anaconda3/envs/fenicsx-0.6.0/lib/python3.10/site-packages/ffcx/ir/representationutils.py:92, in map_integral_points(points, integral_type, cell, entity)
     90     return numpy.asarray(points)
     91 elif entity_dim == tdim - 1:
---> 92     assert points.shape[1] == tdim - 1
     93     return numpy.asarray(map_facet_points(points, entity, cell.cellname()))
     94 elif entity_dim == 0:

Version

0.6.0

DOLFINx git commit

No response

Installation

I used an environment file and run conda env create --file fenicsx-0.6.0.yml

Additional information

No response

[jorgensd]

This is indeed a bug, caused by the lines: https://github.com/FEniCS/ffcx/blob/main/ffcx/analysis.py#L174-L183
which shouldn't be on the form level, but on the integral level.

This would require some refactoring of the code (the analysis data), @mscroggs any idea for a quick fix?

[jhale]

Could we change this code so that it simply verifies that the measure quadrature rule is the same as any quadrature element used within the integral, and if not, throw an error?

[jorgensd]

I dont think that is a good solution, as the form presented by @SuperFranck seems valid to me.

[jhale]

You are right! This is not a fix. But that wasn't quite my point (which I did not put across well).

The current code is trying to be 'too clever' by reconciling integration rules in QuadratureElements with rules in the measure. I would redesign the code to simply assert that there is a good match, then the logic here can be much simpler. It is the user who should write the correct equations, after all.

[a-latyshev]

Dear community,

I can suggest you a temporary solution of this issue.

Actually, if you use variables from vectoral quadrature space only in your forms, your code will work correctly.

In order to deal with scalars I suggest following modifications to the minimal example of the issues described above.

Q2e = ufl.VectorElement("Quadrature", msh.ufl_cell(), degree=1, dim=2, quad_scheme='default')
Quad2 = df.fem.FunctionSpace(msh, Q2e)
u2 = df.fem.Function(Quad2)
new_u = u2[0]
new_u_values = u2.x.array.reshape((-1, 2))[:,0] # For data update
residual = new_u * v * dx_m - v * ds
df.fem.form(residual)

P.S. Tested for fenicsx-0.7.0.

Cheers!

[mscroggs]

I'm hoping this may be fixed by FEniCS/ffcx#622 (I can at least run the minimal example without a failure now)