Skip to content

Per-iteration SNES update callbacks (pressure gauge + boundary-correct generic scatter) - #250

Merged
lmoresi merged 3 commits into
developmentfrom
feature/snes-update-callbacks
Jun 18, 2026
Merged

Per-iteration SNES update callbacks (pressure gauge + boundary-correct generic scatter)#250
lmoresi merged 3 commits into
developmentfrom
feature/snes-update-callbacks

Conversation

@lmoresi

@lmoresi lmoresi commented Jun 18, 2026

Copy link
Copy Markdown
Member

Summary

Adds a user hook that runs at the start of every nonlinear (SNES) iteration via PETSc's SNESSetUpdate, for quantities that must track the current iterate rather than being recomputed once per timestep. Works on any solver.

solver.add_update_callback(fn)   # fn(solver, iteration)

Three parts:

  1. Callback mechanismadd_update_callback(fn) + the SNESSetUpdate dispatcher (scatter iterate→fields, run callbacks, gather fields→iterate, refresh aux vec; also applied once to the converged iterate). With no callback registered the solve path is byte-for-byte unchanged (the hook is only installed when a callback exists).

  2. Pressure-gauge helperSNES_Stokes.set_pressure_gauge(boundary, reference=0.0) removes the surface-mean pressure each iteration, pinning a specific physical gauge on enclosed (pressure-null-space) problems. Converges to machine-zero surface mean.

  3. Boundary-correct, solver-generic field scatter — the scatter that lets a callback read the iterate is now correct everywhere, including driven (non-zero Dirichlet) boundaries, on all solvers.

The scatter fix

A plain globalToLocal fills a field's interior and zero-Dirichlet DOFs but not its non-zero Dirichlet (driven) boundary DOFs — those aren't in the global vector; they're imposed on the local vector by DMPlexSNESComputeBoundaryFEM. Previously the dispatcher's scatter also assumed Stokes-only self.fields/_subdict, so a callback on a scalar/vector solver raised AttributeError.

  • Base class _scatter_global_to_fields / _gather_fields_to_global rewritten to be single-field and boundary-correct (globalToLocalDMPlexSNESComputeBoundaryFEMself.u), serving scalar / vector / multi-component solvers and their subclasses.
  • SNES_Stokes_SaddlePt overrides both to split its multi-field DM (velocity / pressure / block-constraint multipliers); the local index-set build is factored into _ensure_local_field_index_sets(), shared with the post-solve copy-back. Adds _multiplier_is cleanup to the rebuild reset.

This enables the gradient-plasticity / shear-band use case: an auxiliary Helmholtz/Projection smoother fired each iteration whose smoothed strain-rate field stays self-consistent with the velocity (measured error on a driven lid dropped from ~35% to 0).

Tests

tests/test_1016_snes_update_callbacks.py (tier_a / level_1, 7 passed):

  • callback fires; pressure gauge → machine-zero surface mean; no-callback path unchanged;
  • callback reads the driven-boundary velocity correctly (Stokes);
  • callback works on scalar (Poisson) and vector (Vector_Projection) solvers, reading non-zero Dirichlet boundaries correctly (both raised AttributeError before generalisation);
  • Helmholtz/Projection smoother self-consistent at convergence.

The boundary-scatter tests were confirmed to fail with the override disabled (genuine regression guards). test_1010, the Poisson and projection suites, and the constrained-Stokes tests still pass.

Docs

docs/advanced/solver-iteration-callbacks.md documents the generic, boundary-correct behaviour and both use cases.

Notes for review

  • New public API (add_update_callback, set_pressure_gauge) is introduced here on a feature branch; nothing else depends on it yet, so it's folded into this one PR rather than pre-extracted to development.
  • var.data (global-access path) still reads 0 on driven DOFs by design — the global vector never holds them; callbacks read through uw.function.evaluate/the aux-vec, which is now boundary-correct.

Underworld development team with AI support from Claude Code

lmoresi added 3 commits June 17, 2026 21:36
Adds a user-facing hook to run a callback at the start of every nonlinear
(Newton/SNES) iteration, for quantities that must track the current iterate
during the nonlinear solve rather than only per timestep.

API (SolverBaseClass): add_update_callback(fn) where fn(solver, iteration).
A dispatcher scatters the current iterate into the field variables, runs the
callbacks, gathers fields back into the iterate, and refreshes the mesh
auxiliary vector so the next residual sees any changes. Callbacks are also
applied to the final converged iterate (SNESSetUpdate fires only at the start
of each iteration). With no callback registered the solve path is unchanged
(verified bit-for-bit on test_1010; no-op guard).

