Add smooth_mesh_interior — Winslow Jacobi smoother (PR A)#190
Merged
Conversation
Promotes a Winslow-style mesh smoother out of an experimental runner
(docs/developer/design/_phase_i_fs_convection_zoo.py) into the
underworld3.meshing module. Use case: after a deformation that has
crushed some cells (e.g. free-surface convection that has compressed
elements near the surface), apply Jacobi smoothing to interior
vertices to recover triangle quality. Topology, vertex IDs, DOF
mappings, and parallel partitioning are all preserved — only
coordinates change. The final mesh._deform_mesh call fires once at
the end of all sweeps so the rebuild/cache-invalidation cost is paid
once, not per sweep.
API:
from underworld3.meshing import smooth_mesh_interior
smooth_mesh_interior(
mesh,
pinned_labels=None, # None = auto-detect all named boundaries
n_iters=5,
alpha=0.5,
verbose=False,
)
The pinned_labels argument accepts a sequence of label names whose
vertices are held fixed. None (default) inspects mesh.boundaries and
pins every non-sentinel label — the usual case. An empty list pins
nothing (boundary vertices drift; mesh contracts).
Implementation: scipy CSR Mat-Vec on the row-normalised
vertex-vertex adjacency built from DMPlex edge iteration. Adjacency
is cached at module level keyed by (mesh-id, pinned-label-tuple,
topology) so repeated calls on the same mesh skip the rebuild.
Parallel: currently serial-only (raises NotImplementedError under
mpi.size > 1). A future change will swap the scipy Mat-Vec for a
PETSc Mat-Vec with halo exchange between sweeps.
Tests in tests/test_0850_mesh_smoothing.py cover:
- boundary vertices bit-identical before/after
- per-sweep displacement decreases monotonically (graph Laplacian
Jacobi is a contraction on the pinned-boundary problem)
- aspect-ratio quality does not worsen on a perturbed-interior mesh
- explicit pinned_labels list works
- empty pinned_labels list is legal (boundary drifts)
- pinned_labels=None auto-detects mesh.boundaries
Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR adds a new mesh-utility API, underworld3.meshing.smooth_mesh_interior, implementing a Winslow-style Jacobi smoother that relaxes interior vertex coordinates while keeping selected boundary-label vertices pinned. It promotes functionality into the supported meshing module and introduces regression tests for pinning behavior and smoothing outcomes.
Changes:
- Add
smooth_mesh_interior()implementation with adjacency construction and caching. - Export the new API from
underworld3.meshing. - Add a new pytest module covering boundary pinning and basic smoothing behavior.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 7 comments.
| File | Description |
|---|---|
src/underworld3/meshing/smoothing.py |
Implements the serial Winslow/Jacobi interior-vertex smoother with pinned-boundary handling and adjacency caching. |
src/underworld3/meshing/__init__.py |
Re-exports smooth_mesh_interior at module level via __all__. |
tests/test_0850_mesh_smoothing.py |
Adds regression tests for pinning semantics and smoothing behavior/quality checks. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+20
to
+28
|
|
||
| import numpy as np | ||
|
|
||
| import underworld3 as uw | ||
|
|
||
|
|
||
| # Cached adjacency keyed by (mesh-id, pinned-label-tuple, topology). | ||
| # Rebuilt automatically when the mesh topology changes. | ||
| _ADJ_CACHE: dict = {} |
Comment on lines
+185
to
+195
| dm = mesh.dm | ||
| pStart, pEnd = dm.getDepthStratum(0) | ||
| cStart, cEnd = dm.getHeightStratum(0) | ||
| cone_size = dm.getConeSize(cStart) if cEnd > cStart else 0 | ||
| cache_key = (id(mesh), pinned_labels, | ||
| pEnd - pStart, cEnd - cStart, cone_size) | ||
|
|
||
| cache = _ADJ_CACHE.get(cache_key) | ||
| if cache is None: | ||
| A, is_pinned = _build_adjacency(mesh, pinned_labels) | ||
| _ADJ_CACHE[cache_key] = (A, is_pinned) |
Comment on lines
+178
to
+179
| "Parallel (PETSc Mat-Vec) implementation pending — " | ||
| "see docs/developer/subsystems/mesh-smoothing.md.") |
Comment on lines
+101
to
+105
| mesh, | ||
| pinned_labels: Optional[Sequence[str]] = None, | ||
| n_iters: int = 5, | ||
| alpha: float = 0.5, | ||
| verbose: bool = False, |
Comment on lines
+212
to
+214
|
|
||
| # Single DM-coords update at the end: one rebuild, not N. | ||
| mesh._deform_mesh(coords) |
Comment on lines
+126
to
+132
| # Strictly decreasing — graph Laplacian Jacobi is a | ||
| # contraction on the (pinned-boundary) problem. | ||
| for k in range(1, len(disps)): | ||
| assert disps[k] < disps[k - 1], ( | ||
| f"Sweep {k+1} displacement {disps[k]:.3e} not " | ||
| f"less than sweep {k} displacement " | ||
| f"{disps[k-1]:.3e}; full series: {disps}") |
|
|
||
| if pinned_labels is None: | ||
| pinned_labels = _auto_pinned_labels(mesh) | ||
| pinned_labels = tuple(pinned_labels) |
Removes the serial-only restriction on smooth_mesh_interior. The
per-sweep update remains a local scipy CSR Mat-Vec (which already
includes ghost vertices in its rows/columns because DMPlex's local
chart contains the overlap layer), but between sweeps we now do a
PETSc halo exchange so each rank's ghost copies see the updated
owned values from neighbouring ranks.
Per sweep:
1. avg = A @ coords (local Mat-Vec; owned + ghost rows)
2. update owned interior (write to local numpy buffer)
3. localToGlobal INSERT (push owned values out to global)
4. globalToLocal (pull updated owned back into ghosts)
5. read refreshed local Vec into the numpy buffer for next sweep
Ghost-vertex identification via dm.getPointSF().getGraph() — the
SF leaves are the ghost points on this rank; everything else is
owned. The pinned-label set still pins both owned and ghost
boundary vertices (ghost boundary updates are dictated by their
owners anyway).
Parallel regression tests in
tests/parallel/test_0855_mesh_smoothing_parallel.py:
- boundary vertices bit-identical on every rank after smoothing
- ghost-vertex values agree exactly with their owners after the
smoother returns (verified by doing a fresh globalToLocal and
checking the local Vec is unchanged — invariance under halo
refresh = correct consistent state)
- global per-sweep interior displacement decreases monotonically
(matches the serial guarantee using MPI allreduce)
Tests pass under mpirun -np 2 and mpirun -np 4. Serial tests
(tests/test_0850_mesh_smoothing.py) still pass.
Underworld development team with AI support from Claude Code
Reason: UW3 distributes its DMs with cell-overlap-0. For an owned vertex on the rank partition cut, some of its incident edges live in cells owned by another rank that aren't in this rank's local stratum, so the per-rank scipy adjacency was missing entries and the smoother produced visibly wrong updates along every rank seam (np=2 lost 38 edges, np=4 lost 79 in a 16x16 box). Communicating coords between sweeps did not help because the matrix itself was incomplete. Replace the scipy CSR with a parallel PETSc AIJ Mat. Each rank inserts entries for every locally-visible edge using GLOBAL vertex indices, and mat.assemble() routes cross-rank contributions so that owned-vertex rows are complete. Per sweep: A.mult on each coord component into a global Vec, pointwise-divide by the precomputed degree vector, blend, write back to owned interior in the local numpy buffer. Results are bit-identical (one ULP) between serial and parallel runs at any rank count, regardless of where the partition cut lands. Add tests/parallel/test_0855 ::test_parallel_matches_serial_bit _identical which spawns a serial reference subprocess (with MPI env vars stripped so PETSc doesn't try to attach to the parent mpirun) and compares the gathered parallel result to the serial final coords, matching vertices by their pre-perturbation initial coordinate. Underworld development team with AI support from Claude Code
smooth_mesh_interior auto-pinning iterated every label in mesh.boundaries, but Annulus exposes a "Centre" pressure-pin label whose underlying DMLabel has an invalid communicator — any PETSc call on it (getNumValues, getValueIS, view) hard-aborts the interpreter (not a Python-catchable exception). Add "Centre" to the auto-pin skip list (it's a single-point pressure marker, not a geometric boundary) and keep the existing try/except guards in _pinned_mask for any future similar quirks. Annulus + smoother now runs end-to-end. AdvDiffusionSLCN.__init__: forward a new `theta` kwarg to both internally-constructed SemiLagrangian_DDt instances, mirroring the existing monotone_mode forwarding from #189. Lets callers configure the Adams-Moulton θ at construction time instead of patching adv_diff.DuDt.theta / DFDt.theta after the fact. Tested in the free-surface convection zoo: rk4 + theta=1.0 (BE) + monotone_mode="clamp" + smooth_mesh_interior every 2 steps runs 6 steps with T staying in [0,1] and no integrator instability. Underworld development team with AI support from Claude Code
Reason: UW3 mesh generators tag boundaries by EDGE. The vertex stratum of a boundary label sometimes misses 1-2 endpoint vertices at the gmsh seam (observed: θ=0°/180° on the Annulus outer rim — 100 of 102 rim vertices tagged, 2 seam vertices missed; similarly 50 of 52 on the inner rim). Pinning by vertex-stratum-only left those seam vertices free, so the smoother pulled them inward, producing visible dimples in the deformed surface that corrupted free-surface convection runs. Fix: in _pinned_mask, when a label tags an edge, pin both of the edge's endpoint vertices (closure of the tagged edges). Verified on res=16 Annulus: outer rim 102/102 and inner rim 52/52 pinned after the fix (was 100/102, 50/52). Underworld development team with AI support from Claude Code
lmoresi
added a commit
that referenced
this pull request
Jul 6, 2026
…og, value-first call-site sweep (WE-01..03,05,06,08,09,10) (#338) * docs(WE-01): adopt the one-governing-doc-per-topic authority map Repoint CLAUDE.md's Data Access 'Authoritative Reference' from the stale UW3_Style_and_Patterns_Guide.md to subsystems/data-access.md (the guide it crowned teaches patterns the code deprecates at runtime — DOC-04), and record the Style Charter §10 authority table in docs/developer/index.md as the master authority index. The Charter is added to the Getting Started toctree (removes a baseline 'not included in any toctree' warning). Finding: DOC-04 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-02): de-drift the Style Guide's four stale normative sections Rewrites the sections DOC-01 verified as contradicting the settled standards: - Docstring format: the 'Markdown Docstrings for pdoc/pdoc3' section is replaced by the NumPy/Sphinx RST standard (worked example with :math: and Parameters/Returns/Examples/Notes; conversion tracked in docs/plans/docstring-conversion-plan.md), per Style Charter section 6. - Doc file format: Quarto .qmd prescription (zero .qmd files exist in the repo) replaced by MyST .md/Sphinx guidance matching CLAUDE.md; migration table row updated. - Data access examples: 'Preferred' coordinate examples now use the real, runnable API — mesh.X.coords (read), mesh.deform() (coordinate changes), and the swarm.coords getter/setter for particle positions. The previous 'Preferred' example swarm.data += displacement raises AttributeError (getter-only property — SWARM-13 evidence); mesh.data warns at runtime. The private-attribute migration advice (swarm._particle_coordinates, mesh._deform_mesh presented as the NEW pattern) is deleted. - Front matter: the 21-line Quarto YAML header is replaced by a minimal MyST title block, and the guide now states that the UW3 Style Charter is the normative contract and wins on conflict. All replacement examples verified against current source: Swarm.coords setter (swarm.py), Mesh.deform (discretisation_mesh.py:3133), uw.synchronised_array_update / NDArray_With_Callback.delay_callbacks_global. Findings: DOC-01, SWARM-13 (style-guide part). Underworld development team with AI support from Claude Code * docs(WE-03): regenerate the docstring review queue; add the sweep to the release checklist The queue (last generated 2026-01-13, cdf5bb2) misrepresented the codebase both ways: it flagged now-complete items (solve, SNES_Scalar) as missing and contained zero entries for the June 2026 API (DOC-02). Regenerated over src/underworld3/**/*.py + **/*.pyx at the current tip. Two bugs in scripts/docstring_sweep.py's regex-based Cython parser made the regenerated queue lie about .pyx docstrings and are fixed as part of making the regeneration meaningful: - the indent group '(\s*)' with re.MULTILINE consumed preceding blank lines, shifting the computed definition line so the docstring search started ON the def/class line and always missed; - the docstring search started at the definition line rather than after the (possibly multi-line) signature, so long signatures hid their docstrings; - raw-string docstrings (r""", the norm in the solver .pyx) were not recognised. DOC-02 cross-validation on the regenerated queue now passes: solve / SNES_Scalar in the solver pyx are no longer flagged 'none'; the queue contains the June API (add_nitsche_bc, add_rotated_freeslip_bc, boundary_flux, set_custom_fmg, consistent_jacobian: 13 mentions) and flags the DOC-05 targets (Swarm.advection x2, read_timestep, write_proxy) as undocumented. Also adds the sweep to the quarterly release checklist (guides/release-process.md) so the queue cannot go stale unnoticed again. Findings: DOC-02 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-05): backfill the changelog for May - early July 2026; add the changelog sweep to the release checklist The changelog (the quarterly CIG/stakeholder record) ended in April 2026 while ~117 first-parent commits landed May through early July (DOC-03). Backfilled at the existing conceptual granularity — 14 grouped entries, grouped by subsystem rather than by PR, matching the established format (### Title (Month Year), bold lead sentence, hyphen bullets, inline PR references): - New '2026 Q3 (July - September)' section: the July 2026 quality campaign (#309-#313, #317, #322-#326, #329, #334 as grouped entries), rotated strong free-slip / boundary traction / dynamic topography (#293, #294, #298, #306), generalized geometric multigrid via custom prolongation (#290, #297), consistent Jacobian tangent (#258), swarm correctness (#216, #313, #323, #329), numpy 2 support (#301, #305). - Extended '2026 Q2' section with the May-June entries: mesh adaptation movers (#190, #209, #213, #228, #259, #264, #266), moving-mesh field transfer / deform() (#246, #249, #251), semi-Lagrangian accuracy controls (#164, #183, #185-#189, #208, #220), snapshot/checkpoint toolkit (#146, #195, #196, #198), Stokes_Constrained (#224, #229, #240, #265), local-h Nitsche + boundary-slip surfaces (#225, #241, #275), units interoperability (#277, #278, #283, #284), memory/evaluation/solver infrastructure (#161, #177-#179, #181, #182, #222, #237, #250, ...). Every entry is backed by a merged commit on development (verified against git log --first-parent aed517f..3184a40). Also adds a quarterly-changelog sweep step beside the docstring sweep in the release checklist (guides/release-process.md) per DOC-03's proposed fix. Findings: DOC-03 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-06): status headers on the unmarked design docs (per-doc git verification) Adds one-to-three-line Status markers to the 13 design docs that lacked one, following the directory's existing conventions (**Status**: line under the title; status: key inside existing YAML frontmatter for the three frontmatter-only docs), and corrects the stale 'Design Phase' marker on MATHEMATICAL_MIXIN_DESIGN.md (the mixin ships in utilities/mathematical_mixin.py). Every stamp was verified against git history (git log --follow dates) and the current source tree before writing: - Implemented: jacobian-consistent-tangent (PR #258, c63cd70), fmg-checkpoint-hierarchy (3cd73cd), petsc-dmplex-checkpoint-reload-plan (PR #146, write_timestep(petsc_reload=True) in tree), fault-refinement-simplification (smooth_mesh_interior / metric_density_from_gradient / fault_comb_metric all in tree), MATHEMATICAL_MIXIN_DESIGN. - Current reference/contract: mesh-adaptation-formulation, ND_UNITS_BOUNDARY_CONTRACT (PR #278, e0ece9a). - Investigation records (preserved via PR #245, 34a9dd4; production geometric-MG is custom prolongation, PR #290): snesfas-feasibility, snesfas-vanka-feasibility-study. - Design notes / prototypes with honest gaps: in_memory_checkpoint_design (not implemented, per its own trailing Status section), submesh-solver-architecture (extract_region/extract_surface exist; coarsened_companion does not). - Historical: ARCHITECTURE_ANALYSIS (persistence.py layout superseded), COORDINATE_MIGRATION_GUIDE (transition shipped), WHY_UNITS_NOT_DIMENSIONALITY (decision record). The audit's ~16 estimate over-counted: re-derived at this tip, 13 docs were unmarked plus one marked-but-stale (DOC-07). Findings: DOC-07 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-08): convert units.py public docstrings Google -> NumPy style Docstring-only conversion of the 18 public module-level functions that carried Google-style Args:/Returns:/Raises:/Examples: labels (check_units_consistency, get_dimensionality, get_units, non_dimensionalise, show_nondimensional_form, simplify_units, create_quantity, convert_units, to_base_units, to_reduced_units, to_compact, get_scaling_coefficients, set_scaling_coefficients, validate_expression_units, assert_dimensionality, validate_coordinates_dimensionality, enforce_units_consistency, require_units_if_active, convert_angle_to_degrees) to the NumPy/Sphinx standard (Style Charter section 6). dimensionalise was already NumPy style; one-line docstrings and private helpers are untouched. No code, signature, or behaviour changes (verified: every diff hunk is inside a docstring; ast.parse clean). Finding: API-12 (docs/reviews/2026-07/API-CONSISTENCY-REVIEW.md). Underworld development team with AI support from Claude Code * docs(WE-09): sweep call sites of the newer BC methods to value-first (conds, boundary, ...) order Wave C (#334) made the ORIGINAL value-first order canonical for add_nitsche_bc / add_rotated_freeslip_bc / add_constraint_bc (maintainer decisions D2/D3; Style Charter section 6) with deprecation shims for the legacy boundary-first and g= spellings. This sweep updates every call site of those THREE methods to the canonical order so nothing in the repository exercises the shims — 74 sites total: - tests/: 63 call sites across 12 files (test_1017, test_1018, test_1060, test_1061, test_1062, test_1064, test_1065 x2 serial; parallel test_1017, test_1062, test_1063, test_1064). tests/test_0641_wave_c_api_shims.py is deliberately untouched — its legacy-order calls ARE the deprecation contract. - docs/: 7 sites (curved-boundary-conditions.md x4, CONSTRAINED_FREESLIP_MULTIPLIER.md call + signature line, examples/submesh_investigation/test_region_ds_nitsche.py). - .claude/skills/: 3 sites (adapt-on-top-faults x2, free-surface-convection x1). - CLAUDE.md: 1 signature reference (free-slip BC preference section). The ~1,370 legacy-trio (add_dirichlet_bc/add_natural_bc/add_essential_bc) sites already conform and are untouched per the D2 decision. The audit review documents under docs/reviews/2026-07/ record the pre-decision state as evidence and are not swept. Discovered while verifying the swept tests run warning-free: the Wave C zero-datum guard in add_rotated_freeslip_bc rejects FLOAT zero (sympy.sympify(0.0) != 0 is structurally True), so the canonical add_rotated_freeslip_bc(0.0, boundary) raises NotImplementedError while conds=0 works. Filed as issue #336 with a TODO(BUG) marker at the guard (comment-only src touch); the swept call sites use the working integer form add_rotated_freeslip_bc(0, boundary). No fix applied here (Charter section 9 scope discipline). Findings: API-01/API-02 sweep (WE-09, REMEDIATION-WORKLIST.md). Underworld development team with AI support from Claude Code
Open
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.
Summary
Promotes a Winslow-style Jacobi mesh smoother from an experimental runner into the
underworld3.meshingmodule. Use case: after a deformation that has crushed some cells (e.g. free-surface convection compressing elements near the surface), apply Jacobi smoothing to interior vertices to recover triangle quality without changing topology.Topology, vertex IDs, DOF mappings, and parallel partitioning are preserved — only coordinates change. The final
mesh._deform_meshcall fires once at the end of all sweeps so the rebuild / cache-invalidation cost is paid once, not per sweep.API
pinned_labels=None(default): inspectmesh.boundariesand pin every non-sentinel label. The usual case.pinned_labels=['Top', 'Bottom']: pin only those boundaries; others drift.pinned_labels=[]: pin nothing — boundary vertices drift, mesh contracts.Implementation
dm.getPointSF().getGraph()— only owned interior vertices are written.(mesh-id, pinned-label-tuple, topology), rebuilt on topology change.Parallel
Per-sweep halo exchange via
coordDM.localToGlobal(INSERT) followed byglobalToLocal. Each rank computes local Mat-Vec on its own chart (including ghosts populated by the previous sync), updates only owned interior vertices, then pushes those updates out to receivers' ghost copies before the next sweep. Finalmesh._deform_meshinvalidates all the standard caches (already landed in #188).Follow-up PRs
mesh.adaptsemantics. Design notes in session memory; mid-loop swarm migration impractical so the spec accepts a per-sweep metric lag (small for typical N).Test plan
tests/test_0850_mesh_smoothing.py— 6 serial cases:pinned_labelslist workspinned_labelslist is legalNoneauto-detectsmesh.boundariestests/parallel/test_0855_mesh_smoothing_parallel.py— 3 parallel cases (pass on np=2 and np=4):Files changed
src/underworld3/meshing/smoothing.py(new, ~260 lines including parallel halo logic)src/underworld3/meshing/__init__.py(+5 lines: import + export)tests/test_0850_mesh_smoothing.py(new, ~180 lines)tests/parallel/test_0855_mesh_smoothing_parallel.py(new, ~120 lines)Underworld development team with AI support from Claude Code