Skip to content

Mesh shape relaxation, and the exact nested MG transfer for adapt children - #433

Merged
lmoresi merged 16 commits into
developmentfrom
feature/mesh-relax-and-nested-mg
Jul 27, 2026
Merged

Mesh shape relaxation, and the exact nested MG transfer for adapt children#433
lmoresi merged 16 commits into
developmentfrom
feature/mesh-relax-and-nested-mg

Conversation

@lmoresi

@lmoresi lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member

Follows the merged round-3 work (#409). Rebased onto current development — the old feature/adaptivity-round3 branch had gone stale behind #416/#404/#412.

Closes nothing outright; advances #424 and #425, and depends on #420 being merged for the mover to work at all on graded meshes.

What this adds

mesh.relax() — shape relaxation. Refinement picks node positions from combinatorics (which edge the tagging rule nominated), never from geometry, so a refined mesh carries needles that reflect the base mesh rather than the problem. The MMPDE mover could not fix this because its reference frame is the mesh it was handed, making a distorted mesh its own optimum — redistribute_nodes(1) moved exactly nothing. Adding reference="ideal" (a regular simplex, per-cell volume) measures deviation from equilateral-in-the-metric instead.

Exposed as mesh.relax() and adapt(metric, max_levels=N, relax=...) with two placements, neither dominant:

property unrelaxed relax at end relax per generation
P1 interpolation error 2.75e-2 2.56e-2 2.38e-2
99th pct max angle 112 deg 110 deg 96 deg
on-fault size spread 2.80 2.31 2.31
far-field closure halo 42.8 27.7 43.9

Default is at-end (cheaper, cannot change topology). Documented side by side rather than ranked.

Exact nested MG prolongation. adapt() now records the true coarse-to-fine embedding per refinement pass, built from the parent/child relation rather than re-derived by point location. Verified: partition of unity exact, zero zero-columns, exact reproduction of an arbitrary coarse field, 1.0-1.4 nnz/row. Installed for degree-1 fields; higher degree still uses the geometric builder (stage 2, #425).

Orphan repair in the geometric builder (#424). Point location has local support, so moving fine coordinates can leave a coarse DOF with no fine image -> zero column -> singular PtAP. Each orphan now takes its nearest fine DOF as a pure injection, preserving partition of unity and ~d+1 sparsity. 3D relax="per-generation": gamg 23 its -> pc=mg 2 its.

scipy removed from the nested path. nvb.py is now scipy-free; the coordinate lookups were exact-equality identifications wearing a nearest-neighbour search as a disguise, and a kNN search always returns something — it would have silently bound the wrong vertex if its assumption broke.

Testing

41 tests pass across relax, nested prolongation, graded adapt and the custom-MG suites. New: test_0752 (7 cases), test_0753 (11 cases).

Adversarial review to follow as a comment before this is marked ready.

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 27, 2026 01:47

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 27, 2026

Copy link
Copy Markdown
Member Author

Adversarial review (author, pre-review)

Trying to break my own work. Six things a reviewer should push on, roughly in order of how much they would cost if I am wrong.

1. The nested transfer is validated against an oracle I built myself. test_0753::test_reproduces_an_arbitrary_coarse_field computes the reference by locating the coarse edge each vertex lies on and taking the P1 value along it. That is independent of uw.function.evaluate on purpose — evaluate is wrong at exactly these points (#432) — but it shares my assumption that every bisection vertex lies on a coarse edge. The test asserts checked > 0.9 * n_vertices, and in 3D 1-2 vertices per pass are not on any coarse edge and go unchecked. If those are the ones that are wrong, this test cannot see it.

2. I got the correctness of this code wrong twice today, in both directions. I first claimed the transfer was the exact FE embedding on the strength of a linear-reproduction test that cannot distinguish a correct embedding from a mis-attributed one. I then removed a working component after a test using evaluate as the reference showed 2-12 wrong DOFs per 3D pass, and wrote a confident and entirely wrong explanation of why it was unsound. Both were caught only by building a third, independent oracle. Weight my claims here accordingly — the code is now measured, but my track record on reasoning about it unaided is poor.

3. The stage-1 restriction is narrower than it sounds. degree == 1 means P1 Poisson-like problems. Stokes never touches this path — velocity is P2, and instrumentation shows the nested transfer declined 3 times out of 3. Advection-diffusion typically wants P3. So the headline "exact nested transfer" currently does nothing for either production workhorse. A reviewer might reasonably ask whether it is worth merging before stage 2; my answer is that it is the reference implementation the geometric path is priced against, but that is a judgement call.

4. The nested transfer is measurably WORSE on iterations in 3D. Poisson P1, k=8pi, rtol=1e-12: nested 5,6,5 against geometric 4,5,4 — consistent across three configurations, not noise. Likely stencil width (nested is 1.0-1.4 nnz/row, near-injection; barycentric spreads over ~4). Its coarse operators are ~3x sparser so it may still win on wall-clock, but I have not measured wall-clock. If someone merges this expecting a speedup they may be disappointed.

5. The orphan repair is incomplete and I know it. On a deliberately pathological pair (fine cloud covering 60 % of the coarse hull) orphans remain and the RBF retry still fires — and that retry produces a fully dense transfer (nnz/row == n_coarse), which would be ruinous at production size. It does not arise for genuine coarse/fine mesh pairs, where both cover the same domain, but the failure mode is a silent performance cliff rather than an error.

6. Unexercised code paths. The cell-list refinement engine keeps its own edge2mid-based capture and its own relax hook; neither is covered, because both 2D and 3D prefer the native transform when it is built. I previously claimed the 3D test covered the cell-list hook — that was wrong and is corrected in the test docstring, but the coverage gap is real. Also _adapt_parent_cells is recorded but consumed by nothing, and declines on 3 of 6 3D passes.

Dependencies: this needs #420 (per-cell collapse floor) merged, or the mover silently does nothing on any graded mesh and relax() is a no-op.

What I am confident about: the recorded transfers satisfy partition of unity, have zero zero-columns, and reproduce an arbitrary coarse field to ~1e-15 across every pass in 2D and 3D; adapt() output is bit-identical to before this branch when relax is off (checksummed against the pre-change commit); and 41 tests pass.

@lmoresi

lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member Author

Staleness and review follow-ups

Merge hazard found and defused

This branch duplicates PR #420 (the per-cell collapse floor). It was cherry-picked from a branch that already carried that fix, while development still has the old median-based version.

The two copies were functionally identical but differed in em-dashes and a missing See #419., which would have produced a pointless comment-only conflict. The block is now byte-identical to #420 s, so:

Merge order: #420 first, then rebase this branch. The duplicate hunks drop out silently rather than conflicting. Merging this one first would leave #420 as a no-op against a fix it cannot see.

Staleness, for the record

The original feature/adaptivity-round3 branch is stale and should not be used. PR #409 was squash-merged on 2026-07-24, and development has since moved on via #416, #404 and #412 — a PR from that branch would have shown ~3445 lines of spurious deletions. All 14 commits were cherry-picked onto current development here and re-verified (built, 41 tests). That branch can be deleted; I have not done so, since it is not mine to remove.

Review item 5 addressed

The RBF retry now says it is dense. It rescues correctness when the orphan repair cannot place every coarse DOF, but returns nnz/row == n_coarse, making the Galerkin coarse operators dense — harmless on the meshes where it has fired, a silent performance cliff at production size, exactly the failure mode the GAMG fallback had. The warning now states that and says to fix the cause rather than accept the symptom.

Review items NOT addressed, and why

  • (1) Oracle gap. 1-2 vertices per 3D pass do not lie on a coarse edge and go unchecked by test_reproduces_an_arbitrary_coarse_field. Closing it needs a reference for points interior to a coarse face/cell — real work, and it is the gap most likely to hide a genuine error. Worth doing before this is relied on.
  • (4) Wall-clock unmeasured. The nested transfer is one iteration worse in 3D but ~3x sparser; whether it wins on time is unknown. It should not be merged as a claimed speedup.
  • (6) Cell-list engine paths still uncovered. Both its prolongation capture and its relax hook are unexercised, because both 2D and 3D prefer the native transform when built. Forcing that path in a test is straightforward (stub _nvbx) and has not been done.

23 tests re-run green after these changes.

lmoresi added 15 commits July 27, 2026 13:56
… 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
…h' 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
…-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
#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
…ustom-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
…ping 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
…eady 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
…ass (#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
…nsfer 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
)

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
…der (#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
…tion, 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
… 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
…CE 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
…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 force-pushed the feature/mesh-relax-and-nested-mg branch from 0581d8d to 380afcb Compare July 27, 2026 03:59
@lmoresi
lmoresi changed the base branch from development to bugfix/mmpde-graded-area-floor July 27, 2026 03:59
@lmoresi

lmoresi commented Jul 27, 2026

Copy link
Copy Markdown
Member Author

Restacked — the chain is no longer fragile

The duplication is gone, not merely aligned. This PR is now stacked properly:

#420  bugfix/mmpde-graded-area-floor -> development   [MERGEABLE]
#433  feature/mesh-relax-and-nested-mg -> #420        [MERGEABLE]

Verified: git diff #420...HEAD -- mmpde.py contains zero references to a0_own / _min_area_frac. The collapse-floor fix now exists in exactly one place — #420 — and this branch carries only its own changes.

What that buys:

Rebuilt and re-tested on the new base: 37 tests pass (relax, nested prolongation, graded adapt, custom-MG Stokes). Pre-restack state is preserved locally as backup/mesh-relax-pre-restack in case the history needs checking.

Remaining chain hygiene, not done by me

feature/adaptivity-round3 is dead and should be deleted — #409 was squash-merged on 2026-07-24 and development has moved on via #416/#404/#412. Everything from it lives here. I have left the deletion to you rather than removing a branch I do not own.

Still open from the adversarial review

Unchanged by the restack, and listed again so they are not lost:

  1. Oracle gap — 1-2 vertices per 3D pass lie on no coarse edge and go unchecked. The likeliest place for a real error to hide.
  2. Wall-clock unmeasured — nested is one iteration worse in 3D but ~3x sparser. Do not merge as a claimed speedup.
  3. Cell-list engine paths uncovered — both its capture and its relax hook; forcing that path needs _nvbx stubbed in a test.

@lmoresi
lmoresi deleted the branch development July 27, 2026 04:31
@lmoresi lmoresi closed this Jul 27, 2026
@lmoresi lmoresi reopened this Jul 27, 2026
@lmoresi
lmoresi changed the base branch from bugfix/mmpde-graded-area-floor to development July 27, 2026 04:31
@lmoresi
lmoresi merged commit cbf61d2 into development Jul 27, 2026
2 of 3 checks passed
@lmoresi
lmoresi deleted the feature/mesh-relax-and-nested-mg branch July 27, 2026 04:32
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