Anisotropic mesh-adaptation movers + mesh-owned boundary tangent-slip#228
Open
lmoresi wants to merge 34 commits into
Open
Anisotropic mesh-adaptation movers + mesh-owned boundary tangent-slip#228lmoresi wants to merge 34 commits into
lmoresi wants to merge 34 commits into
Conversation
…etric
Add the elliptic-MA mover and the arc-length monitor as options for the
two-knob follow_metric API, and unify the boundary slip across movers.
- follow_metric gains 'mover' ("anisotropic" default | "ma") and
'boundary_slip' (default True, used by the MA path). mover="ma" routes to
the Benamou-Froese-Oberman elliptic Monge-Ampere solver (_winslow_elliptic)
with n_outer=1 — one Caffarelli-clean convex-potential map, untangled by
construction, no Jacobi polish. The metric contrast (peak rho = refinement^d)
sets refinement; refinement <~ 3 is the capture-stable regime, higher is the
opt-in riskier path where the map stops tracking the metric.
- metric_density_from_gradient gains metric_choice="arc-length": a smooth
sqrt(1+(A*ghat)^2) monitor (no percentile-window kink, parameter-light),
capped at the refinement envelope. Pairs best with the MA mover — better
metric capture than the clipped front-following / gradient-uniform choices.
- Unify the boundary slip of BOTH movers (_winslow_equidistribute and
_winslow_elliptic) onto shared _ot_adapt._resolve_slip /
_build_slip_projector helpers using mesh.Gamma_P1 — radial-gated tangential
slide with on-surface snap-back; Cartesian boundaries pin. Replaces
_winslow_elliptic's old box/ring slip code.
- tests/test_0750: lock the arc-length envelope, MA+arc-length cleanliness
(signed-area untangled) + metric capture, MA boundary-slides-on-circle, and
the invalid-mover guard.
On the step-100 mode-1 field, follow_metric(mover="ma", metric="arc-length")
is untangled (0 inverted) and captures the metric better than front-following
(r~0.88 vs 0.78). Candidate clean default; the gradient-flow OT mover stays the
aggressive opt-in. follow_metric tests 14/14.
Underworld development team with AI support from Claude Code
…ility
The fault r-adapt work picked up from feature/elliptic-ma, end-to-end.
One coherent change set; tier-A 278 -> 295 passing, zero failures.
API additions (src/underworld3/meshing, exported via meshing/__init__.py)
--------------------------------------------------------------------
* fault_metric_tensor(mesh, faults, refinement=R, width=W) -- analytic
Eulerian normal-aligned anisotropic metric tensor M(x) for the
supplied-tensor mover; centres two close faults (gap 0.06) at
+/-0.0016 (beats the h-adapt reference -0.0037 by ~2.3x). Accepts a
Surface, polyline array, or a list mixing those.
* fault_comb_metric(mesh, faults, cell_size, n_across=4) -- scalar
"comb" density (teeth at d = k*cell_size from each fault) for the
isotropic MA mover. Equidistribution drops a node row at each tooth
-> evenly-spaced rows -> a roughly-uniform band of cell_size, with
the k=0 tooth pinning a row on the fault line (centring ~0.0002).
Handles curved/polyline faults via per-fault min-distance (bands
follow the curve).
* fault_metric(mesh, faults, method=..., cell_size=..., n_across=...) --
intent-based facade dispatching to the right representation per
method: scalar comb (ma), 2x2 tensor (anisotropic), h^-2 MeshVariable
(adapt / MMG). Documents that cell_size is *exact* for adapt (adds
nodes) and a *target* for the r-adapt methods (fixed budget).
* compose_metrics(metrics, compose="max") + smooth_mesh_interior now
accepts metric=[m1, m2, ...] or [(m, w), ...]. Internal weighted-max
on the excess (rho_combined = 1 + max(w_i * (rho_i - 1))). User just
hands over a list of feature metrics -- the routine composes them.
Mover changes
-------------
* _winslow_anisotropic (smoothing.py): supplied_D entry point on the
tensor mover -- explicit normal-aligned M used directly as D,
bypassing the grad-rho eigenframe derivation that mis-centres
codim-1 features. Cartesian box-face slip branch (geometric,
axis-aligned) for boundary-reaching faults. Validated gradient block
is byte-identical modulo indentation (verified via -w diff).
* _winslow_elliptic (smoothing.py): dimension-general -- now works in
3D. Three fixes, bit-identical at cdim==2:
1. Equidistribution normalisation c: was 2D-specific
1/<b^-1/2>^2; for the 3D simple-Picard source the leading term
is g (not 2 sqrt(g)) so c = 1/<b^-1>. With the wrong c the
source has nonzero mean, the pure-Neumann phi-Poisson is
unsolvable (DIVERGED_LINEAR_SOLVE -- the constant nullspace
fixes solution ambiguity, not RHS inconsistency), which was
the actual cause of 3D failure.
2. 3D source: was (g-1) - det(H), dropping the 2x2 principal
minors of det(I+H). Replaced (cdim!=2 branch) with the
dimension-general simple-Picard form
f_src = tr(Hs) + g - det(I + Hs)
(Hs symmetrised). Reduces to the old 2D else branch exactly.
3. Tet signed-volume backtrack: _tri_cells returns None for tets,
so 3D previously had no anti-tangle guard. Added _tet_cells +
_signed_volumes and a tet backtrack branch in the move.
Validated on a 3D slab and spherical-shell adapt -- refines toward
the feature, zero inverted tets. tier-A 290/0 (2D unchanged).
* smooth_mesh_interior: detects tensor-valued metric and routes to the
anisotropic supplied_D entry; new compose kwarg (default "max")
drives the list-of-metrics composition.
Boundary slip
-------------
* _ot_adapt: new _build_box_slip_projector helper (geometric
axis-aligned face slip for Cartesian boxes -- corners/edges pinned,
face nodes slide along the plane). Wired into _winslow_anisotropic
for the Cartesian case; _resolve_slip / _build_slip_projector
themselves are unchanged (the OT/MA shared projector stays
radial-only to preserve OT_adapt's validated pinned-box behaviour).
Surface fix (carried from prior session)
----------------------------------------
* Surface._symbol default name[0].upper() -> name (full, unique).
Two faults sharing a first letter were silently aliasing onto the
same distance varsymbol; the analytic Eulerian metric path needed
distinct symbols.
Tests
-----
tests/test_0762_fault_metric_tensor.py -- 17 tier-A locking the new
API: tensor structure / Surface-equivalence / centring / non-2D
raise; comb teeth structure / two-fault band / curved (arc) /
non-2D raise; facade ma==comb / anisotropic-is-tensor / adapt-h^2
field / unknown-method-raise; compose passthrough / equal-weights ==
sympy.Max / weighted-excess scaling / tensor-rejection /
list-of-metrics through smooth_mesh_interior. 17/17 pass.
Underworld development team with AI support from Claude Code
Replace the radial-only Gamma_P1 slip and the geometric box-face slip with ONE topology-based projector that works on any 2D/3D simplicial mesh (Cartesian box, annulus, sphere, polyhedron, curved surface). Root cause of the previous Gamma_P1 path failing on Cartesian: _update_projected_normals (discretisation_mesh.py:1927) evaluates the PETSc petsc_n quadrature symbol at vertices, which is undefined off boundary integration. Radial mesh classes redefine Gamma as an analytic radial unit vector (vertex-evaluable) and so the path appeared to work, but Cartesian got garbage normals -- which is why the old code gated slip to radial-only and pinned every Cartesian boundary. New helpers (_ot_adapt.py) -------------------------- * _boundary_facets(mesh, cdim) -- per-cell facets (2D edges / 3D triangles) hashed; a boundary facet is one that occurs in exactly ONE cell. Returns the facet vertex tuples plus the opposite cell vertex (used to orient the outward normal). * _boundary_vertex_normals(mesh, parallel_tol_deg=15.0) -- per-facet outward normal computed GEOMETRICALLY from the cell coordinates and oriented away from the opposite vertex; area-weighted averaged per boundary vertex. A vertex is classified "face-slip" iff every incident facet normal lies within parallel_tol_deg of the average (one smooth face / single tangent plane); otherwise it is pinned (corner, or 3D edge between two faces). Verified directly: box -> 60 face-slip + 4 corners pinned, normals exactly (-1,0)/(1,0)/(0,-1)/(0,1); annulus -> 96/96 face-slip, normals match r-hat to 5e-10. Wiring ------ * _resolve_slip: drop the radial-only gate. Generic slip works for any geometry; the helper now just maps the boundary_slip argument to a bool. * _build_slip_projector: uses _boundary_vertex_normals. Tangential slide at face-slip nodes, pin at corners/edges; radial snap-back layered on top for curved boundaries (kept). * _winslow_anisotropic: replaces its inline ring + geometric box-face branches with a call to the shared _build_slip_projector. One slip path across _winslow_anisotropic, _winslow_elliptic, _winslow_equidistribute, and _ot_adapt. * Dead helpers removed: _slip_normals (Gamma_P1-based, supplanted by topology normals) and _build_box_slip_projector (geometric, now subsumed by the generic projector). Behaviour change ---------------- OT_adapt on a Cartesian box now SLIDES boundary face nodes tangentially (it used to pin every Cartesian boundary node). Corners stay pinned, box shape is preserved (boundary nodes stay on the bounding planes). Updated test_ot_adapt_box_moves_interior_pins_boundary -> test_ot_adapt_box_moves_interior_slides_boundary_on_faces to assert the new (and intended) behaviour: corners present, every boundary node on an axis-aligned bounding plane, FE-remap pattern preserved (tolerance loosened 5e-2 -> 1.5e-1 -- face slip mildly degrades the remap of a sharp tanh feature at coarse base mesh, still well- controlled relative to the BL amplitude of 1). Validation ---------- * fault_boundary_slip.py harness on a full-span box-crossing fault gives the SAME boundary refinement as the previous geometric box-slip (min boundary node gap 0.0200 -> 0.0044, faces flat, corners pinned, 0 inverted). * The annulus radial-slip test (test_follow_metric_ma_boundary_slides_on_circle) is unchanged. * tier-A 295 passed / 0 failed. Underworld development team with AI support from Claude Code
Captures the convergence after the fault-meshing work (c977e48 + f2d4050): one mover (single-shot Monge-Ampere), one metric form (scalar density), one composition (weighted max on excess), one slip (topology-based vertex normals). Works in 2D and 3D, on Cartesian boxes, annulus, sphere, polyhedra, curved surfaces. The note covers the recipe, why each piece, what this collapses (the anisotropic tensor mover, fault_metric_tensor, the per-segment analytic machinery, the ring-projection / geometric box-slip branches), the honest limits (budget cap, multi-feature composition, n_outer>1 metric convection, 3D MA simple-Picard fragility on strong metrics), and a migration table for users of the deprecated paths. Underworld development team with AI support from Claude Code
Side-by-side comparison (fault_compose_demo2.py, cs=1/100, only boundary_slip flipped) shows that for any composed case where a feature TOUCHES the boundary -- thermal BL running full width, fault reaching the wall -- pinning the boundary wastes the budget at the edges: the refined band visibly fades as it approaches the wall. With the generic topology slip (f2d4050), boundary face nodes slide along the face to cluster where the metric demands them, and refinement runs uniformly to the wall. Corners pinned, box shape preserved. So boundary_slip=True is the default for the recommended recipe, not an optional tweak. Underworld development team with AI support from Claude Code
…de_rho note Pinning down when n_outer>1 helps and when it breaks down. The MA mover re-queries the target metric on the deformed mesh between outer iters. For Eulerian metrics (analytic sympy of X -- the comb from segments, the anisotropic supplied D) this gives sharper, more on-line refinement and n_outer=2 + target_side_rho=True is the recommended setting. For Lagrangian field-backed metrics (metric_density_from_gradient on a MeshVariable T, the comb on a Surface.distance field) the underlying values ride with the mesh -- the feature itself convects, the mover chases a moving target, bands smear instead of tightening. Stick with n_outer=1 in that case. Adds the safety table for which metric is safe at n_outer>1. Underworld development team with AI support from Claude Code
…onestly
After instrumenting the Picard loop with a per-iter probe, target_side_rho
turns out to be the wrong lever to recommend:
- The Picard fixed-point coupling between phi and gradphi is much tighter
than the source-side iteration. The default n_picard=25 is fine for
source-side (converges in ~13 iters) but wildly under-converged for
target-side -- needs ~100+ iters for moderate demand, ~200 for strong.
At n_picard=25 with target_side_rho=True the iteration is cut off
mid-transient: |gradφ|_max is still climbing monotonically, |Δφ| is
many orders of magnitude above the convergence floor, the resulting
displacements are larger than the converged solution and the
backtrack scales them down.
- Even fully converged (n_picard=200), target_side_rho does NOT produce
sharper realised refinement than iterating source-side with multiple
smooth_mesh_interior calls.
- The mechanism IS implemented correctly (verified via the two-pass
diagnostic: target_side_rho pass-2 motion is ~10% of pass-1 because
the fixed point has been reached; source-side pass-2 motion equals
pass-1 because each call ignores prior adaptation). But the
user-facing cost/benefit is bad.
n_outer characterised honestly:
- n_outer>1 is the patch-volume-aware MA composition: each outer step
divides the metric by the current deformed mesh's per-vertex patch
volume, so it's CONSERVATIVE about already-achieved refinement
("the mesh is partly adapted; don't over-pull"). It is NOT
"compose maps until convergent" in the naive sense -- on a comb
metric, n_outer=2 plateaus at far/band ratio ~2.6x and doesn't move
further at n_outer=4.
- Calling smooth_mesh_interior repeatedly with n_outer=1 each time
is more aggressive (each call sees patch=ones and re-pulls as
though uniform) and reaches ~3.7x far/band ratio for the same
metric. Not strictly the Caffarelli optimum, but stable and
predictable.
Both behaviours are useful; the design note now says so plainly and
stops misdirecting users to target_side_rho.
Underworld development team with AI support from Claude Code
In _winslow_elliptic, the source-side density V at each vertex was
either patch=ones (at n_outer=1, the default) or the lumped mass
diagonal M^lumped_ii = Σ |T| / 3 used as a density (at n_outer>1).
Both were wrong, in different ways:
* patch=ones assumed the input mesh is uniform regardless of how
it actually was, so repeated calls didn't compose. Calling the
mover a second time produced the same displacement the first
call would have produced from cold, applied on top of the
existing deformation — biases compounded across calls instead
of correcting.
* M^lumped_ii is an integral (∫ ψ_i dx, units of area), not a
density. On unstructured Delaunay with valence varying 5..7,
using it as V_i added a ~30% spurious source non-uniformity
from FE bookkeeping with no relation to actual mesh deformation.
Replace with the lumped L2 projection of the cell-wise V_T = |T|
into the P1 vol_field:
V_i = (Σ_{T ∋ i} |T|² / k) / (Σ_{T ∋ i} |T| / k)
= Σ |T|² / Σ |T|
— the area-weighted average of incident cell volumes. Strictly local,
no neighbour mixing (zero kernel scale), valence-independent on
uniform meshes (= |T_0| exactly when all |T| equal). Same path in
2D and 3D, with the cell corner count k branching automatically.
Effect: the mover is now properly composable. Repeated calls produce
displacements from the actual current mesh state toward the fixed
point, |Δo| decreases monotonically, no compounding.
This is what unlocks the pre-placement + redistribution recipe (a
wide max-of-Gaussians metric clusters cells around the whole feature
system, then a narrow stage localizes them onto individual peaks).
For two close faults at gap 0.060 with target band width 0.015, that
recipe reaches f0/f1 = -0.0005/-0.0014 (bands within 1/12 of one mesh
cell of the actual lines), where cold-narrow-from-uniform settles at
the centroid-bias floor of f0/f1 = -0.0109/+0.0103.
An intermediate attempt used the consistent-mass uw.systems.Projection
to project V_T → vol_field; that introduced an intrinsic L2 kernel
of ~one element which smoothed cell-density signals narrower than the
kernel into a halo and made iteration regressive (refinement undone
on the second pass). The lumped form has no kernel scale.
Tier-A 64/0 on level_1 mesh/fault subset. Design doc updated with the
diagnosis, the fix, the pre-placement recipe, and a width-vs-separation
sweet-spot table.
TODO: the np.add.at accumulation is rank-local. At MPI partition
boundaries, vertices owned by one rank under-count off-rank incident
cells. Same parallel deficit as the old _patch_volumes had. Fix is to
assemble num/den into PETSc Vecs with ADD_VALUES so assembly handles
ghost summation. Required before parallel use on adapted meshes.
Underworld development team with AI support from Claude Code
Diagnostic experiment (diagnose_convergence.py): same MAX recipe
(w=0.120 ×4 → w=0.015 ×8) run at n_picard ∈ {25, 50, 200}, tracking
cv(ρV) = std(ρ·|T|)/mean(ρ·|T|) per outer call as the equation-natural
equidistribution residual.
* n_picard=50 and n_picard=200 are bit-identical (inner Picard fully
converged at 50), and the recipe locks into the centroid local
minimum at cv(ρV) = 1.07 and never escapes. Bands stuck at f0/f1
= -0.0041/+0.0060.
* n_picard=25 escapes the centroid local minimum at outer iter 7
via numerical annealing: the under-converged inner Picard leaves
a residual perturbation each call, and accumulated perturbations
eventually kick the system into a deeper minimum at cv(ρV) = 0.79.
Bands settle at f0/f1 = -0.0005/-0.0014 (on the lines, ≤ 1/12
mesh cell).
The "outer iterations" weren't compensating for inner non-convergence
in the standard sense — they were *using* it to escape shallow local
minima. Tightening n_picard locks the recipe into the wrong basin.
Practical stopping condition: cv(ρV) plateau after evidence of a
significant drop. Document the rule and call out that n_picard=25 is
deliberate.
Geometric |Δo| reads ≈ 0 in the centroid minimum — useless as a
stopping signal there. cv(ρV) distinguishes the two basins.
Underworld development team with AI support from Claude Code
This session's exploration superseded the earlier "n_picard=25 +
Anderson" framing. The new section reflects the validated current
state:
* The mover catastrophically misses one of two close faults from
cold-start without a low-amplitude smooth Gaussian added on top
of the sharp narrow one. The smooth aid provides non-trivial
∇ρ everywhere in the fault neighbourhood — essential.
* Plain Picard is geometrically equivariant. Anderson finds a
deeper basin per-call (3× faster) on a fixed geometry but a
*different* basin on a translated geometry: a uniform shift
of all faults that should give a uniformly shifted solution
does not under Anderson. cv at the same geometry can be
0.39 (Anderson, deep) or 1.08 (Anderson on shifted, stuck).
Plain Picard hits the same cv ≈ 0.57 regardless of fault
position — recipe is the reliable default.
* |ΔX|/(√N · h) is the natural fixed-point convergence signal.
cv(ρV) is a quality measure (lower = deeper basin) and
distinguishes which basin you landed in but is NOT a
convergence test.
* For moving-fault use case: build a fat refined band wide
enough (~ v_fault × Δt_remesh) to contain the fault over
multiple timesteps, re-mesh only when fault exits.
~3× speedup vs re-meshing every step. Warm-start from prior
timestep's mesh does NOT track — cells inherit old fixed
point and plain Picard from that state finds suboptimal
basin, not the moving one.
* Next major efficiency lever: SNES wrap with approximate
Jacobian on F(X) = X - mover(X). Folding mesh deformation
into a Newton step inside the linear solve is the natural
next move. Left as a follow-up session.
The earlier sections (n_outer single-shot, scalar comb metric,
composable list, generic slip, dimension-general MA, lumped V_T
projection) remain valid — this update adds the user-facing recipe
on top.
Underworld development team with AI support from Claude Code
…sign doc Finalises the variational anisotropic moving-mesh adaptation (Huang-Kamenski MMPDE, JCP 301 (2015) 322 / arXiv:1410.7872) landed as smooth_mesh_interior(metric=M_tensor, method="mmpde"). - Document the per-node step cap as the sole serial/parallel divergence (~0.006%): the velocity assembly via localToGlobal(ADD_VALUES) is bit-identical serial-vs-parallel; only the rank-local min-incident-edge cap differs at partition boundaries. Accepted as a known small non-reproducibility (below the move's own non-determinism); the ghost-complete MIN reduction route is unavailable (PETSc localToGlobal MAX_VALUES errors on the coord DM). - Add the design doc docs/developer/design/anisotropic-mmpde-mover.md: dimension-general (2D/3D) formulation with q=d*p/2, the essential dG/dM metric-variation term, energy line-search backtrack, boundary slip, tunability/cap characterisation, PETSc-parallel-safe port plan, and the non-folding rationale. Validated: 2D band 0.27 on-fault 0 crushed; 3D on-plane 1.00 0 crushed (serial + np2); uniform metric exact no-op; composable (settles, never explodes); tunable (across-ratio, p); grad-T composition by tensor addition. MA and spring movers unaffected. Underworld development team with AI support from Claude Code
The metric is a guide field (not a solved quantity), so exact FE/PETSc interpolation is unnecessary. Bake the analytic metric ONCE onto the fixed reference cloud, then Shepard/k-NN interpolate to the moving centroids each step (metric_eval="rbf", default). 80-step annulus fault adapt: 72.9s -> 19.8s; cross-fault refinement 2.5x -> 2.3x (unchanged), 0 crushed, min-angle 18->20 (RBF smoothing of the analytic endpoint "elbow" kink helps). metric_eval="analytic" retains the exact path. KNOWN PARALLEL LIMITATION (serial unaffected): the RBF reference cloud is rank-local; a mesh node that drifts past the halo's spatial coverage gets nearest-neighbours from distant on-rank points (silently wrong metric at partition boundaries). The robust fix is a migrating SWARM for the Eulerian metric cloud (fixed coordinates, ownership migrates to track the moving decomposition) - flagged for parallel hardening. Underworld development team with AI support from Claude Code
… bug Harden the metric-driven adaptivity strategy: - Lock to SIMPLEX meshes. Non-simplex (quad/hex) meshes previously hit a silent no-op (_tri_cells/_tet_cells returned None, mover returned early). Now _assert_simplex_for_adaptivity() raises a clear NotImplementedError at the smooth_mesh_interior front door for any metric-driven call; _winslow_mmpde also raises defensively. Simplex test: first cell cones to cdim+1 facets (3 edges / 4 faces). - Add _simplex_cells(dm, cdim) — single dimension-general entry point (triangles in 2D, tetrahedra in 3D) so connectivity consumers don't each branch on cdim. - FIX latent 3D bug: _winslow_mmpde used _signed_areas (z-projection) for the tet fold-detection backtrack instead of _signed_volumes. Earlier 3D tests passed only because they were non-folding cases; now corrected. Verified: 3D simplex 0 crushed, quad raises, 2D unaffected. Underworld development team with AI support from Claude Code
The variational MMPDE mover is demonstrably the most capable and most straightforward of the metric movers, so make it the default for smooth_mesh_interior (method="mmpde", was "spring"). Safe: the method dispatch only runs when metric is not None, and every internal metric-bearing caller specifies method= explicitly; the metric=None graph-Laplacian Jacobi path is unchanged. Verified all movers still dispatch (ma/ot/anisotropic/spring) and a bare metric call now runs mmpde; scalar metric -> isotropic tensor, d×d tensor -> supplied path. Also folds in the mover work from this session (previously uncommitted): - tol: scale-relative convergence exit (dmax < tol·h0), so adapts exit on residual not the n_outer cap; the old absolute outer_tol never fired. - metric_eval="rbf" (default): bake the analytic metric once, Shepard- interpolate during iteration — ~3.7x faster, same mesh, smooths the analytic endpoint kink. metric_eval="analytic" keeps the exact path. - area_floor_frac, pernode_backtrack knobs. - simplex lock-out (_assert_simplex_for_adaptivity) + dimension-general _simplex_cells; fixed a latent 3D bug (signed_areas->signed_volumes in the tet fold-detection backtrack). Docs: docstring method= section now leads with "mmpde" (default) + paper refs (JCP 301 (2015) 322; Math. Comp. 87 (2018) 1887); docs/advanced/ mesh-adaptation.md updated to present mmpde as default with scalar/tensor guidance and the legacy movers as opt-in. Known follow-ups (next): named-surface tangent slip via Gamma_P1 (_build_slip_projector still all-or-nothing); resume-from-deformed-mesh quality; optional ALE single-interpolation variant. Underworld development team with AI support from Claude Code
Replace the all-or-nothing boolean boundary_slip with a per-named-surface
slip mechanism shared by every metric mover via _build_slip_projector.
Key fix: slip-vs-pin is now LABEL-DRIVEN, not normal-agreement. A boundary
vertex slips iff it lies on exactly one slip surface; vertices on a non-slip
boundary or at a junction of two slip surfaces (e.g. a box corner, where the
normal is ambiguous) pin. The old topology classifier used a 15-degree
incident-facet-agreement test that spuriously pinned the coarse-but-smooth
annulus inner ring (a low-resolution polygon's adjacent facet normals diverge
>15 degrees yet it is no corner) — inner ring now slips 28/28 (was 4/12),
while box corners correctly pin (4 pinned, edges slide along their line).
- Tangential slide uses the projected P1 boundary normal (mesh.Gamma_P1),
smooth and consistently oriented on curved boundaries.
- Return-to-bounds has two modes: an exact analytic |r| snap for radial
geometries (annulus/sphere/cylinder), which is concave-safe; and a
geometry-general nearest-reference-facet snap as the fallback for
non-analytic surfaces. Free surfaces (dict value False) skip snapping.
- slip_surfaces accepts True / label / [labels] / {label: snap_bool};
boundary_slip is kept as a deprecated alias. The resolved spec is threaded
to the movers unchanged so dict no-snap survives, and Gamma_P1 is
pre-touched before any solver DM is built (footgun guard).
A TODO(watch) marks the known concave-non-analytic-surface bias of the facet
snap (radial geometries are immune via the analytic branch); the deferred
cure is a smoothness / mean-preservation constraint.
tests/test_0855_slip_surfaces.py locks the API (tier_a, level_1).
Underworld development team with AI support from Claude Code
The slip rewrite builds tangent normals from mesh.Gamma_P1, which lazily creates the _n_proj MeshVariable. When that creation happened inside the mover (via _resolve_slip), it restructured the DM and invalidated the JIT/interpolation state needed to evaluate a MeshVariable-valued metric — a hard, uncatchable abort. Pure-sympy metrics never tripped it (they touch no DM interpolation), so it surfaced only when composing a real density field (e.g. metric_density_from_gradient + fault_metric). smooth_mesh_interior now pre-creates Gamma_P1 once at the top, before dispatching to the mover and before any metric is evaluated; the in-mover _resolve_slip touch then finds it already built. Regression test added (mmpde + slip + MeshVariable metric). Underworld development team with AI support from Claude Code
Add a generic opt-in _pre_solve_hook on the base SNES solver and a utilities/gmg_geometric_interpolation module that re-targets the finest multigrid prolongation to current node positions each setup (recompute- nested-values: overwrite the existing interp Mat values in place with the minimal correction that reproduces the moved positions). Keeps velocity-block GMG iteration-flat on mover-adapted meshes without mutating mesh.dm. Underworld development team with AI support from Claude Code
…t-in stol stagnation exit A degenerate/near-inverted cell gives an inf gradient; the per-node cap then makes step = inf*0 = NaN, which produces a NaN trial whose centroid kd-tree query crashes _energy on a subset of ranks and deadlocks the parallel job. Zero non-finite steps at source + reject non-finite trials in the line-search. Also add an opt-in stol residual-stagnation exit to _winslow_mmpde. Underworld development team with AI support from Claude Code
…ries by default
`Mesh._test_if_points_in_cells_internal` used a strict `> 0` test on the
squared-distance difference between mirrored inner/outer control points
placed ±1e-3 along each face normal. A query exactly on a cell face has
zero distance difference and was rejected. That meant mesh vertices —
which sit on the faces of every cell containing them in their closure —
failed the in-cell test for every candidate cell, and
`_get_closest_local_cells_internal` returned -1 for them.
For evaluation use cases this is the wrong semantics: a closed domain Ω
includes ∂Ω, and the FE basis at a shared vertex / face is consistent
across the adjacent cells (a DM consistency requirement of FE assembly).
So "any cell whose closure contains the point" is the right notion.
Routing those queries through RBF Shepard extrapolation (the legacy
behaviour of `uw.function.evaluate` for points classified outside the
domain) is less accurate than evaluating via FE on any adjacent cell.
Add `on_boundary: bool = True` kwarg to the in-cell test and forward
through the call chain:
- `_test_if_points_in_cells_internal(on_boundary=True)` — core test
- `_get_closest_local_cells_internal(on_boundary=True)` — forwards
- `test_if_points_in_cells(on_boundary=True)` — public wrapper
- `get_closest_local_cells(on_boundary=True)` — public wrapper
With on_boundary=True (the default) the comparison is `>= -1e-12` (well
below the 1e-3 control-point offset, well above 64-bit float roundoff);
with on_boundary=False the historical `> 0` strict-inside semantics is
preserved. Callers that need uniqueness (a point claimed by exactly one
cell, never a shared-face point claimed by both adjacent cells) can opt
back into strict.
Observable consequences:
- `_get_closest_local_cells_internal` is now an authoritative cell-id
source: returns a cell whose closure contains the query, or -1 if no
local cell does.
- `points_in_domain` reports boundary-vertex queries as in the domain
(matches mathematical convention; a point on ∂Ω is in Ω).
- `uw.function.evaluate` now routes boundary-vertex queries through FE
instead of RBF Shepard — the more accurate path. This shifts the
output of `evaluate` at mesh vertices on the domain boundary.
Test impact:
- `tests/test_0820_in_cell_test_loose_semantics.py` (new, 7 tests):
locks the on_boundary=True default — every vertex of 2D simplex, 3D
simplex, and 2D quad meshes resolves to a containing cell; strict
mode still rejects on-face queries; loose mode returns cells whose
closure genuinely contains the query.
- `tests/test_0820_deform_mesh_solver_rebuild_regression.py` continues
to pass (it was the motivating regression for the PR #203 bypass
work that depends on this).
- `tests/test_1100_AdvDiffCartesian.py::test_advDiff_boxmesh[mesh0]`
marked xfail (strict=True). The file header already calls this "not
a great test"; its atol=0.05 tolerance previously aligned only
because boundary-vertex queries went through RBF Shepard's
smoothing, and the test was tuned to that result. Needs reworking
(smoother IC, larger transport distance) to pass under the more
accurate FE-evaluate path. The two simplex variants of the same
test (mesh1, mesh2) continue to pass.
Also addresses Copilot review feedback:
- docstring on `_test_if_points_in_cells_internal` corrected to say
model-coords input (the public wrapper handles units).
- removed an unused `inside` initialisation from the refactor.
- public `test_if_points_in_cells` now coerces `cells` to a 1-D numpy
array so list/tuple input works as the docstring promises.
This unblocks the bypass design in PR #203
(`feature/dminterp-bypass-element-check`), which needs an authoritative
per-rank cell-id source. With the new default,
`mesh._get_closest_local_cells_internal(coords)` directly gives the cell
id whose closure contains each query — exactly what the bypass requires.
Underworld development team with AI support from Claude Code (claude.com/claude-code)
The mover-adapted parallel (np>1) convection runs grew spurious "cold finger" temperature artifacts at partition seams that never appeared in serial. Root cause: this branch (elliptic-ma base) forked from development before PR #213 ("parallel adaptive seam-spike fix") and so lacked the global_evaluate field-transfer + robust parallel locator that carry fields onto mover-moved nodes correctly across rank boundaries. A fixed-dt serial-vs-parallel diff confirmed it: on a fixed mesh serial==parallel to roundoff; with adaptation on, pre-merge bulk rms divergence ~0.04 (widespread ±0.3 dipoles); post-merge ~0.007 (the residual is mesh-discretisation noise between two valid adapted meshes). Resolution (two divergent mover lineages): - Took development's discretisation_mesh.py / _ot_adapt.py / test wholesale (locator hardening, global_evaluate remap) — our on_boundary cherry-pick is subsumed. - Took development's parallel-correct smooth_mesh_interior architecture (deform-back/global_evaluate/deform-forward field-transfer wrap) and grafted our mmpde mover (+3 helpers), the arc-length metric_choice, and the slip_surfaces alias + Gamma_P1 guard back onto it. - Grafted 8 slip helpers into _ot_adapt.py that the mmpde mover imports. Serial performance fix (_dminterp_wrapper.pyx): - Development gated the DMLocatePoints-bypass hint behind _eval_use_robust_location() (parallel-only), which regressed serial FE evaluate ~7x (advect 8s->58s, adapt 6.5s->38s). The bypass is safe whenever the hint is authoritative (simplex/manifold) independent of rank count, so gate use_hint on that instead. Restores serial speed (advect 2.8s, adapt 7.8s); the hint *source* still follows the parallel robust path. (PR #203 fast trace-back, serial.) Underworld development team with AI support from Claude Code
_winslow_mmpde required a full d×d tensor (sympy.Matrix(metric) on a scalar raised TypeError). The ma/ot/anisotropic movers all accept a scalar density, so mmpde now coerces a scalar expression or a 1×1 (scalar-MeshVariable) metric to rho*I. Fixes test_0855 test_mmpde_slip_with_meshvariable_metric (passes metric=f.sym[0]). Underworld development team with AI support from Claude Code
After merging development's parallel-correct mover architecture (keeping mmpde + the arc-length metric_choice, dropping the elliptic-ma MA mover / follow_metric mover= / list-compose), six tier-A tests exercised dropped or moved capabilities. Aligned them with the supported API: - test_0762 centres_two_close_faults: tensor metrics are moved by mmpde (the production fault path), not dev's scalar-density anisotropic mover. Switched method="anisotropic" -> "mmpde"; relaxed the band-count proxy 20 -> 15 (mmpde concentrates ~19 vs the anisotropic mover's ~20; the centring assertion — the rigorous check — is unchanged and passes). - test_0750 follow_metric tests: dropped the removed mover=/boundary_slip= kwargs; arc-length survives as metric="arc-length". invalid_mover_raises -> invalid_metric_raises. The two MA-mover-specific behaviours (strong arclength alignment>0.6, boundary slip) are xfail(strict=False) with pointers to where the capability now lives (OT mover; mmpde slip_surfaces). - test_0762 list_of_metrics: xfail(strict=False) — list composition dropped; compose via fault_metric_tensor / fault_comb_metric. Full tier-A (level_1): 200 passed, 3 xfailed, 0 failed. Underworld development team with AI support from Claude Code
The geometry-aware multigrid prolongation experiment (commit 486d0ae) was parked when we settled on the default Galerkin-nested velocity multigrid for mover-adapted meshes (see project_stokes_gmg_velocity_block). Remove the unused machinery: - src/underworld3/utilities/gmg_geometric_interpolation.py (the parked custom-interpolator build; injection was never unblocked) - the opt-in `_pre_solve_hook` seam in petsc_generic_snes_solvers.pyx, whose only setter was that file. The hook was a no-op for every shipped solver, so removal is behaviour-neutral. The work remains in history (486d0ae) if the geometric-interpolation route is revisited. Underworld development team with AI support from Claude Code
…slip leak) On a freshly mover-adapted annulus run in parallel, the free-slip penalty BC appeared to "leak": nodal v.n at the walls reached ~100 (vrms ~50) and advection then deposited hot T on the cold Dirichlet wall. The 2x2 factorial pinned it to parallel x adapt x boundary-slip. Root cause is NOT the BC and NOT DM-state corruption (coords, boundary labels, facet normals and the point SF are all bit-clean on the leaking mesh; a cold-start Stokes solve on the final deformed mesh is always clean -- post-solve nodal v.n ~1.5). The corruption is purely in the adapt field-transfer: `_remap_var_set` evaluates the OLD P2/P3 field at the NEW boundary DOF coordinates, which after tangential slip sit a sagitta OUTSIDE the OLD boundary cell (arc vs chord). In parallel the swarm migrate lands those points in a containing cell on another rank, where the FE basis OVERSHOOTS -- and the result is delivered as a "valid" (un-flagged) value. It is the same FE-overshoot family as the SemiLagrangian trace-back bug. Fix: evaluate the transfer with `monotone="clamp"`, bounding each resampled value to its k-NN source-nodal range. This is bit-identical to plain FE in smooth regions and parallel-safe (rank-local), and it kills the overshoot. Validation (np=5 fault repro): wall nodal max|v.n| 118 -> 3.0, hot-wall node count 171 -> 0. Isolated deterministic-rotation transfer gate: serial boundary error bit-identical with/without the clamp (no serial regression); parallel boundary error 1.36 -> 0.045. tier-A (level_1) 200 passed; parallel swarm tests green. Underworld development team with AI support from Claude Code
Lets the monotone clamp in _remap_var_set be switched off for diagnostics (default stays 'clamp'). Keeps the working tree clean across session resumes during the parallel-transfer bisection. Underworld development team with AI support from Claude Code
_winslow_mmpde is first-order steepest descent of Huang's functional (velocity
= -grad, with the energy/min-area line-search as a fold guard). On a stiff
radial-equidistribution map it converges only linearly and the first adapt hits
the n_outer cap. Add an opt-in acceleration of the search DIRECTION; the
line-search stays the guard, so acceleration overshoot is backtracked and never
tangles (verified fold-proof even at step_frac=2):
MMPDE_ACCEL=heavyball step += beta * prev accepted displacement (Polyak;
beta from MMPDE_MOMENTUM)
MMPDE_ACCEL=hb-restart heavyball + gradient restart
MMPDE_ACCEL=cg nonlinear conjugate gradient (Polak-Ribiere+),
parameter-free
First radial adapt: 202 -> 64 iters (heavy-ball 0.99), 202 -> 15 (CG); a
cold->developed-plume adapt: 51 -> 14 (CG). CG gives the best mesh quality and
needs no tuning, making adapt-every-step affordable. Default off (none).
Underworld development team with AI support from Claude Code
Harness updates for the adaptive-convection mover work: - MOVER env dispatch (mmpde [default] / anisotropic / ma / ot) + MOVER_SLIP (ring) + MMPDE_ACCEL pass-through; --res resolution knob. - FMG/GMG velocity block (REFINE/PCVEL/MG_TYPE) on the Annulus dm_hierarchy. - median + direction-aware dt (--dt-cell-percentile) so sliver cells don't collapse dt; collective Tmax overshoot guard + collective vrms (were rank-local). - per-step run_log.txt in the run dir: Nu/vrms/solver-iterations + per-phase wall (Stokes/advection/adaptation) + metric mismatch; _RUN_DONE sentinel for the live render watcher. Depends on estimate_dt(percentile=) and the swarm empty-partition fix (PRs to development) — merge development in once those land. Underworld development team with AI support from Claude Code
…M env) The nonlinear-CG / heavy-ball acceleration of the steepest-descent direction in _winslow_mmpde was selected by the MMPDE_ACCEL / MMPDE_MOMENTUM environment variables. Promote them to real function kwargs (accel="cg", momentum=0.0), removing the os.environ reads from the library. accel="cg" (parameter-free nonlinear conjugate gradient) is the default — the production choice. Unknown accel values now raise ValueError. The stagnant_lid_adapt_loop harness forwards accel/momentum through method_kwargs (script-level env reads are fine; the library no longer reads the environment). Underworld development team with AI support from Claude Code
`refinement` (R) was accepted but unused — a silent no-op; grading was set only by amp/power, so all earlier R=1.4/3/5 comparisons used an identical metric. Make R the target COARSEST:FINEST edge-length ratio: the mover equidistributes rho=(1+amp*t)^power so cell edge h ∝ rho^(-1/d), giving h_max/h_min = (1+amp)^(power/d) over t∈[0,1]. Invert to set amp from R, so R is a predictable, mesh-independent resolution knob (overrides strategy amp). The stagnant_lid_adapt_loop harness gains --resolution-ratio R to drive the metric through this live knob while keeping the MMPDE mover. Underworld development team with AI support from Claude Code
…c-mover-adapt-transfer
Adds the `facet` kind to BoundingSurface: nearest-point restore onto reference facets (segments in 2D, triangles in 3D) via _ot_adapt._nearest_on_facets_*. This is the piece that will let mesh.boundary_slip reach full parity with _build_slip_projector for non-analytic boundaries (instead of pinning them). WIP — not yet wired: mesh.boundary_slip still pins unregistered labels; the transient-facet-surface construction, the call-site swap, and the harness A/B validation are the remaining step-2 work. Additive and inert for now (nothing calls facet restore yet; test_0762 11/11 still green). Underworld development team with AI support from Claude Code
…lip engine (step 2) Swap all five metric movers (mmpde, spring, ma, ot, anisotropic) from their private / inline boundary tangent-slip onto the mesh-owned mesh.boundary_slip contract (step 1, PR #225) and delete the duplicates. Net -317 lines. mesh.boundary_slip: - Facet fallback: a slip label with no registered analytic surface now builds a transient `facet` BoundingSurface from the reference facets (was: pinned). - Normals computed ONCE per build (not per project() call) to avoid a Gamma_P1 re-solve inside the movers' line-search backtrack; a DESIGN NOTE on BoundingSurface.normals records the re-solve-vs-cached trade-off + escape hatch. - Projects only OWNED slip vertices (parallel safety; ghosts via the movers' halo-sync). BoundingSurface rejects degenerate plane normals / bad radii. Movers: - mmpde / spring / ma / ot / anisotropic call mesh.boundary_slip(...) with the current coords as reference. spring/ma/anisotropic drop their per-ring rotation anchor (mmpde never had one; signed-area guards prevent tangle — only cosmetic azimuthal re-parameterisation). OT keeps its radial-only slip gating and now emits a DeprecationWarning (incomplete; superseded by mmpde + a scalar metric). Deletions: _build_slip_projector, _gamma_p1_at_vertices, _label_vertex_mask, _boundary_centre (no remaining callers). Validation: parity vs the deleted engine was machine precision (1.57e-16 annulus, 0.0 box); convection-harness A/B (serial 40-step + np=5) physics unchanged; 22 slip tests green. Documented follow-up (inline TODO): _pinned_mask misses 3D face-only labels (markVertices=False) -- a shared-helper limitation. Tests: test_0855 migrated to mesh.boundary_slip; test_0762 covers the facet fallback + degenerate-geometry rejection; new test_0763_boundary_slip_correctness locks exact-on-surface landing. Also fixes a latent missing `import warnings`. Underworld development team with AI support from Claude Code
…ropic-mover-adapt-transfer # Conflicts: # src/underworld3/discretisation/discretisation_mesh.py # src/underworld3/meshing/bounding_surface.py
Contributor
There was a problem hiding this comment.
Pull request overview
This PR integrates validated anisotropic mesh-adaptation work and hardens parallel correctness by (1) making the variational MMPDE-based mover the default path for adaptive node movement, (2) refactoring boundary tangential slip so movers consume the mesh-owned mesh.boundary_slip(...) contract (with a facet fallback), and (3) fixing a parallel-only adapt field-transfer overshoot by enabling monotone clamping during global_evaluate remap.
Changes:
- Refactor mover boundary slip to use
mesh.boundary_slip(...)(including transientfacetrestore for unregistered labels) and add/realign tests around slip semantics. - Add fault-refinement metric builders (
fault_metric_tensor,fault_comb_metric,fault_metric,compose_metrics) plus new tests locking behavior and APIs. - Harden remesh field-transfer under MPI by clamping remapped values (
monotone="clamp") to prevent FE overshoot at newly-slipped boundary DOFs.
Reviewed changes
Copilot reviewed 18 out of 18 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
src/underworld3/meshing/smoothing.py |
Switch mover slip handling to mesh.boundary_slip, add slip_surfaces alias/deprecation plumbing, and introduce the mmpde mover implementation. |
src/underworld3/discretisation/discretisation_mesh.py |
Extend mesh.boundary_slip to build transient facet surfaces from reference facets and cache normals per build for performance. |
src/underworld3/meshing/bounding_surface.py |
Add reference_facets support and implement facet restore via nearest-point-on-facets. |
src/underworld3/meshing/_ot_adapt.py |
Add slip/facet helpers (_boundary_facets, _nearest_on_facets_*, _resolve_slip, _all_boundary_labels) used by mesh-owned slip and tests. |
src/underworld3/discretisation/remesh.py |
Enable monotone clamping during adapt field remap to prevent parallel boundary overshoot corruption. |
src/underworld3/function/_dminterp_wrapper.pyx |
Broaden the cell-hint DMLocatePoints bypass eligibility to simplex/manifold meshes (including serial) for performance/correctness consistency. |
src/underworld3/meshing/surfaces.py |
Fix Surface symbol collision and add fault metric builders + metric composition utilities. |
src/underworld3/meshing/__init__.py |
Export new fault metric and composition helpers. |
tests/test_0855_slip_surfaces.py |
New tests locking slip-surface semantics and mesh.boundary_slip API forms. |
tests/test_0763_boundary_slip_correctness.py |
New direct correctness tests for mesh-owned slip (radial/plane/facet/free behaviors). |
tests/test_0762_bounding_surfaces.py |
Update existing tests to reflect facet fallback behavior (unregistered labels slip rather than pin). |
tests/test_0762_fault_metric_tensor.py |
New test suite locking new fault metric builders and facade behavior. |
tests/test_0760_mesh_ot_adapt.py |
Update expectations: box boundaries now slide on faces (corners pinned) and relax error tolerance accordingly. |
tests/test_0750_meshing_follow_metric.py |
Add arc-length metric-choice tests and mark known follow-metric behavior differences as xfail. |
scripts/stagnant_lid_adapt_loop.py |
Update adaptation harness options (dt percentile, resolution ratio, GMG setup, collective safety fixes, richer logging). |
docs/developer/design/fault-refinement-simplification.md |
New design note documenting the consolidated fault refinement approach and tradeoffs. |
docs/developer/design/anisotropic-mmpde-mover.md |
New design document describing the MMPDE mover formulation and implementation details. |
docs/advanced/mesh-adaptation.md |
Update user docs to reflect method="mmpde" as the default mover and explain usage/knobs. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+3203
to
+3207
| if parallel: | ||
| vloc.array[:] = vel_loc.ravel() | ||
| coord_dm.localToGlobal(vloc, vglob, addv=True) | ||
| coord_dm.globalToLocal(vglob, vloc) | ||
| vel = np.asarray(vloc.array).reshape(-1, cdim).copy() |
Comment on lines
+2487
to
+2497
| from underworld3.meshing.smoothing import _edge_pairs | ||
| ep = _edge_pairs(mesh.dm) | ||
| Xc = np.asarray(mesh.X.coords) | ||
| if ep.shape[0]: | ||
| h0 = float(np.linalg.norm( | ||
| Xc[ep[:, 1]] - Xc[ep[:, 0]], axis=1).mean()) | ||
| else: | ||
| h0 = 1.0 | ||
| if uw.mpi.size > 1: | ||
| h0 = uw.mpi.comm.allreduce(h0) / uw.mpi.size | ||
| W = h0 / 6.0 |
Comment on lines
+2713
to
+2722
| from underworld3.meshing.smoothing import _edge_pairs | ||
| ep = _edge_pairs(mesh.dm) | ||
| Xc = np.asarray(mesh.X.coords) | ||
| if ep.shape[0]: | ||
| h0 = float(np.linalg.norm(Xc[ep[:, 1]] - Xc[ep[:, 0]], axis=1).mean()) | ||
| else: | ||
| h0 = 1.0 | ||
| if uw.mpi.size > 1: | ||
| h0 = uw.mpi.comm.allreduce(h0) / uw.mpi.size | ||
| return h0 |
Extend boundary-slip-strategy.md with the design we settled on for growing the tangent-slip contract from "the outer boundary" to "any surface the mesh must preserve under redistribution" (mesh-redistributor work, with codim-1 submesh extraction as the horizon): - per-surface declaration over the mesh's own topology — never an alternative topology; the mesh decides what its labels represent; - geometry is per-surface, never per-mesh (regional spherical: radial caps + plane great-circle sides — no mesh-level coordinate frame to inherit); - a submesh declares its OWN surfaces (borrow by reference, never re-home); - geometry-kind orthogonal to capabilities (tangent_moving / extractable); - tangent-movement is the broad preservation requirement — generalises the slip gate from is_bnd to "any tangent_moving surface" (internal interfaces need it or adaptation destroys them); the first concrete extension; - internal interfaces yes, faults no; - HDF5 persistence of analytic surface metadata for checkpoint roundtrip; - discoverability by inspection. Not implemented — agreed direction and the constraints it must honour. Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Anisotropic mesh-adaptation movers + parallel field-transfer correctness
This branch lands the validated anisotropic mesh-adaptation work and fixes a
parallel-only corruption in the adapt field-transfer. Three independent
pieces:
(a) Recover the PR #213 parallel-adaptation fix on this lineage
Merges
origin/developmentso the branch carries the #213parallel-singular-corruption fix (the robust barycentric in-cell locator /
DMLocatePointsbypass) and the serial hint-bypass perf fix(advect 58 s → 3 s). No behaviour change beyond what #213 already shipped.
(b) mmpde mover as default + arc-length metric
method="mmpde"anisotropic MMPDE mover, made the default; RBF-baked metricevaluation (~3.7× faster, same quality); simplex lock-out + dimension-general
cells (3D signed-volume fix).
follow_metric;fault-refinement layer (anisotropic + comb, composable).
slip_surfaces=) across the movers; scalardensity metric accepted as isotropic
ρ·I; non-finite line-search-step guard(a parallel deadlock fix) with opt-in stagnation exit;
on_boundarykwarg onthe in-cell test (accept on-face queries by default).
(c) Fix the parallel free-slip "leak": FE overshoot in the adapt field-transfer
On a parallel mover-adapted annulus the free-slip penalty BC appeared to leak
— nodal
v·nat the walls reached ~100 (vrms ~50) and advection then depositedhot T on the cold Dirichlet wall. This is not a BC bug and not DM-state
corruption. One-line fix: the adapt field-transfer (
remesh.py_remap_var_set) now evaluates withmonotone="clamp".What this bug looks like (read this before debugging the same symptom in free-surface work)
The free-surface / ALE path remaps fields through the same
global_evaluatetransfer machinery, so this exact failure can reappear there.The recognisable signature:
deformed mesh is always clean (post-solve nodal
v·n ≈ 1.5). The blow-upappears only after the adapt's field-transfer, never after a solve — so the
checkpoint (written post-adapt) reads ~100 while the run is dynamically fine
under cold-start.
clean: a boundary integral (
BdIntegralof(v·n)²) reads ~0.2, while thenodal
v.data(and the checkpoint) reads ~100. When an integral and anodal probe of the same field disagree by 100×, the corruption is in the
transferred nodal values, not in the field.
normals (
dm.computeCellGeometryFVM(face)) deviate ≤ 1e-2 from radial, labelcounts/checksums are stable, coordinate halos are bit-consistent
(
haloErr = 0), and the point SF is unchanged. (Selecting boundary faces bygetSupportSize == 1instead of by the DMLabel is a probe trap — it alsocatches parallel partition-seam interior faces and shows a phantom ~1.4
"normal deviation".)
parallel × adapt × boundary-sliponly. Any one factor off → clean.Mechanism
The transfer evaluates the OLD P2/P3 field at the NEW boundary DOF coordinates.
After tangential slip a new boundary node sits on the circle at a new angle, a
sagitta (~O(h²)) outside the OLD boundary cell (arc vs chord). In parallel
the swarm migrate lands that point in a containing cell on another rank,
where the FE basis overshoots — and the result is delivered as a "valid"
(un-flagged, not extrapolation) value. Same FE-overshoot family as the
SemiLagrangian trace-back bug.
Fix and validation
monotone="clamp"bounds each resampled value to its k-NN source-nodal range:bit-identical to plain FE in smooth regions, parallel-safe (rank-local).
max|v·n|118 → 3.0, hot-wall node count171 → 0.
bit-identical with/without the clamp (no serial regression); parallel
boundary error 1.36 → 0.045.
Also removes the abandoned GMG geometric-interpolation experiment
(
utilities/gmg_geometric_interpolation.py+ the no-op_pre_solve_hookseam);we settled on the default Galerkin-nested velocity multigrid. The work remains
in history if that route is revisited.
(d) Boundary tangent-slip: movers consume the mesh-owned
mesh.boundary_slip(step 2)The metric movers each carried their own boundary tangent-slip — duplicated and
geometry-coupled (guessing radial-ness at runtime, recovering the snap centre by
a per-call
allreduce). This factors it onto the mesh-ownedmesh.boundary_slipcontract that landed via #225 (
BoundingSurfaceobjects:radial/plane/facet/free, registered by the analytic constructors).mmpde/spring/ma/ot/anisotropic) now callmesh.boundary_slip(...); the private_build_slip_projector(+ the_gamma_p1_at_vertices/_label_vertex_mask/_boundary_centrehelpers)and the inline per-ring / box slip blocks are deleted (net −317 lines).
mesh.boundary_slipgains: afacetfallback (a transientfacetsurfacebuilt from the reference facets for labels with no analytic surface — instead
of pinning); normals computed once per build (avoids a
Gamma_P1re-solveinside the movers' line-search backtrack; a
DESIGN NOTEonBoundingSurface.normalsrecords the re-solve-vs-cached trade-off + escapehatch); and owned-only projection (parallel safety — ghosts via halo-sync).
spring/ma/anisotropicdrop their per-ring rotation anchor (mmpdenever had one; the signed-area guards prevent tangle).
otkeeps itsradial-only slip gating and now emits a
DeprecationWarning(incomplete;expected to be superseded by
mmpdewith a scalar metric).Validation: parity vs the deleted engine was machine precision (1.57e-16
annulus, 0.0 box); convection-harness A/B (serial 40-step + np=5) physics
unchanged. Tests: new
test_0763_boundary_slip_correctness;test_0855migratedonto
mesh.boundary_slip;test_0762covers the facet fallback + degenerategeometry; 3D radial registration moved to
test_0002_bounding_surface_3d.py(itmust run in the early CI batch —
SphericalShellconstruction is fragile underthe heavy batch's accumulated PETSc state, a pre-existing mesh-lifecycle issue
unrelated to slip).
The
origin/developmentmerge also brings this lineage current with #220(
estimate_dtpercentile/median cell reduction) and #225 (boundary-slipstep-1), plus the
solCx/ analytic-solution work.Underworld development team with AI support from Claude Code