Sum factorisation

FFCx enables sum factorisation by default wherever the elements and quadrature rule support it. The sum_factorization code-generation option defaults to True; applications do not need to set it explicitly.

For eligible cell integrals, FFCx evaluates coefficients, the spatial coordinate and the cell Jacobian on the tensor-product quadrature rule, then contracts the argument tables one direction at a time. This reduces the work for a bilinear form of degree \(p\) in dimension \(d\) from \(O(p^{3d})\) to \(O(p^{2d+1})\) when the number of quadrature points per direction grows proportionally to \(p\).

Element and quadrature requirements

Directional contractions require cell integrals on quadrilaterals or hexahedra, with argument elements whose basis tables admit a supported tensor-product representation and a tensor-product quadrature rule. FFCx constructs tensor-product rules for the supported quadrature schemes on these cells.

FFCx preserves the supplied element’s degree-of-freedom (DOF) ordering. Use Basix’s create_tp_element to create an element with tensor-product DOF ordering, then wrap it for UFL:

import basix
import basix.ufl

element = basix.ufl.wrap_element(
    basix.create_tp_element(
        basix.ElementFamily.P,
        basix.CellType.hexahedron,
        degree=4,
        lagrange_variant=basix.LagrangeVariant.gll_warped,
    )
)

Use basix.CellType.quadrilateral for a quad mesh. Vector and tensor spaces can be built from this scalar element using basix.ufl.blocked_element.

In DOLFINx, create the function space using this element and compile the form normally. Given a hexahedral mesh mesh:

import ufl
from dolfinx import fem

V = fem.functionspace(mesh, element)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
compiled_a = fem.form(a)

With the direct FFCx JIT interface, the same default applies:

from ffcx.codegeneration.jit import compile_forms

compiled_forms, module, code = compile_forms([a])

Automatic fallback

FFCx uses ordinary quadrature when the integral cannot use directional contractions. This includes other cell types, facet integrals, argument tables without supported tensor factors, assembly of only the diagonal (part="diagonal"), and custom or vertex quadrature rules without tensor-product information. These cases do not require disabling the option. A coefficient or geometry table without tensor factors can still be evaluated normally inside a supported contraction. FFCx also uses ordinary quadrature when a kernel’s arrays would not fit the stack limit described below.

Stack usage

A sum-factorised kernel evaluates the integrand, coefficients and geometry at every point of the quadrature rule and keeps these values in arrays on the stack. The space they take grows with the number of quadrature points, so cubically with the degree on hexahedra, and with the number of terms in the integrand. For example, an elasticity form on hexahedra takes about 126 KiB of stack per kernel call at degree 4 and about 510 KiB at degree 8.

The main thread of a program usually has 8 MiB of stack, but other threads can have much less: 512 KiB for secondary threads on macOS, and 128 KiB for threads on musl-based systems such as Alpine Linux. A kernel that exceeds its thread’s stack crashes the program.

The sum_factorization_stack_limit option bounds this space, in bytes. It defaults to 256 KiB. When a kernel’s arrays would exceed it, FFCx splits the terms of the integrand into batches and evaluates and contracts one batch at a time, so that each batch’s arrays fit. This costs some speed: the integrand is evaluated once per batch, and the element tensor is updated once per batch. If even a single term does not fit, for example with very large quadrature rules, FFCx uses ordinary quadrature for that integral. Set the option to 0 for no limit.

With the FFCx JIT:

compiled_forms, module, code = compile_forms(
    [a], options={"sum_factorization_stack_limit": 96 * 1024}
)

With DOLFINx:

compiled_a = fem.form(
    a, form_compiler_options={"sum_factorization_stack_limit": 96 * 1024}
)

From the command line:

ffcx --sum_factorization_stack_limit 98304 forms.py

Lower the limit if kernels run on threads with small stacks: about 96 KiB is safe for 128 KiB threads. Raise it, or set it to 0, if kernels only run on threads with large stacks, such as the main thread, and the fastest kernels at high degree matter more.

Disabling sum factorisation

To use ordinary quadrature explicitly with the FFCx JIT:

compiled_forms, module, code = compile_forms(
    [a], options={"sum_factorization": False}
)

With DOLFINx, set the corresponding form compiler option:

compiled_a = fem.form(
    a, form_compiler_options={"sum_factorization": False}
)

From the command line:

ffcx --no-sum_factorization forms.py

This option controls kernel generation. C compiler flags are configured separately; see Compiler flags for suggested optimization and FMA settings.