Skip to content

fix(mmpde): per-cell collapse floor — mover was a silent no-op on graded meshes (#419) - #420

Merged
lmoresi merged 1 commit into
developmentfrom
bugfix/mmpde-graded-area-floor
Jul 27, 2026
Merged

fix(mmpde): per-cell collapse floor — mover was a silent no-op on graded meshes (#419)#420
lmoresi merged 1 commit into
developmentfrom
bugfix/mmpde-graded-area-floor

Conversation

@lmoresi

@lmoresi lmoresi commented Jul 25, 2026

Copy link
Copy Markdown
Member

Closes #419.

The bug

The MMPDE line-search accept test used a single absolute collapse floor derived from the median cell volume:

a0_own_med = _global_mean(float(np.median(np.abs(signed_vol(coords, cells_own)))))
a_min_floor = float(area_floor_frac) * a0_own_med
...
if np.all(np.isfinite(trial)) and _min_area(trial) > a_min_floor:

1. It disabled the mover on any graded mesh. The floor is 1% of the median cell volume. Anything out of mesh.adapt has finest cells orders of magnitude below the median, so they start already under the floor. Every trial step was rejected, the line search backtracked to scale=0, and the mover returned having moved nothing — silently:

mmpde outer 1/8: I=1.991126e+00 dI=+0.00e+00 scale=0.000 max|dx|=0.00e+00

2. It was partition-dependent. A _global_mean of per-rank medians is not the global median, so which cells fell foul of the floor depended on how cells were spread across ranks. The same mesh moved at np=1 and stood completely still at np=2.

This affected redistribute_nodes / node_redistribution generally, not just new code.

The fix

Floor each cell against its own starting volume: no cell may shrink below area_floor_frac of what it started as. Scale-free, grading-free, partition-independent, and still a strict no-fold certificate — a positive ratio certifies positive volume. area_floor_frac keeps its name and default; only its reference changes.

Tests

New tests/test_0751_mmpde_graded_area_floor.py (level_1, tier_a):

  • test_mover_moves_on_a_graded_mesh — asserts the mesh is genuinely graded first, so the test cannot pass without exercising the floor, then asserts the mover actually moves.
  • test_collapse_guard_still_rejects_folding — the guard must still do its job: no cell folds, and none collapses below the documented fraction of its own start.

54 existing tests pass across the mover, follow_metric, graded adapt and mesh smoothing.

Reviewer notes

  • The min(initial=np.inf) in _min_area_frac is deliberate: a rank owning no cells must contribute the identity for the global MIN rather than raising on an empty reduction.
  • Worth checking: area_floor_frac=0.01 now means something stricter in the fine region than it used to (1% of a small cell, not 1% of the median), so a previously-accepted step that shrank a fine cell hard will now be backtracked. That is the intended reading, but it is a behaviour change for well-graded meshes that were already moving.
  • Not addressed here: the mover still reports no warning when the line search backtracks to scale=0 for any other reason. A silent no-op is hard to diagnose — probably worth a follow-up warning.

Underworld development team with AI support from Claude Code

…ded meshes (#419)

The MMPDE line-search accept test used a single ABSOLUTE collapse floor,
`area_floor_frac x median cell volume`, with the median taken as a
`_global_mean` of per-rank medians. Two defects followed.

Any graded mesh disabled the mover outright. Anything out of `mesh.adapt`
has finest cells orders of magnitude below the median, so those cells
start already under the floor; every trial step failed the test, the line
search backtracked to `scale=0`, and the mover returned having moved
nothing — with no warning:

    mmpde outer 1/8: I=1.991126e+00 dI=+0.00e+00 scale=0.000 max|dx|=0.00e+00

And the behaviour depended on the partition, since a mean of per-rank
medians is not the global median: the same mesh moved at np=1 and stood
completely still at np=2.

Floor each cell against its OWN starting volume instead — no cell may
shrink below `area_floor_frac` of what it started as. Scale-free,
grading-free, partition-independent, and still a strict no-fold
certificate (a positive ratio certifies positive volume). The kwarg
keeps its name and default; only its reference changes.

This affected `redistribute_nodes` / `node_redistribution` generally,
not just new code. Found while building shape relaxation on the mover.

Underworld development team with AI support from Claude Code
Copilot AI review requested due to automatic review settings July 25, 2026 10:27

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

@lmoresi

lmoresi commented Jul 25, 2026

Copy link
Copy Markdown
Member Author

Adversarial review (author, pre-review)

Trying to break my own fix. Four things a reviewer should push on.

1. area_floor_frac silently changes meaning, and 0.01 may now be too strict. Previously it was 1% of the median cell volume — for a fine cell in a graded mesh that was a very loose bound (the cell could shrink to nothing and still clear it). Now it is 1% of that cell's own volume. Any workflow that legitimately compresses a cell by more than 100x within a single mover call will now backtrack where it previously did not. Free-surface collapse and strong re-grading are the plausible victims. I have not tested those paths. If someone is relying on large single-call compression, this is a regression for them and the default may need to be looser (1e-3?).

2. It caps how much redistribution one call can do. A cell cannot shrink below 1% of its entry volume in this call, no matter what the metric asks. Repeated calls compose fine, but a single redistribute_nodes onto a very aggressive metric is now bounded in a way it was not. Arguably correct (that is what a collapse guard is for) but it is a new constraint, not a restoration of intended behaviour.

3. Degenerate input cells. np.maximum(a0_own, 1e-300) means a cell that arrives with ~zero volume gets an essentially infinite ratio, so the floor never binds for it. Folding is still caught (the ratio goes negative), but a sliver that arrives degenerate is effectively exempt from the collapse half of the guard. That is probably the right call — you cannot floor against zero — but it should be a conscious decision, not an accident of the clamp.

4. The test does not actually reproduce the original parallel failure. The sharpest form of this bug was partition-dependent: np=1 moved, np=2 did not. My test is serial and only asserts "a graded mesh moves". It would catch a reintroduced median-based floor, but it would not catch a re-introduction of a partition-dependent reduction. A proper regression would compare max|dx| at np=1 vs np=2. I did verify that by hand (both now move, energies agree to ~3 digits) but did not automate it, and the repo has form for rank-local diagnostics slipping back in.

Also not fixed here: the mover still emits no warning when the line search backtracks to scale=0. That silence is what made this bug expensive to find — it presented as "adaptation does nothing" with no diagnostic. Worth a follow-up.

What I am confident about: the no-fold guarantee is preserved (a positive volume ratio certifies positive volume), the reduction is partition-independent by construction, and the empty-rank path is handled (min(initial=np.inf)).

lmoresi added a commit that referenced this pull request Jul 27, 2026
…llback loud

Two review follow-ups.

STALENESS / MERGE HAZARD. This branch was cherry-picked from a branch that
already carried the per-cell collapse floor, so it duplicates PR #420. The
two copies were functionally identical but differed in em-dashes and a
missing issue reference, which would have produced a pointless
comment-only conflict when #420 lands. The block is now byte-identical to
#420's, so the rebase after #420 merges is a clean no-op.

Merge order is therefore: #420 first, then rebase this branch — the
duplicate hunks drop out silently rather than conflicting.

DENSE FALLBACK NOW SAYS SO. The RBF retry rescues correctness when the
barycentric orphan repair cannot place every coarse DOF, but it returns a
FULLY DENSE transfer (nnz/row == n_coarse), which makes the Galerkin
coarse operators dense. On the small meshes where it has fired that is
harmless; at production size it is a silent performance cliff of exactly
the kind the GAMG fallback was. The warning now says this outright and
tells the reader to fix the cause rather than accept the symptom.

Underworld development team with AI support from Claude Code
lmoresi added a commit that referenced this pull request Jul 27, 2026
…llback loud

Two review follow-ups.

STALENESS / MERGE HAZARD. This branch was cherry-picked from a branch that
already carried the per-cell collapse floor, so it duplicates PR #420. The
two copies were functionally identical but differed in em-dashes and a
missing issue reference, which would have produced a pointless
comment-only conflict when #420 lands. The block is now byte-identical to
#420's, so the rebase after #420 merges is a clean no-op.

Merge order is therefore: #420 first, then rebase this branch — the
duplicate hunks drop out silently rather than conflicting.

DENSE FALLBACK NOW SAYS SO. The RBF retry rescues correctness when the
barycentric orphan repair cannot place every coarse DOF, but it returns a
FULLY DENSE transfer (nnz/row == n_coarse), which makes the Galerkin
coarse operators dense. On the small meshes where it has fired that is
harmless; at production size it is a silent performance cliff of exactly
the kind the GAMG fallback was. The warning now says this outright and
tells the reader to fix the cause rather than accept the symptom.

Underworld development team with AI support from Claude Code
@lmoresi
lmoresi merged commit 2d5e5d3 into development Jul 27, 2026
2 checks passed
@lmoresi
lmoresi deleted the bugfix/mmpde-graded-area-floor branch July 27, 2026 04:31
@lmoresi
lmoresi restored the bugfix/mmpde-graded-area-floor branch July 27, 2026 04:31
@lmoresi
lmoresi deleted the bugfix/mmpde-graded-area-floor branch July 27, 2026 04:32
lmoresi added a commit that referenced this pull request Jul 27, 2026
…ldren (#433)

* fix(mmpde): per-cell collapse floor — mover was a silent no-op on graded meshes (#419)

The MMPDE line-search accept test used a single ABSOLUTE collapse floor,
`area_floor_frac x median cell volume`, with the median taken as a
`_global_mean` of per-rank medians. Two defects followed.

Any graded mesh disabled the mover outright. Anything out of `mesh.adapt`
has finest cells orders of magnitude below the median, so those cells
start already under the floor; every trial step failed the test, the line
search backtracked to `scale=0`, and the mover returned having moved
nothing — with no warning:

    mmpde outer 1/8: I=1.991126e+00 dI=+0.00e+00 scale=0.000 max|dx|=0.00e+00

And the behaviour depended on the partition, since a mean of per-rank
medians is not the global median: the same mesh moved at np=1 and stood
completely still at np=2.

Floor each cell against its OWN starting volume instead — no cell may
shrink below `area_floor_frac` of what it started as. Scale-free,
grading-free, partition-independent, and still a strict no-fold
certificate (a positive ratio certifies positive volume). The kwarg
keeps its name and default; only its reference changes.

This affected `redistribute_nodes` / `node_redistribution` generally,
not just new code. Found while building shape relaxation on the mover.

Underworld development team with AI support from Claude Code

* feat(meshing): mesh.relax() — shape relaxation in the ideal reference frame

Refinement chooses node positions from combinatorics — which edge the
tagging rule nominated, or the cell centroid — never from geometry. A
refined mesh therefore carries needles and slivers that reflect the base
mesh's arbitrary choices rather than the problem.

The MMPDE mover could not fix this, and the reason was its REFERENCE
FRAME, not its metric: the reference is the mesh at first call, so
"distortion" meant "deviation from the mesh I was handed" and a distorted
mesh was its own optimum. `redistribute_nodes(1)` moved exactly nothing.

Add `reference=` to the mover:

  "mesh"          (default, unchanged) the redistribution frame
  "ideal"         one regular simplex, scaled per cell to that cell's own
                  volume => r = 1 at entry, so size is held and only shape
                  moves. Pure shape repair.
  "ideal-metric"  uniform reference volume; the metric sets size. Needed
                  when the sizes we were handed are themselves the problem.

Deviation is then measured from equilateral-in-the-metric, so distortion
the mesh arrived with is no longer self-justifying. Same functional, same
non-folding guarantee, same parallel assembly.

Exposed as `mesh.relax()` (no metric => "ideal"; with one =>
"ideal-metric"), and as `adapt(metric, max_levels=N, relax=True)`, which
relaxes each generation INSIDE the loop — before the next pass marks from
it, which is where the gain comes from.

Measured, 2D, gmsh base + NVB adapt, fault-localised P1 interpolation
error (Stokes velocity block on custom-P FMG in every case):

  adapt only               832 cells  2.746e-2        4 vel-KSP its
  + relax at end           832 cells  2.559e-2 (-7%)  3 vel-KSP its
  + relax per generation   856 cells  2.384e-2 (-13%) 4 vel-KSP its

Relax-per-generation wins on error and on error-per-cell for 3% more
cells. MG is unaffected: full hierarchy, no GAMG fallback, solutions
agree to 4 s.f. — the operators need only the topological hierarchy, and
the coarse levels keep their own coordinates.

3D IS A WASH (+0.5%, unchanged at 120 iterations) and is documented as
such; this is a 2D win only.

Underworld development team with AI support from Claude Code

* docs(relax): 3D improves QUALITY not error — correct the earlier 'wash' claim

Measured on an adapted 3D mesh: cells with q<0.1 halve (3.6% -> 1.8%),
median quality 0.32 -> 0.39, p99 dihedral 153 -> 146 degrees. The
earlier note said 3D was 'a wash' on the strength of an unchanged
interpolation error, conflating accuracy with mesh quality. The error
really is unchanged, and that is consistent — relaxation holds the size
distribution, and the 2D error gain came from cells aligning onto the
feature, which an isotropic 3D metric gives no reason to do.

Also refutes the base-provenance hypothesis: relaxation repairs the
PETSc uniform 1:8 split's damage well (q p10 0.449 -> 0.526 on a
uniformly split base), so the fixed reference-octahedron diagonal is
not what was limiting it.

Underworld development team with AI support from Claude Code

* feat(adapt): relax= takes a placement — True (at end, default) or per-generation

Both placements beat no relaxation, and neither dominates the other: at-end
keeps the cell count identical to the unrelaxed mesh and cuts the far-field
closure halo most (42.8 -> 27.7 excess cells), per-generation gives the
cleanest element shapes (99th pct max angle 112 -> 96 deg) and the lowest
interpolation error, at ~3% more cells. Rather than manufacture a ranking
the evidence does not support, both are offered and documented side by
side, with the cheaper, less invasive one as the default.

  adapt(metric, max_levels=N, relax=True)              # relax once, at end
  adapt(metric, max_levels=N, relax="per-generation")  # relax in the loop

Also fixes a latent bug from when `relax` was a boolean: both in-loop call
sites passed `None if relax is True else relax` as the mover's METRIC, so
`relax="per-generation"` fed the mode string to sympify. `relax` is a mode,
never a metric; the metric is now resolved once (a callable metric cannot
go to the mover, so those relax in the pure shape frame at fixed size).

Underworld development team with AI support from Claude Code

* test+docs(relax): cover both relax hooks, document the frame trade-off (#422)

Tests (tests/test_0752, level_1 tier_a, 7 cases): topology and no-fold
invariants, grading preservation, the max-angle-tail improvement, both
`adapt(relax=...)` placements, and the argument guards. Includes a 3D case
because 3D goes through the CELL-LIST refinement engine — the native uwnvb
transform is 2D-only — so that second relax hook was otherwise untested.

Writing the tests turned up a real distinction I had been eliding: relax()
and relax(metric) are different operations. The metric-free call is a
shape guarantee; passing a metric switches to the ideal-METRIC frame,
which re-grades as well as reshapes and will trade element shape away to
chase the size field. Measured on a 4-level graded box, 99th-percentile
max angle:

    no metric      117.9 -> 113.8   (improves)
    with metric    117.9 -> 127.4   (WORSE)
    per-generation 117.9 -> 101.2   (best)

The first version of the test asserted the improvement while passing a
metric, and failed. Documented in the docstring and the new subsystem
page rather than papered over.

New docs/developer/subsystems/mesh-shape-relaxation.md covers the
reference-frame explanation, both placements with the measured table, the
MG result, and an explicit "what relaxation cannot do" section (does not
rescue centroid refinement, cannot flatten discrete bisection size steps,
does not concentrate cells). 3D is documented as a quality win only, with
the anisotropic-metric question tracked in #422.

Underworld development team with AI support from Claude Code

* docs(relax): correct the "MG is unaffected" claim — 3D can fail the custom-P build (#424)

Re-running the adaptivity demos with relaxation turned on exposed a
failure I had ruled out on 2D-only evidence. In 3D, relaxing an adapted
mesh can make the mesh-owned custom-P FMG hierarchy fail to build:

    custom_mg: mesh-owned FMG build failed (transfer 6->7 has 26 zero
    columns (coarse DOFs with no fine image) ...); using the solver's
    default preconditioner.

and the solver then SILENTLY falls back to GAMG. Measured on a 3D
Poisson gate: adapt-only and adapt+relax-at-end gave pc=mg at 2 its,
relax="per-generation" gave pc=gamg at 23; in a larger demo relax-at-end
failed too, so both placements are affected and it appears to depend on
mesh size and depth. 2D is unaffected throughout (Stokes and Poisson).

The prolongation locates fine DOFs inside coarse cells; relaxation moves
the fine coordinates and some coarse DOF ends up with no fine image ->
zero column -> singular coarse operator.

My earlier claim came from measuring only 2D and generalising. The docs
and the Mesh.relax docstring now state the 3D caveat and point at #424.

Underworld development team with AI support from Claude Code

* fix(custom-mg): retry with the global-support RBF builder before dropping to GAMG (#424)

The custom-P prolongation is built GEOMETRICALLY, not from the refinement
relation: `barycentric_prolongation` re-triangulates the coarse DOF point
cloud (`Delaunay(coarse_coords)`) and locates each fine DOF in one simplex
of it. That is deliberate — it is what lets custom-P accept non-nested
coarse/fine pairs — but it means a preserved 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 simplex touching it.

Move the fine coordinates (which is exactly what mesh.relax does) and a
coarse DOF can lose every fine image. Its column in P goes to zero, the
Galerkin coarse operator PᵀAP acquires a zero row and column, and it is
singular. The build then abandoned geometric MG entirely and fell back to
the solver's default preconditioner — silently, so it presented only as a
slow solve. 3D is the exposed case: Delaunay of a graded 3D cloud is
sliver-ridden, each simplex has four vertices, more DOFs sit near the hull.

`rbf_prolongation` has GLOBAL support — every coarse DOF is reached through
the RBF solve — so it cannot produce a zero column. Retry with it before
giving up, and warn that we did rather than degrade in silence.

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

This is strictly better for ANY cause of zero columns, not just relaxation:
dropping to GAMG discarded the whole hierarchy when a cheaper degradation
was available. Docs and the Mesh.relax docstring now carry the mechanism,
replacing my earlier "MG is unaffected" claim (which generalised from
2D-only measurements).

Underworld development team with AI support from Claude Code

* docs(design): two MG transfer paths — use the nested hierarchy we already pay for

The adapt hierarchy is maintained exactly (that is what NVB conforming
closure and its halo cost buy) and the MG transfer then discards it,
re-deriving an approximation by Delaunay point location on the coarse DOF
cloud. That is the direct cause of #424: the barycentric builder has only
LOCAL support, so moving the fine coordinates can leave a coarse DOF with
no fine image and a singular coarse operator.

The nested transfer is trivial and cannot go singular: every fine vertex
is either an inherited coarse vertex (weight 1) or the midpoint of a
coarse edge (weights 1/2, 1/2). Structurally full rank, immune to node
motion, 1-2 nonzeros per row. Both engines already record it (edge2mid in
nvb.py; the native transform's split-edge bookkeeping), and our own
petsc_dm_set_regular_refinement wrapper — the flag for PETSc's exact
nested path — is never called.

Keeps the geometric builder for the genuinely non-nested pairs it was
asked for (moved base vs child, unrelated meshes, the planned
seed-point variant).

Records the ordering trap: the engine-id -> DM-point bridge is coordinate
matching, which is exact only before snapping/relaxation moves midpoints
off their parents. The native path has a clean window; the cell-list path
does not, so to_dm should return its id map instead.

Flags the accuracy-vs-robustness trade as something to MEASURE on curved
cases rather than assume.

Underworld development team with AI support from Claude Code

* feat(adapt): record the EXACT nested MG prolongation per refinement pass (#425)

Step 1 of using the hierarchy we already pay to maintain. adapt() now
records, for every refinement pass, the true P1 embedding of the coarse
space in the fine one — built from the parent/child relation rather than
re-derived by Delaunay point location on the coarse DOF cloud.

Every fine vertex is either an inherited coarse vertex (weight 1 on itself)
or the midpoint of a coarse edge (1/2, 1/2). Measured on the recorded
transfers, 2D and 3D:

  partition of unity   rowsum in [1.000000000000, 1.000000000000]
  zero columns         0 in every pass  <- makes #424 impossible here
  linear reproduction  max|P x_c - x_f| = 0 (worst 1.1e-16)
  sparsity             1.0-1.4 nnz/row vs dim+1 for point location

Two things this had to get right:

A refine PASS is not one bisection LEVEL. The conforming closure cascades,
so a vertex can be the midpoint of an already-split half-edge and sit at the
QUARTER point of the original coarse edge — it matches neither a coarse
vertex nor a coarse edge midpoint (observed: 3.3e-2 off, a quarter edge).
Resolution is therefore iterative: seed from coarse vertices and coarse edge
midpoints, then repeatedly resolve any vertex that is the midpoint of two
already-resolved fine neighbours, composing weights. Exact at any depth.

The capture must happen in the pristine window. The engine-id/DM-point
bridge is coordinate matching, exact only before snapping or relaxation
moves midpoints off their parents. On the native path that window is
between _nvbx.refine() and snap_level_boundaries(); the capture goes there.
The cell-list engine records its map in to_dm() instead, where the two
numberings meet.

Also corrects a claim I committed earlier: 3D does NOT go through the
cell-list engine. Both 2D and 3D prefer the native uwnvb transform when it
is built (cell-list is the np=1 fallback), so the 3D relax test does not
cover the cell-list hook as its docstring claimed. That hook is still
uncovered; the docstring now says so.

Not yet wired into custom_mg — consuming these as a "topological" builder
is the next step, and the accuracy-vs-robustness question on curved
boundaries still needs measuring before it becomes the default.

Underworld development team with AI support from Claude Code

* fix(custom-mg): repair orphaned coarse DOFs in place, keeping the transfer sparse (#424)

The RBF retry added earlier rescues correctness but is not a scalable fix,
and I shipped it without checking its density. Measured:

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

nnz/row == n_coarse exactly: the RBF prolongation is FULLY DENSE, so the
Galerkin coarse operator PᵀAP is dense too. It passed the tests only
because those meshes are small; at production size it would be ruinous.
That is a silent performance cliff of the same kind as the GAMG fallback
it was meant to cure.

The actual defect is narrower than "wrong builder". 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 fixing, not the whole transfer.

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

  3D adapt + relax="per-generation":  gamg 23 its -> pc=mg 2 its, and the
  RBF retry no longer fires at all.

The retry stays as a last resort for pathological pairs (a synthetic cloud
covering 60% of the coarse hull still leaves orphans), with the docs now
stating plainly that it is a rescue and not a scalable transfer.

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

Underworld development team with AI support from Claude Code

* feat(custom-mg): use the EXACT nested transfer for degree-1 fields (#425)

Closes the gap where adapt() maintained an exact hierarchy and the MG
transfer then re-derived an approximation to it by Delaunay point location.
For the adapt generations the recorded parent/child relation is now
installed directly; the uniform coarse tail beneath them is not a bisection
hierarchy, so it keeps using the geometric builder.

  2D adapt          nested 2, builder 2   pc=mg its=3  err 1.6e-13
  3D adapt          nested 3, builder 2   pc=mg its=3  err 9.1e-13
  3D adapt + relax  nested 3, builder 2   pc=mg its=4  err 4.5e-12

On this path the zero-column failure (#424) cannot occur at all: every
coarse DOF carries weight 1 into its own inherited fine DOF, so the
transfer is structurally full rank. It is also sparser — 1-1.4 nnz/row
against d+1 for point location — and immune to node motion, so relaxation
and snapping cannot disturb it.

Verified before relying on it: `_get_coords_for_basis(1, True)` returns
coordinates in DM vertex-stratum order, identical arrays in 2D and 3D, base
and child. Had that not held, the recorded indices would have installed a
silently WRONG transfer rather than failing.

RESTRICTED TO degree == 1. The recorded relation is vertex-level; a P2 field
also has edge DOFs it says nothing about, so higher degree still goes
through the builder, and a test pins that it declines rather than
half-applies. Stage 2 (#425) is the any-degree extension: record the parent
CELL, then evaluate the coarse basis at the fine DOF reference coordinates —
exact for any degree and still orphan-free, since every coarse cell has
children. Locating cell CENTROIDS is robust where locating DOFs is not:
centroids are strictly interior, DOFs sit on cell boundaries and drift.

Shapes are checked rather than assumed at every level; a mismatch falls back
to the builder rather than installing a wrong transfer.

59 tests pass across the four MG suites, relax, nested prolongation, graded
adapt and variable transfer.

Underworld development team with AI support from Claude Code

* docs(custom-mg): correct the "FE-exact" claim on the barycentric builder (#425)

The docstring said each fine DOF is interpolated by "the coarse element that
contains it". That is only true when Delaunay of the coarse DOF cloud
coincides with the coarse mesh. Measured agreement:

  2D uniform base, P1     168 cells    168 Delaunay    100 %
  2D adapt child, P1      238          238            90.8 %
  3D uniform base, P1     800          812            58.8 %
  3D adapt child, P1     4685         5201            17.1 %
  3D adapt child, P2     4685        39597              --

It triangulates the coarse DOF POINT CLOUD from scratch and never consults
the mesh cells, so it is linear interpolation over an ad-hoc triangulation,
not the FE embedding. Linear reproduction is why multigrid converges well;
but the coarse P2 space is not represented exactly, and the 58.8 % on a
UNIFORM base reflects PETSc 1:8 tetrahedral splits not being Delaunay — the
same fixed-diagonal property as the 148 degree dihedral measured earlier.

Also records a predicted (not yet measured) risk: Delaunay simplices do not
respect mesh topology, so on a non-convex domain or one with an internal
interface a simplex can bridge the discontinuity and smear the coarse
correction across it.

Underworld development team with AI support from Claude Code

* refactor(adapt): remove scipy from the nested path — exact identification, no search

Maintainer ruling: the nested-mesh case must not use the general geometric
mapping, and must not pull in meshing/interpolation dependencies. "If you are
using scipy it is probably wrong."

I had introduced two violations myself:

  nvb.nested_prolongation_from_dms   scipy.spatial.cKDTree + a tolerance
  _adapt_nested._relax_generation    scipy.spatial.cKDTree nearest-neighbour

Both are now exact byte-equality lookups. This is not a weaker version of the
search — it is the correct operation. Every fine vertex of a bisection pass is
an inherited coarse vertex or the exact average of two already-known vertices,
computed by the same float arithmetic on both sides, so equality is exact.
Measured: 100% of 2D fine vertices match a coarse vertex or coarse edge
midpoint bit-for-bit; the 3D remainder (27%) are the closure-cascade quarter
points, which resolve exactly against already-resolved FINE neighbours by the
same test. No spatial index, no tolerance parameter, no nearest-neighbour
guess that could silently pick the wrong point.

The `tol=` argument is gone with it — it was a fragility, not a tunable.

Unchanged and re-verified after the refactor: 2D 2 passes, 3D 6 passes, all
captured, partition of unity exact, zero zero-columns, 1.0-1.4 nnz/row; MG
still installs the nested transfer (2D 2, 3D 3) with pc=mg and errors at
round-off. 46 tests pass.

Still using scipy, PRE-EXISTING and reported not silently fixed:
  nvb.py:596, :992                  cKDTree in to_dm (defensive; the code
                                    already documents the map as identity)
  discretisation_mesh.py:2499,:2824 cKDTree
  custom_mg.py                      Delaunay + cKDTree (the geometric path)
UW3 has its own uw.kdtree.KDTree, so scipy is not needed even where a real
spatial search is wanted.

Underworld development team with AI support from Claude Code

* refactor(nvb): drop scipy from to_dm — exact coordinate lookup, not a spatial query

`to_dm` used a scipy cKDTree nearest-neighbour search to map DM vertices back
to engine vertex ids. The DM is built by createFromCellList from those very
coordinates, so the two arrays hold bit-identical values: equality is the
correct test and a spatial query was standing in for an identity lookup.

Worse, it failed in the wrong direction. A nearest-neighbour search always
returns something, so if the assumption it was guarding ever broke it would
silently bind the WRONG vertex. The replacement raises instead, naming the
cause (coordinates modified between building the DM and reading it back).

nvb.py is now free of scipy entirely.

Also records the MG convergence check that preceded this cleanup. Poisson,
manufactured sin solution (the earlier linear-solution test was worthless —
the coarse space represents a linear field exactly, so every transfer
converged in 1-3 iterations regardless of quality):

  2D, 303 -> 7713 dofs      its 3, 3, 2, 2   nested == geometric
  3D, k=8pi, rtol=1e-12     its 5, 6, 5      geometric 4, 5, 4

2D is mesh-independent and identical. 3D costs the nested transfer ONE extra
iteration, consistently across three configurations — not noise. The likely
cause is stencil width: the nested transfer is near-injection (1.0-1.4
nnz/row) where barycentric spreads over ~4, and prolongation width affects
how well the coarse correction damps high-frequency error. Exactness is not
the same as effectiveness. Partly offset by the nested coarse operators being
~3x sparser, so each iteration is cheaper; net wall-clock not measured.

55 tests pass across all NVB engines (2D/3D, serial/parallel), nested
prolongation, relax and SBR adapt.

Underworld development team with AI support from Claude Code

* fix(adapt): restore round-2 composition — it was correct; the REFERENCE was wrong

I removed the cascade-vertex composition from nested_prolongation_from_dms on
the strength of a test showing 2-12 wrong DOFs per 3D pass with O(1) errors,
and wrote a confident commit message explaining why it was unsound. Both the
evidence and the reasoning were wrong.

The reference was uw.function.evaluate, which returns incorrect values at
points lying exactly on cell edges (#432) — precisely where every bisection
vertex sits. Checked against an INDEPENDENT reference instead (the P1 value
along the coarse edge each vertex lies on, computed directly):

  2D  3 passes   111/138/151 vertices   0 wrong   worst 2.3e-15
  3D  6 passes   334..2225 vertices     0 wrong   worst 2.9e-15

Every vertex, including every cascade generation. The transfer is the coarse
P1 embedding and always was.

For each disagreement the point lay on exactly ONE coarse edge, and the
edge-truth matched the transfer, not evaluate:

  dof  504   transfer +0.700344   evaluate +0.923709   edge-truth +0.700344
  dof 2177   transfer -2.122083   evaluate -6.128067   edge-truth -2.122083

Adds test_reproduces_an_arbitrary_coarse_field, which is the test that should
have existed from the start. Reproducing a globally LINEAR field is far too
weak — any local averaging passes it — so the original test could not
distinguish a correct embedding from a mis-attributed one in either
direction. A random coarse nodal field can, and the reference is computed
independently rather than via evaluate.

Also records the parent-CELL map where it is unambiguous (foundation for the
any-degree transfer, #425). It declines on 3 of 6 3D passes: value-correctness
of the vertex transfer does not imply unique cell attribution, so stage 2
still needs the true parent relation from the refinement transform.

Underworld development team with AI support from Claude Code

* chore: align the duplicated floor fix with #420 and make the dense fallback loud

Two review follow-ups.

STALENESS / MERGE HAZARD. This branch was cherry-picked from a branch that
already carried the per-cell collapse floor, so it duplicates PR #420. The
two copies were functionally identical but differed in em-dashes and a
missing issue reference, which would have produced a pointless
comment-only conflict when #420 lands. The block is now byte-identical to
#420's, so the rebase after #420 merges is a clean no-op.

Merge order is therefore: #420 first, then rebase this branch — the
duplicate hunks drop out silently rather than conflicting.

DENSE FALLBACK NOW SAYS SO. The RBF retry rescues correctness when the
barycentric orphan repair cannot place every coarse DOF, but it returns a
FULLY DENSE transfer (nnz/row == n_coarse), which makes the Galerkin
coarse operators dense. On the small meshes where it has fired that is
harmless; at production size it is a silent performance cliff of exactly
the kind the GAMG fallback was. The warning now says this outright and
tells the reader to fix the cause rather than accept the symptom.

Underworld development team with AI support from Claude Code
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants