Boundary flux on degree-3 traces: per-slot DOF identity, true edge-node coordinates, and the consistent line mass (#459)

[lmoresi]

2D boundary-flux recovery: degree-3 traces returned silently wrong flux — key by interpolation node, not by DMPlex point

Fixes#459

Mechanism

A degree-3 trace carries TWO interpolation nodes per boundary edge point. The
CBF node gather (_boundary_field_nodes) emitted one entry per DMPlex point,
with the point's single coordinate (the edge midpoint from the coordinate
section). Three defects followed, none of them raising:

1. boundary_flux read only the first DOF at each point — one of the two
   edge-interior reactions was never read at all;
2. the coordinate keys for the two nodes collapse onto one dictionary entry in
   the de-smear's nodeR, so even per-DOF reads would silently overwrite;
3. the 2D trace-order detection sniffed "midpoint key present → P2" and
   de-smeared the P3 reactions with the P2 line mass.

Measured (8×8 StructuredQuadBox, T = 1 − y, exact unit flux at every degree):
the degree-3 wall read flux 0.565–0.739 against an exact 1.0, from 17 of the
25 trace nodes — matching the numbers in the issue. Degrees 1–2 were correct.

Fix (src/underworld3/utilities/boundary_flux.py)

Key by interpolation node, never by point:

• _boundary_field_nodes now emits (point, slot, coord) per NODE. True
  edge-interior node coordinates come from a new _trace_interior_coords:
  the mesh coordinate field interpolated into a matching-degree space, built
  exactly as UW3 builds every field FE (createDefault,
  node_endpoints=False, the Mesh._get_coords_for_basis construction), so
  the per-point slot ordering matches the field section by construction.
  Verified directly: the constrained DOF values of T = xy on the wall equal
  the claimed coordinates to machine precision. The build is collective;
  whether it is needed is agreed by allreduce so ranks owning none of the
  boundary participate (np2/np4 tested).
• boundary_flux reads the reaction per (point, slot); vector reactions
  read node-major. _node_normals keys per node.
• _desmear (2D) derives the trace order structurally from the field section
  — no more coordinate-key sniffing — and assembles a general Lagrange line
  mass from the measured node parameters (_line_mass_1d, Gauss–Legendre,
  exact). P1/P2 keep their exact hardcoded matrices, byte-for-byte: the
  degree-1/2 fluxes are bit-identical before/after.
• The de-smear REFUSES a key collapse exactly (len(nodeR) != len(xs), both
  the 2D and 3D branches, shared _node_reactions helper) — the cheap exact
  detection the issue proposed, retained as a permanent guard.

mass="auto" at degree ≥ 3 selects the consistent solve

PETSc places the cubic edge-interior nodes asymmetrically (Gauss–Jacobi,
t ≈ 0.330/0.670, measured — not equispaced). There, row-sum lumping is only
O(h) pointwise for a varying flux: with the fix but lumped mass, a linear flux
(T = xy) shows a flat 1.6e-2 error at every interior node (res 8), scaling
O(h). The consistent line mass is exact away from the (documented) corner
mixing. P1/P2 keep lumped, which is exact for linear flux on their symmetric
layouts — this is why the defect class never appeared at low degree. The
solver docstring (petsc_generic_snes_solvers.pyx) records the new auto
rule.

Measured flux error, before → after (unit-flux oracle, res 8, mass="auto")

trace degree	before	after	nodes reporting
1	8.0e-5	8.0e-5 (bit-identical)	9/9
2	1.2e-4	1.2e-4 (bit-identical)	17/17
3	4.35e-1	3.1e-4	17/25 → 25/25

Placement/ordering oracle (T = xy, flux = x on Top, res 16, mid-wall):
5.8e-4 against ~2.1e-2 for a hypothetical slot swap (36× separation).

Pattern audit (the issue class, not just the instance)