Use cases landed:
  - set_pressure_gauge(boundary) on SNES_Stokes: removes the surface-mean
    pressure each iteration so the gauge is pinned (e.g. zero mean pressure on
    "Top"). Validated to machine precision on an enclosed lid-driven cavity.
  - firing a Helmholtz/Projection smoother each iteration (gradient-plasticity
    / shear-band stabilisation) — documented pattern.

Tests: tests/test_1016_snes_update_callbacks.py (callback fires; pressure gauge
zero-mean; no-callback path converges). Docs: docs/advanced/solver-iteration-callbacks.md.

KNOWN LIMITATION (tracked for follow-up): the dispatcher's field scatter uses
subdm.globalToLocal, which fills interior + zero-Dirichlet DOFs but NOT
non-zero Dirichlet (driven) boundary DOFs. A callback that reads the velocity
on a driven boundary therefore sees stale values there. The fix is to mirror
the post-solve copy-back (globalToLocal -> DMPlexSNESComputeBoundaryFEM ->
per-field local-IS extraction). The pressure gauge is unaffected (pressure has
no Dirichlet DOFs); the Helmholtz use case needs this scatter upgrade for
driven-boundary problems.

Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
Generalises the per-iteration SNES update-callback machinery so a callback can
read the solver's field(s) at the current Newton iterate on ANY solver, with the
values correct everywhere including on driven (non-zero Dirichlet) boundaries.

Before this change the dispatcher's scatter used Stokes-only `self.fields` /
`self._subdict` (so a callback on a scalar/vector solver raised AttributeError),
and even on Stokes it used a plain `subdm.globalToLocal` that does not carry
non-zero Dirichlet boundary DOFs — those are imposed on the LOCAL vector by
DMPlexSNESComputeBoundaryFEM. A smoother reading velocity on a driven lid saw a
stale value (~35% error in the Helmholtz self-consistency check).

- Base-class `_scatter_global_to_fields` / `_gather_fields_to_global` rewritten
  to be single-field and boundary-correct (globalToLocal -> boundary FEM ->
  `self.u`), serving scalar / vector / multi-component solvers.
- `SNES_Stokes_SaddlePt` overrides both to split its multi-field DM
  (velocity / pressure / block-constraint multipliers); the local index-set
  build is factored into `_ensure_local_field_index_sets()` shared with the
  post-solve copy-back. Adds `_multiplier_is` cleanup to the rebuild reset.
- No-callback path unchanged: the dispatcher is only installed when a callback
  is registered, so solves without callbacks are bit-identical.

Tests (tests/test_1016_snes_update_callbacks.py, tier_a/level_1):
- callback sees the driven-boundary velocity (Stokes);
- callback works on scalar (Poisson) and vector (Vector_Projection) solvers,
  reading non-zero Dirichlet boundaries correctly (both raised AttributeError
  before generalisation);
- Helmholtz/Projection smoother self-consistent at convergence (rel 0.35 -> 0).
All confirmed to fail with the scatter override disabled. test_1010, Poisson and
projection suites, and constrained-Stokes tests still pass.

Docs: docs/advanced/solver-iteration-callbacks.md updated for the generic,
boundary-correct behaviour; the resolved follow-up design brief removed.

Underworld development team with AI support from Claude Code
Copilot AI review requested due to automatic review settings June 18, 2026 06:50

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.

Pull request overview

Adds a solver-generic mechanism to run user-defined callbacks at the start of each PETSc SNES (Newton) iteration, enabling per-iterate updates such as pressure gauge fixing and auxiliary “inner” solves, and includes a boundary-correct iterate→field scatter so callbacks can reliably read driven Dirichlet boundary values.

Changes:

  • Introduces add_update_callback() + SNESSetUpdate dispatcher, plus boundary-correct scatter/gather utilities (with Stokes saddle-point specialisation).
  • Adds SNES_Stokes.set_pressure_gauge(boundary, reference=0.0) helper to enforce a surface-mean pressure gauge per iteration.
  • Adds tests and documentation covering callback firing, boundary correctness, scalar/vector solver support, and a Helmholtz/Projection smoother use case.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 3 comments.

