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

[garth-wells]

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:

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.

[jorgensd]

This is not something there is any point of pursuing in my opinion, as we can and should get the right forms if we use ufl.MixedFunctionSpace and our blocked assemblers.

[garth-wells]

Fix posted: FEniCS/ufl#503

Root cause was in SumDegreeEstimator.argument()/.coefficient() (ufl/algorithms/estimate_degrees.py) -- both return the whole mixed element's embedded_superdegree for any component, rather than the specific sub-element's degree, exactly as the existing # FIXME: Use component to improve accuracy for mixed elements comment flagged.

Fixed in indexed(): when the indexed operand is an Argument/Coefficient and the multi-index is fully fixed (not a free index), flatten it into the operand's component space, find which sub-element that component belongs to, and return that sub-element's own degree. Falls back to the original whole-element degree otherwise.

Verified:

• stokes_th_tet's estimated degree drops from 3 to 2, giving the minimal 4-point rule instead of 5 points (with a negative weight) -- confirmed numerically identical (max abs diff 1.5e-15) since the true degree was 2 all along.
• Measured 1.23x faster on the affected kernel (921ns -> 749ns, -O2), with every other kernel in this survey unaffected.
• Full UFL test suite (973 tests) and FFCx's own suite (212 tests) both pass -- two pre-existing UFL test assertions had # TODO: This should be N, but M is expected behaviour now comments documenting this exact bug; both now assert the corrected value.

PR: FEniCS/ufl#503