CBF boundary_flux: _assemble_volume_reaction is wrong on inhomogeneous Dirichlet boundaries (2D and 3D)

[lmoresi]

Found by the adversarial review of PR #404 (2026-07-24) — pre-existing, dimension-independent, NOT introduced by that PR.

_assemble_volume_reaction (the residual-based reaction used by boundary_flux) is exact on g=0 Dirichlet boundaries but produces garbage when the boundary value is non-zero:

• P1 tet box, unit conduction: Top (g=0) boundary reaction sums to −1.000000 (exact); Bottom (g=1) sums to −4.41 (P1) / −11.59 (P2) against a true flux of ±1.
• 2D box, Bottom g=1: sumR = −27.27, recovered 'flux' ≈ −29.5 for a true value of 1.
• Shell case: interior residual is O(1) even though the KSP converged to rnorm 1e-13 — the Dirichlet boundary VALUES are evidently not inserted into the residual evaluation state before DMPlexSNESComputeResidualFEM.

Existing regressions (test_1019 etc.) pass only because they measure homogeneous boundaries. Every published CBF heat-flux / Nusselt recipe that reads the flux on a non-zero-temperature boundary is affected.

Likely fix: scatter the constrained boundary values into the local state vector (boundary-FEM-correct scatter, the same class as the SNES-callback scatter fix in #? — see test_1016's driven-boundary machinery) before computing the residual-based reaction.

Underworld development team with AI support from Claude Code

[lmoresi]

Scope note: which reaction paths are affected (checked at code level, 2026-07-24)

NOT affected — the rotated-constraint path (rotated_bc.boundary_normal_traction, and therefore dynamic_topography and the FreeSurface manager in PR #375): its reaction is built from the ASSEMBLED rotated operator — A·u−b (linear) or the stashed final computeFunction residual (nonlinear). Rotated constraints are imposed by matrix modification in the rotated frame, so the constrained DOFs remain in the global system and U carries their true values — there is no essential-constraint stripping anywhere in that path. This stays true when the rotated datum becomes non-zero (#403's driven-normal extension): the datum lives in U and b, and A·u−b remains the exact reaction. (The #404 review incidentally validated this path's np4 behaviour to 1e-8.)

Affected — the residual path (_assemble_volume_reaction → boundary_flux CBF): any reading on or NEAR a classical non-zero Dirichlet boundary. The corruption is local: the stripped values sit at the g≠0 wall's nodes, so the residual is wrong only in cells touching that wall. Consequences for mixed models (e.g. free-surface boxes with plate-driven side walls, heated-base convection):

• flux/traction read ON the driven wall: garbage (the reported case);
• flux read on an adjacent boundary is fine EXCEPT at nodes whose support touches the driven wall — i.e. the corner columns where a free surface meets a kinematically driven side wall are polluted;
• natural/Nitsche boundaries are not constrained, so nothing is stripped at their own nodes — same corner rule applies if a g≠0 Dirichlet wall is adjacent.

Regression plan should therefore include, beyond the g=0/g=1 two-wall box: a corner case (driven side wall + free top, assert the top flux including its end nodes) — the end-node assertion is the one the wall-only tests miss.

Underworld development team with AI support from Claude Code

[lmoresi]

Fixed by PR #412, merged to development on 2026-07-24 after a public adversarial review round. (Manual close — PRs target development.)

Underworld development team with AI support from Claude Code