Quadrature degree over-estimated for mixed-element form components (stokes_th_tet: 5-point rule instead of the minimal 4-point rule)
- 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