FEniCS / FEniCS/ffcx

Quadrature degree over-estimated for mixed-element form components (stokes_th_tet: 5-point rule instead of the minimal 4-point rule)

Open
#866 2 comments 0 reactions 0 assignees View on GitHub
Dominant language
Python
Stars
192
Forks
45
Avg merge
1d 11h
Merged PRs (30d)
16

Description

## Mixed-element degree estimation over-counts quadrature degree for Taylor-Hood-style forms (stokes_th_tet: 5-point rule instead of the minimal 4-point rule)

Investigating `stokes_th_tet`'s quadrature cost as part of ongoing performance work on `garth/perf-stack`.

**Root cause, confirmed empirically**: FFCx's generated kernel for `(inner(grad(u), grad(v)) - div(v)*p + div(u)*q) * dx` (Taylor-Hood Stokes, P2 vector / P1 scalar mixed space) uses a 5-point degree-3 quadrature rule (`weights_175[5] = {-0.133..., 0.075, 0.075, 0.075, 0.075}`, note the negative centroid weight). The true minimal exact degree for this form is 2, which only needs 4 points with strictly positive weights — confirmed by checking each term's UFL-estimated degree separately:

```
inner(grad(u),grad(v)) degree=2 (correct)
div(v)*p degree=3 (should be 2)
div(u)*q degree=3 (should be 2)
```

`div(v)` and `div(u)` correctly estimate to degree 1. The problem is `p`/`q` (the *scalar*, degree-1 sub-arguments extracted via `TrialFunctions`/`TestFunctions` on a `MixedFunctionSpace`) themselves estimate to **degree 2** — the mixed space's overall max degree, not the P1 subelement's actual degree of 1.

**Located precisely** in `ufl/algorithms/estimate_degrees.py`'s `SumDegreeEstimator`:

```python
def argument(self, v):
return (
v.ufl_element().embedded_superdegree
) # FIXME: Use component to improve accuracy for mixed elements

def coefficient(self, v):
...
d = e.embedded_superdegree # FIXME: Use component to improve accuracy for mixed elements
```

Both handlers operate on the whole `Argument`/`Coefficient` node's element (the full `MixedElement`) — the FIXME comment is already there, written by UFL's own maintainers, acknowledging this exact limitation. Compounding it: `indexed(self, v, A, ii)` just returns `A` unchanged — the component index `ii` that would let a fix identify which sub-element is actually being accessed is discarded before it would ever reach `argument()`/`coefficient()` (which fire on the base node, not the `Indexed` wrapper).

**Impact**: this isn't `stokes_th`-specific — any mixed-space form combining subelements of different polynomial degree (any Taylor-Hood-like Stokes/Navier-Stokes discretization, mixed Poisson, poroelasticity, multi-field problems generally) will over-estimate quadrature degree by however much the mixed space's max-degree subelement exceeds the specific component's actual degree, whenever that component is used un-differentiated (differentiation already correctly reduces degree via `_reduce_degree`, which is why `div(v)`/`div(u)` were fine — it's the *undifferentiated* `p`/`q` that leak the wrong degree). A 5-vs-4-point difference is a 25% cut in quadrature-point count, and hence proportionally less per-cell work, for any affected form.

**Fix scope**: in UFL (`ufl/algorithms/estimate_degrees.py`), not FFCx-specific — benefits every downstream form compiler using UFL. Needs the `indexed()`/`component_tensor()` handlers to track which flattened component is being accessed and thread that through to `argument()`/`coefficient()` so they can look up the correct sub-element's `embedded_superdegree` via the mixed element's component-to-subelement map (already available, e.g. via `sub_elements()` and basix.ufl's own component bookkeeping) instead of the whole mixed element's.

Not yet attempted — flagging as a fresh, concrete, high-value lead distinct from #864.

Contributor guide

No contributing guide indexed for this repository

Research direction

Start in ufl/algorithms/estimate_degrees.py, reading SumDegreeEstimator’s argument, coefficient, indexed, and component_tensor handlers. Reproduce the stokes_th_tet form and compare the estimated degrees of the mixed components. Done means undifferentiated mixed subcomponents use their own degree and the affected form no longer selects an unnecessarily higher-degree quadrature rule.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
compilers, performance
Issue type
Bug
Difficulty
4/5
Estimated time
3-5 days
Activity status
Active
Clarity
Mostly clear
Newbie friendliness
48/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.