3D: shape relaxation can break the custom-P FMG build (zero-column transfer) — silent GAMG fallback

[lmoresi]

Found while re-running the adaptivity demos with the new mesh.relax() / adapt(..., relax=...).

In 3D, relaxing an adapted mesh can make the mesh-owned custom-P geometric-MG hierarchy fail to build:

custom_mg: mesh-owned FMG build failed (transfer 6->7 has 26 zero columns
(coarse DOFs with no fine image) — BC-per-level reduction failed; coarse
operator would be singular.); using the solver default preconditioner.

The solver then silently falls back to GAMG — correct answers, but much slower, and nothing tells the user their geometric hierarchy is gone.

Measured

Poisson MG gate, 3D box, refinement=2:

config	cells	pc	its
adapt only	24594	mg	2
adapt + relax (at end)	24594	mg	2
adapt + relax "per-generation"	24160	gamg	23

and in the larger combine_mmpde_nvb_3d.py demo (max_levels=2, ~36k cells), relax-at-end also failed — transfer 6->7 has 26 zero columns, pc=gamg, 26 its. So both placements are affected; whether it trips appears to depend on mesh size / depth rather than on the placement.

2D is unaffected in everything tested: pc=mg throughout, Stokes velocity block and Poisson, 2 its.

Likely mechanism

The custom-P prolongation is built by locating fine DOFs inside coarse cells (barycentric). Relaxation moves the fine level coordinates; some coarse DOF then has no fine DOF mapping to it, giving a zero column and a singular coarse operator. 3D is presumably more exposed because tetrahedral point location is less forgiving than triangular.

Why it matters

adapt(..., relax=True) is documented as safe for MG on the strength of 2D measurements. That claim is wrong in 3D and the docs are being corrected alongside this issue.

Suggested directions

1. The fallback should be loud — an adapted+relaxed mesh silently losing its geometric hierarchy is the kind of thing that shows up as "why is my 3D solve slow" months later.
2. The transfer builder could tolerate zero columns by falling back to nearest-neighbour / RBF for the orphaned coarse DOFs rather than rejecting the whole level (an RBF builder already exists).
3. Or relaxation could be constrained so no fine DOF leaves its coarse parent cell — a step cap tied to the coarse cell size.

Reproducer: ~/+Simulations/nvb_3d_adapt_evaluation/combine_mmpde_nvb_3d.py (case C), and the smaller isolation above.

[lmoresi]

Root-caused and fixed in feature/adaptivity-round3 (commit 5db0320).

Root cause — it is not really about relaxation

The custom-P prolongation is built geometrically, not topologically:

tri  = Delaunay(coarse_coords)        # re-triangulate the coarse DOF CLOUD
simp = tri.find_simplex(fine_coords)  # locate each fine DOF in one simplex

barycentric_prolongation discards the coarse mesh cells and never consults the parent/child refinement relation. That is deliberate — it is precisely what lets custom-P accept non-nested coarse/fine pairs — but it means a preserved hierarchy topology buys the transfer nothing, and the builder has only LOCAL support: a coarse DOF is reached only if some fine DOF lands in a Delaunay simplex touching it.

Move the fine coordinates and a coarse DOF can lose every fine image -> zero column in P -> PᵀAP has a zero row and column -> singular -> guard fires. Relaxation is simply the operation that moves them; anything else that perturbs node positions would do the same.

3D is the exposed case because a Delaunay triangulation of a graded 3D point cloud is sliver-ridden, each simplex has only four vertices (sparse per-DOF coverage), and proportionally more DOFs sit near the convex hull (where find_simplex returns -1 and the nearest-DOF fallback concentrates coverage further).

Fix

rbf_prolongation has GLOBAL support — every coarse DOF is reached through the RBF solve — so it cannot produce a zero column. The FMG build now retries with rbf before abandoning geometric MG, and warns that it did rather than degrading in silence.

case	before	after
3D adapt + relax="per-generation" (24k cells)	gamg, 23 its	pc=mg, 2 its
3D combine demo, transfer 6->7, 26 zero cols (~36k cells)	gamg, 26 its	pc=mg, 2 its

This is strictly better for any cause of zero columns, not just relaxation: the old path discarded the entire hierarchy when a cheaper degradation was available.

45 tests pass across test_1014/1015/1017/1020 (MG), test_0752 (relax) and test_0836 (graded adapt).

Still open

• Suggestion 1 from the original report — making the final fallback loud — is only partly addressed. The retry warns, and the terminal GAMG fallback still warns, but a user not reading warnings still gets a silent performance cliff. Worth considering whether losing geometric MG should be an error under an opt-in strict mode.
• No regression test yet pins the retry itself (the 45 above cover the surrounding behaviour). A test that constructs a zero-column transfer deliberately and asserts pc == "mg" would lock it in.

[lmoresi]

Superseded my own first fix — the RBF retry was a silent scalability cliff, which I should have checked before shipping it.

The RBF prolongation is fully dense

	barycentric nnz/row	rbf nnz/row
nc=200, nf=800	2.9–3.4	200.00
nc=800, nf=3200	3.0–3.8	800.00

nnz/row == n_coarse exactly — every fine DOF couples to every coarse DOF, so the Galerkin coarse operator PᵀAP is dense as well. It passed the tests only because those meshes are small. At production size that is ruinous, and silently so — the same class of failure as the GAMG fallback it was introduced to cure.

The narrower fix

The defect is not "wrong builder", it is a handful of orphaned columns. Point location has local support, so a coarse DOF is reached only if some fine DOF lands in a simplex touching it; move the fine coordinates and a coarse DOF can lose every fine image. Only those orphans need repair.

Each orphan now gets its nearest fine DOF as a pure injection (weight 1) — the same fallback the builder already applies to fine points outside the coarse hull. Partition of unity holds (row still sums to 1), sparsity holds (~d+1 nnz/row), the column is no longer empty, and distinct rows are used so two orphans cannot claim the same fine DOF.

Result on the case that motivated this issue:

	before	after
3D adapt + relax="per-generation"	gamg, 23 its	pc=mg, 2 its

and the RBF retry no longer fires at all on real adapt children.

What remains

• The retry is kept as a genuine last resort. A deliberately pathological synthetic pair (fine cloud covering 60% of the coarse hull) still leaves orphans the repair cannot place; that does not arise for a real coarse/fine mesh pair, where both cover the same domain.
• Docs now state plainly that the RBF path is a correctness rescue, not a scalable transfer.
• The principled answer for nested hierarchies remains Use the nested hierarchy for MG transfers on adapt children (topological prolongation) #425: use the exact parent/child prolongation and never do point location at all. That transfer is 1–1.4 nnz/row and structurally cannot produce a zero column.

52 tests pass across the four MG suites, relax, nested prolongation and graded adapt.