Show a summary per file
File Description
tests/test_1016_snes_update_callbacks.py New regression tests for per-iteration callbacks, pressure gauge, and driven-boundary correctness across solvers.
src/underworld3/systems/solvers.py Adds set_pressure_gauge() helper on SNES_Stokes.
src/underworld3/cython/petsc_generic_snes_solvers.pyx Implements callback registry + SNES update hook, boundary-correct scatter/gather, and Stokes saddle-point overrides.
docs/advanced/solver-iteration-callbacks.md New documentation for the callback API, boundary-correct semantics, and example use cases.
docs/advanced/index.md Adds the new callbacks doc page to the advanced docs index.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Comment on lines +112 to +118
Callbacks run in registration order. Registering one forces a re-setup so
the PETSc ``SNESSetUpdate`` hook is attached. With no callbacks registered
no hook is installed and the solve path is byte-for-byte unchanged.
"""
self._snes_update_callbacks.append(callback)
self._needs_function_rewire = True
return callback
Comment on lines +862 to 867
# Attach the per-iteration callback dispatcher here (after all
# setFromOptions in the solve path). No-op when no callbacks registered.
self._maybe_install_snes_update()
self.snes.solve(None, gvec)
if divergence_retries <= 0:
return
Comment on lines +6984 to +6990
# SNESSetUpdate fires only at the START of each iteration, so the final
# converged iterate is otherwise un-hooked. Apply the callbacks once more
# to it (e.g. so a pressure gauge pins the FINAL pressure, and an
# auxiliary field is consistent with the converged solution).
if self._snes_update_callbacks:
self._dispatch_snes_update(self.snes, -1)

@lmoresi

lmoresi commented Jun 18, 2026

Copy link
Copy Markdown
Member Author

Review — approved ✅

Reviewed the core .pyx change (+234/−53) closely since it touches the solver:

  • The −53 lines are a clean extraction, not a behaviour change. The inline get_local_field_is block that built _pressure_is/_velocity_is/_multiplier_is is now the cached _ensure_local_field_index_sets helper with byte-identical logic, called from the post-solve copy-back (same site as before) and reused by the new mid-solve scatter. The existing Stokes copy-back path is preserved.
  • No-op when unused is guaranteed by construction: _snes_update_callbacks is empty by default → _maybe_install_snes_update attaches no SNESSetUpdate hook → solve path byte-for-byte unchanged. Verified the docstring claim against the code.
  • Boundary-correct scatter is right: _scatter_global_to_fields mirrors the post-solve copy-back (globalToLocalDMPlexSNESComputeBoundaryFEM → copy into the field), so a callback reading v/p on a driven (non-zero Dirichlet) boundary sees the imposed value, not a stale one. Stokes correctly overrides both scatter/gather for its multi-field DM.
  • set_pressure_gauge cleanly layers on add_update_callback (surface-mean removal each iteration) — the right altitude.
  • CI green incl. the 2 new tier_a tests in test_1016 (which fail with the override off). Docs added under docs/advanced/.

Merging for the release.

Underworld development team with AI support from Claude Code

@lmoresi
lmoresi merged commit 54b815c into development Jun 18, 2026
2 checks passed
@lmoresi

lmoresi commented Jun 18, 2026

Copy link
Copy Markdown
Member Author

Thanks for the review — all three points addressed in b82acea:

  1. Spurious rewire on add_update_callback — removed _needs_function_rewire = True. The SNESSetUpdate hook is already (re)attached at the start of every solve in _snes_solve_with_retries, so registering a callback no longer forces a JIT/pointwise-function rewire.

  2. Final-iterate dispatch only on Stokes — centralised the final _dispatch_snes_update(self.snes, -1) in _snes_solve_with_retries, so every solver (scalar/vector/multi-component as well as Stokes) applies callbacks once to the converged iterate, matching the docs. Added a test asserting the iteration == -1 dispatch fires on a scalar (Poisson) solver.

  3. Double-fire risk — removed the per-class dispatch from SNES_Stokes_SaddlePt.solve(); each solve calls the helper exactly once, so callbacks fire once on the final iterate.

Also renamed _maybe_install_snes_update_attach_snes_update_hook to match the existing _attach_constant_nullspace / _attach_stokes_nullspace convention. test_1016 (8) and test_1010 (6) pass.

gthyagi added a commit to gthyagi/underworld3 that referenced this pull request Jun 19, 2026
Revert the local Stokes_Constrained.topography default change. The solver API again defaults to raw multiplier topography and requires reference='mean' explicitly for enclosed gauge-fixed topography.

This matches issue underworldcode#254 and PR underworldcode#250: the gauge fix already exists through set_pressure_gauge(...) and explicit mean-stripped topography; benchmark scripts should enable those paths explicitly rather than changing solver defaults.

Validation:
- ./uw build
- python3 -m py_compile src/underworld3/systems/solvers.py tests/parallel/test_1063_constrained_freeslip_parallel.py
- pixi run -e amr-dev pytest -q tests/test_1061_constrained_freeslip.py::test_multiplier_and_topography_api
- pixi run -e amr-dev mpirun -np 2 pytest -q --with-mpi tests/parallel/test_1063_constrained_freeslip_parallel.py
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
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