Investigate whether the matrix-free operator algebra extends to discontinuous Galerkin

[kylebeggs]

From Dennis Ogiermann, 2026-08-14: "This can also be interesting for discontinuous Galerkin methods, but they have some pain-points (surface integration)."

Filing this for the record with the lowest priority of the five suggestions in that message. The honest default answer is "not in this package", and the point of the issue is to establish that with an argument rather than a shrug — or to find out it is wrong.

The blocker, stated precisely

The package's solver boundary is a flat vector over interior cell DOFs. size(L) spans interior cell DOFs (DESIGN.md:312) and mul! is an interior-only flat-vector map (DESIGN.md:1037-1044), with CLAUDE.md restating it as an invariant: "Flat (Krylov) vectors span interior DOFs only." DG carries Np DOFs per element — modal or nodal — with no continuity between elements. Those two statements are not compatible as written, and every prepared-operator path, flatten/flat_to_interior!, and operator_diagonal sits on the cell-DOF side of it.

There is also no discretization abstraction to hang an alternative on. The complete list of abstract types is AbstractGrid, AbstractField, AbstractBlockField, BlockLayout, AbstractOperator, AbstractBC, AbstractJVPBackend — no AbstractDiscretization, no scheme trait, no basis, quadrature, or element concept. Field{Center} is the only location. Several operators hard-guard on the concrete grid type anyway (src/operators/diffusion.jl:234, src/operators/diagonal.jl:42, src/operators/diagonal.jl:126, src/multigrid.jl:367, src/distributed.jl:57).

The one angle that is genuinely interesting

The package's "element type carries tensor rank" invariant means a Field is an array of SVectors and operators have no rank parameter. In principle the same trick carries per-element DG DOFs: a Field over Array{SVector{Np,T}}, one static vector of modal coefficients per element, no new rank machinery. That is a real structural affordance and it is why the question is not absurd.

What it collides with immediately:

• Gradient / Divergence are the declared rank-changers. Under a DG reading an SVector means modal coefficients, not vector components, and those two operators would mean something entirely different on the same type. Nothing distinguishes the two interpretations today.
• _scalar_eltype (src/Fields.jl:134-135) unwraps SVector to get the Krylov scalar type. It would have to tell "3 components of a vector field" apart from "6 modal coefficients of a P2 element."
• operator_diagonal and the multigrid smoothers assume one scalar per cell.

The surface-integration pain point Dennis named

It is worse here than in a generic DG code, and for a specific reason: halo_update! ships cell values. DG face terms need face traces u⁻, u⁺ and a numerical flux evaluated on the face. So DG needs a face-trace exchange, which is a third distinct exchange shape alongside the two already in flight — the corner-aware exchange a cross-derivative tensor leaf needs, and the "face-flux descriptor" seam #54 raises as an open question against its FluxDivergence. Whether those three are one mechanism or three is the most useful thing this investigation could actually determine.

Proposed investigation

1. Write down what a DG residual application needs from a grid, a field, and an exchange, in the vocabulary the package already has.
2. Decide whether the flat-vector contract can be widened to Np DOFs per cell without breaking prepare, boundary_rhs, the declared adjoints, or the multigrid path — or whether widening it is exactly the kind of speculative abstraction DESIGN.md:862-870 forbids ("deferred until a second … actually exists … Do not introduce it speculatively").
3. Check whether feat(fv): add conservative finite-volume operators without sacrificing matrix-free AMR #54's numerical-flux-callable FluxDivergence subsumes the face machinery, in which case the answer changes from "different package" to "a leaf on top of the FV seam."
4. Recommend. A clear "no, and here is the boundary" is a successful outcome and should be written into DESIGN.md §10 so it stops being re-asked.

Out of scope

Implementing anything. No basis/quadrature layer, no mortar or hp machinery, no limiters, and specifically no AbstractDiscretization introduced on the strength of one hypothetical use case.

---

_{🤖 Claude investigated whether the cells could hold more than one number each, and mostly found reasons they shouldn't.}

[kylebeggs]

Cross-reference: step 3 of the proposed investigation ("does #54's numerical-flux-callable FluxDivergence subsume the face machinery?") depends on an open decision in #54, so this issue should not start before that one answers whether FluxDivergence becomes the public primitive. The diffusion slice that would have been the first flux to generalize from (#56) is now complete, so #54 can take that question up.

---

_{🤖 Beep boop — Claude linked the surface integrals; Kyle is not integrating himself by parts.}