formoniq: stabilizing the weak Lie derivative, and the top grade that transports nothing

[luiswirth]

The weak Lie derivative (LieDerivativeElmat) is assembled centrally: each cell
integrates its own trace over a shared facet, no numerical flux is chosen, and
the neighbors' disagreement is the whole of the coupling.
That is consistent and, at two grades, exactly conservative.
At the grades between it drifts, and at the top grade it does nothing at all.

This issue records what fails, why, and what the two candidate fixes each buy,
so the choice is made against the diagnosis rather than by analogy to the scalar
case.

What fails

cargo run --release --example advection, transporting a bump 0.2 at unit speed:

dim	grade	drift
2	0	+0.17%
2	1	+14.1%
2	2	0
3	0	+1.0%
3	1	+21.2%
3	2	+5.8%
3	3	0

The pattern is structural in the grade, not incidental.

Why

The antisymmetry defect of the element matrix is exactly a boundary integral,

$$a_K(\omega,\eta) + a_K(\eta,\omega) = \int_{\partial K} \langle\omega,\eta\rangle, \iota_v \mathrm{vol},$$

asserted by the_lie_derivative_is_skew_up_to_its_boundary_term.
It vanishes wherever the facet terms of neighboring cells cancel,
and that happens at both ends of the grade range and nowhere between:

• At $k = 0$ the shape functions are continuous, so $\langle\omega,\eta\rangle$
  is single-valued on a facet and the two sides cancel.
• At $k = n$ they are constant per cell, so $\langle\omega,\eta\rangle$ leaves
  the integral and what remains is $\int_K \operatorname{div} v = 0$.
• In between only the tangential part is single-valued, the normal part jumps,
  and nothing cancels.

Closing the gap in between would need a Killing field,
and a piecewise-flat manifold generally has none:
a Killing field must preserve curvature, hence fix every hinge of nonzero angle
defect, which on a simplicial sphere leaves only the zero field.
So this is not a field that can be chosen better.
The velocity construction in studio already takes the smoothest closed field
the mesh admits, which is as close as the manifold allows, and grade 1 still
drifts.

The top grade is worse than inaccurate

It is the zero operator. Measured on the flat box with a divergence-free field,
the assembled Lie derivative has $|A|_\infty$ of 3.6e-14 in two dimensions and
9.1e-13 in three: the top-grade field does not move at all.

The reason is exact rather than numerical. At top grade $d\omega = 0$ kills the
volume term, and with $\omega$, $\eta$ constant per cell the boundary term is
$f g \oint_{\partial K}\sigma = f g,(d\sigma)_K = 0$ for a divergence-free
field. This is the classical degeneracy of a central flux on piecewise
constants, and it is the strongest argument in this issue: at $k = n$ the scheme
is not imprecise, it is vacuous.

The two candidates

Upwinding. Replace the own-side trace on a facet with one selected by the
flow: the interior trace on outflow, the neighbor's on inflow.
The facet term is then symmetric negative semi-definite, so it damps and never
injects.

It needs machinery that does not exist. assemble_galmat iterates cells and
scatters a square block over a cell's own faces; an upwind term writes into the
neighbor's degrees of freedom, so it needs a facet assembler, and the two
cells' shape functions must be brought into a common frame through
Transition::differential. It also costs the exact conservation at $k = 0$ and
$k = n$, where the central scheme is already right, and it does not fix a
velocity whose divergence is nonzero, which is physics rather than instability.

What it does fix that nothing else does is the top grade, where a non-central
flux is the only thing that makes the operator nonzero at all.

Galerkin least squares. Add

$$\delta \sum_K \int_K \langle \mathcal{L}_v\omega, \mathcal{L}_v\eta\rangle,$$

symmetric positive semi-definite, dissipation aligned with the streamline,
entirely inside the cell. It reduces to classical SUPG at grade 0 and needs no
facet assembler, no neighbor, no Transition — the element matrix machinery
already carries it. But it damps rather than transports, so it is unlikely to
rescue the top grade, where the operator has nothing to stabilize.

Suggested order

1. GLS first, since it needs nothing structurally new, and measure whether it
   closes the intermediate-grade drift.
2. The facet assembler after, for upwinding, whose real target is the top grade.
   Its first validation is free: a facet-assembled central operator must equal
   the current element matrix's boundary half exactly.

Boundary conditions come with step 2 rather than separately. The operator's own
facet terms act on a domain boundary today, which is outflow-like and not a
prescribed inflow, so a solution transported into the boundary is not
meaningful; the inflow case of the upwind switch is where the data belongs.

What already holds and should keep holding

• Conservation at $k = 0$ and $k = n$ on any closed mesh, exactly.
  grade_zero_transport_conserves_on_a_closed_mesh asserts it below 1e-10.
• The antisymmetry identity above, at every grade and dimension.
• The velocity's flux is a cocycle, so the discrete divergence vanishes to the
  digit, and its facet flux is single-valued.

A stabilization that breaks the first is trading away something the central
scheme gets right, and that trade should be deliberate.