• 3D branch: one node per point holds for its supported P1/P2 traces; 3D
  P3 was and remains rejected legibly (facet/edge DOF checks). The Free surface in 3D: spherical shell support, dimension-general surface machinery, and the 2D compatibility-floor fix #469
  mass="p1" fold is untouched. The collapse guard now also protects the 3D
  nodeR.
• rotated_bc.boundary_normal_traction builds its own per-point node list
  (P2 velocity: one node per edge — sound). A degree-3 velocity would need
  the same per-node treatment there; with this PR the shared _desmear now
  detects that case structurally and refuses with a legible error ("trace
  nodes have no reaction entry") instead of computing garbage. Filed as
  follow-up work rather than changed here — the rotated-BC reaction path has
  its own conventions and P3 velocity is not a current use case.
• write_boundary_scalar_field / boundary_flux_field: coordinate-keyed
  hand-off is per-node and collision-free after the fix; a P1 target field
  consumes vertex keys only (vertex recovery is exact under both masses).
• 2D corner semantics (pre-existing, unchanged): with the consistent mass,
  the documented corner mixing decays over ~one element instead of staying
  strictly at the corner nodes as with lumping. Visible at degree 2 with
  mass="consistent" on all-wall-driven problems on development today; not
  a regression from this PR.

Tests

• tests/test_1019_boundary_flux.py
  • test_boundary_flux_degree_sweep_2d[1|2|3] — unit-flux oracle, every
    node exact, every node reporting. Fail-before validated on the unfixed
    build: degree 3 failed, degrees 1–2 passed.
  • test_boundary_flux_p3_interior_node_placement — linear-flux oracle;
    catches slot swaps/misplacement the constant-flux oracle cannot see.
  • test_boundary_flux_p3_collapse_guard — negative control: monkeypatch
    the coordinate build back to the pre-fix collapse; the guard must raise.
• tests/parallel/test_1065_boundary_flux_parallel.py
  • test_boundary_flux_degree3_partition_independent — np2 and np4 (the
    per-node coordinate build is collective; empty ranks must participate).

Verification

• tests/test_1019_boundary_flux.py: 14 passed (serial).
• tests/parallel/test_1065_boundary_flux_parallel.py: 4 passed at np2, 4 at np4.
• tests/test_1018_rotated_freeslip.py (σ_nn reaction consumers): 19 passed.
• Full gate pytest tests -m "level_1 and tier_a" --ignore=tests/test_0050_utils.py: 601 passed, 0 failed (17 skipped, 1 xfailed).

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes incorrect (but silently finite) 2D boundary-flux recovery for degree-3 traces by switching from “one coordinate per DMPlex point” to “one coordinate per interpolation node/slot”, preventing coordinate-key collisions and enabling correct de-smearing with an order-appropriate line mass. It addresses Issue #459 by ensuring all trace nodes are reported and by making mass="auto" choose a consistent line-mass solve for degree ≥ 3 traces.

Changes:

• Reworked boundary node gathering to emit (point, slot, coord) per interpolation node, including true edge-interior node coordinates built via coordinate-field interpolation.
• Updated de-smearing to (a) structurally determine trace order from the section, (b) assemble a general consistent 1D line mass for higher degrees, and (c) hard-refuse coordinate-key collapse.
• Added serial + parallel regression tests covering degree sweep, P3 interior-node placement, collapse-guard behavior, and collective participation on ranks without boundary ownership.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 1 comment.

File	Description
src/underworld3/utilities/boundary_flux.py	Core fix: per-node identity/coordinates, collapse guard, and consistent line-mass support for degree ≥ 3 traces.
src/underworld3/cython/petsc_generic_snes_solvers.pyx	Updated boundary_flux docstring to reflect new mass="auto" behavior for degree ≥ 3 2D traces.
tests/test_1019_boundary_flux.py	Added regression tests for degree sweep correctness, P3 interior-node placement, and collapse guard.
tests/parallel/test_1065_boundary_flux_parallel.py	Added partition-independence test for degree-3 trace flux recovery and collective coordinate build.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.