Avoid replicated 3D boundary-flux recovery

[gthyagi]

Summary

Fix the MPI memory scaling of 3D boundary-flux recovery without changing its reaction, trace-mass, gauge-removal, or local-return mathematics.

Problem

The existing 3D path uses MPI.allgather for nodal reactions and P1/P2 trace elements. Every rank therefore retains the complete global boundary mesh, assembles the same sparse surface-mass matrix, and runs the same consistent-mass solve. Memory grows with both global boundary size and rank count.

This became the limiting post-processing stage for high-resolution spherical rotated-free-slip models: the Stokes solve completed, but consistent P2 normal-traction recovery exhausted aggregate job memory.

Change

• gather reactions, trace elements, and each rank requested coordinate keys to rank zero;
• assemble and solve the global 3D boundary-mass system once;
• scatter only the recovered values requested by each rank, preserving local xs order;
• preallocate sparse COO arrays instead of appending Python scalar entries;
• broadcast root-side errors before scatter so all MPI ranks fail collectively;
• retain the existing 2D path unchanged.

The behavior of mass={"auto","lumped","consistent","p1"}, partial-reaction sum/overwrite semantics, mean removal, and P2 midpoint interpolation is unchanged.

Validation

• tests/test_1019_boundary_flux.py -k 3d: 3 passed.
• Four-rank tests/parallel/test_1065_boundary_flux_parallel.py -k 3d: 2 passed on every rank.
• Previous production validation of this exact commit completed the Zhong cellsize=1/64, 192-rank post-processing path in 7:29 with 301.77 GB peak memory; the replicated path had exhausted 576 GB.

Scope

This PR contains only the general boundary_flux gather/solve/scatter correction and its rotated_bc documentation update. Zhong geoid mathematics and benchmark-facing APIs remain isolated in #591.