From 9ee07f9f577f95605fe4cba6577627ff5ad12d31 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 22:58:50 +1000 Subject: [PATCH 01/17] feat(rotated-bc): prescribed-normal datum through the nonlinear SNES rotated path (#438) - v_n = u_n now carried through the manual Newton/Picard loop: accepted iterates are kept feasible (affine snap in _impose_normal_constraint), so increments stay homogeneous (n.delta = 0) and the datum never touches the tangent, the increment RHS, the custom-FMG prolongation or the nullspace handling. - A COLD start imposes the datum through the FIRST increment's affine lift (zeroRowsColumns x/b at the rest-state tangent, exactly the linear one-shot's treatment): snapping the zero state onto the datum manufactures a boundary-strip strain state whose shear-thinning tangent cannot descend (measured: first Newton step 4 orders too large, no line-searched alpha improves). Warm starts take the exact affine snap directly. - Collective datum-activity flag in BOTH solve paths: branching on the rank-local datum_map desyncs the ranks' collective sequences at np>1 (the two zeroRowsColumns variants scatter differently; the lift skips the line-search collectives). The linear path carried the same latent pattern, tolerated only because annulus partitions happen to give every rank a piece of the outer arc. - Regression: nonlinear power-law annulus with a u.n = cos(theta) datum imposed to machine precision through genuine Newton iteration (test_1018). Tests: test_1018 serial 18/18; test_1066 np2; test_1070 serial + np2 (7); test_1071 serial. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 121 ++++++++++++++++++------ tests/test_1018_rotated_freeslip.py | 44 +++++++++ 2 files changed, 134 insertions(+), 31 deletions(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index c0eb262f..c552e15e 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -16,6 +16,7 @@ """ import numpy as np import sympy +from mpi4py import MPI from petsc4py import PETSc # Rank-safe printing only (mpi.pprint). Safe at module top: underworld3.mpi is a @@ -371,6 +372,11 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo datum_specs = getattr(solver, "_rotated_freeslip_datum", None) Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) + # COLLECTIVE datum-activity flag: datum_map is RANK-LOCAL (only owners of + # datum boundary nodes hold entries), so branching on it directly desyncs the + # ranks' collective sequences (the two zeroRowsColumns variants below scatter + # differently) — the np>1 deadlock class. One reduction decides for everyone. + datum_active = dm.comm.tompi4py().allreduce(len(datum_map), MPI.SUM) > 0 # rotate:  = Q A Qᵀ, b̂ = Q b Ahat = Aorig.ptap(Qt) @@ -384,7 +390,7 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo # diagonal-based Schur approximations and MG smoothing in the boundary strip. Any # positive diagonal is exact — the solution rows are set explicitly. diag = _velocity_diag_scale(Ahat, solver) - if datum_map: + if datum_active: # v_n = ũ_n: pass the datum (rotated-frame normal component = v_n) as x so # zeroRowsColumns(rows, diag, x, b) sets b[rows]=diag·x[rows] AND subtracts the # eliminated columns' contribution A[:,cols]·x from the other rows. x=0 @@ -451,7 +457,9 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo # (the constrained matrix already drove the interior — only these rows were reset). # Restored BEFORE the residual report below, so the reported norm describes the # solution actually returned (the datum rows are consistent with b̂ by construction). - if datum_map: + # (_set_rows_local is purely local surgery — datum_active keeps the branch + # rank-consistent for readability, an empty local map is simply a no-op.) + if datum_active: _set_rows_local(Uhat, datum_map) # True rotated residual ‖·û − b̂‖ for the solve report (one matvec). Computed @@ -547,25 +555,32 @@ def _gather_fields_to_global(solver): return U -def _project_out_normal_component(u, Q, Qt, normal_rows): - """Impose the strong ``v_n = 0`` constraint exactly on the composite vector - ``u``, in place: rotate to the boundary frame (``û = Q u``), zero the - constrained normal rows, rotate back (``u = Qᵀ û``). Q is orthogonal, so this - is the exact projection onto the constraint-satisfying subspace.""" +def _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map=None): + """Impose the strong affine constraint ``v_n = ũ_n`` exactly on the composite + vector ``u``, in place: rotate to the boundary frame (``û = Q u``), zero the + constrained normal rows, set the prescribed-datum rows to ``sgn·ũ_n`` + (``datum_map``; empty/None ⇒ pure free-slip ``v_n = 0``), rotate back + (``u = Qᵀ û``). Q is orthogonal, so this is the exact affine projection onto + the constraint-satisfying set. Keeping every Newton iterate feasible this way + is what lets the increment problem stay HOMOGENEOUS (``n̂·δ = 0``) — the + datum never touches the tangent or the increment RHS.""" uh = u.duplicate() # transient projection buffer Q.mult(u, uh) _zero_rows_local(uh, normal_rows) + if datum_map: + _set_rows_local(uh, datum_map) Qt.mult(uh, u) uh.destroy() def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, - max_halvings=8): + datum_map=None, max_halvings=8): """Backtracking line search on ‖F̂‖ for the rotated Newton/Picard update ``u + α d`` (full step α=1 first, halved on failure). Cheap insurance far from the solution; α=1 is accepted immediately near it. Every trial iterate is - projected back onto the strong v_n=0 constraint before its residual is - measured. + snapped back onto the strong ``v_n = ũ_n`` constraint (``datum_map``; free-slip + when empty) before its residual is measured, so ‖F̂‖ is always evaluated at a + feasible point. Owns all its temporaries' destroys: on acceptance the input ``u`` is destroyed and replaced by the accepted iterate. Returns ``(u, improved)``; ``improved`` @@ -576,7 +591,7 @@ def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, for _ls in range(max_halvings): utry = u.copy() utry.axpy(alpha, d) - _project_out_normal_component(utry, Q, Qt, normal_rows) + _impose_normal_constraint(utry, Q, Qt, normal_rows, datum_map) Ftry = rotated_residual(utry) fnorm = Ftry.norm() Ftry.destroy() @@ -592,8 +607,10 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T verbose=False, zero_init_guess=True, picard=0, rtol=None, atol=1.0e-11, stol=1.0e-8, max_it=50): """Nonlinear rotated strong-free-slip solve: a manual outer Newton/Picard loop - that rotates the residual F(u), the Jacobian J(u) and the v_n=0 constraint EVERY - iteration, reusing the validated self-contained rotated fieldsplit-Schur solve + that rotates the residual F(u), the Jacobian J(u) and the strong ``v_n = ũ_n`` + constraint (``ũ_n = 0`` for pure free-slip; a prescribed wall-normal datum via + ``solver._rotated_freeslip_datum``) EVERY iteration, reusing the validated + self-contained rotated fieldsplit-Schur solve (``_solve_rotated_iterative``, incl. custom geometric FMG / GAMG velocity block and the rotated coupled null space) for each Newton increment. @@ -605,14 +622,20 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T strong constraint exactly at every iterate. Each iteration (unknown carried in the CARTESIAN frame ``u``; the increment is - solved in the rotated frame ``û = Q u``): + solved in the rotated frame ``û = Q u``). Every ACCEPTED iterate is FEASIBLE + (``v_n = ũ_n`` imposed exactly), so increments satisfy the homogeneous + constraint ``n̂·δ = 0`` and the datum never enters the tangent. The one + exception is a cold start with a non-zero datum, where the first increment + carries the datum jump through the affine lift (see the initial-guess comment + in the body — snapping the zero state onto the datum creates a boundary-strip + strain state whose nonlinear tangent is unusable): * ``F = computeFunction(u)`` → Cartesian residual (native essential BCs on other boundaries already applied by the DM); * ``F̂ = Q F``, zero ``F̂`` at the constrained normal rows (the constraint - residual is 0 there); converge on ‖F̂‖; + residual is 0 there — the iterate is feasible); converge on ‖F̂‖; * ``Ĵ = Q J(u) Qᵀ`` with ``zeroRowsColumns(normal_rows)``; * solve ``Ĵ δ̂ = −F̂``, ``δ = Qᵀ δ̂``, with a ‖F̂‖ backtracking line search; - * ``u += α δ`` and re-impose ``v_n = 0`` exactly. + * ``u += α δ`` and re-impose ``v_n = ũ_n`` exactly. The tangent used by ``computeJacobian`` is the solver's own (``consistent_jacobian`` → Picard / Newton / continuation), so the rotated loop inherits the same tangent @@ -683,26 +706,42 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T if continuation: solver._set_newton_alpha(0.0) # start in the Picard phase - # Q, the custom-FMG prolongation and the coupled null space depend only on the - # geometry / normals (NOT the solution), so build them ONCE and reuse each step. - if getattr(solver, "_rotated_freeslip_datum", None): - raise NotImplementedError( - "a prescribed non-zero wall-normal datum (u.n = ũ_n) on rotated free-slip is " - "implemented only for the LINEAR solve path (solve_rotated_freeslip); the " - "nonlinear path still imposes u.n = 0. Use a linear model for a prescribed " - "normal velocity, or extend this path.") - Q, Qt, normal_rows, _ = build_rotation(solver, boundaries) + # Q, the custom-FMG prolongation, the coupled null space and the datum values + # depend only on the geometry / normals / prescribed surface field (NOT the + # solution), so build them ONCE and reuse each step. The datum ũ_n enters the + # iteration ONLY through the feasible-iterate projection below: every iterate + # carries v_n = ũ_n exactly, so the Newton increment satisfies the HOMOGENEOUS + # constraint n̂·δ = 0 and the per-iteration operator treatment (zeroRowsColumns + # with no lift, tangent transparency, custom_Pl, nullspace) is unchanged from + # pure free-slip. + datum_specs = getattr(solver, "_rotated_freeslip_datum", None) + Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) nsp = _rotated_nullspace(solver, Q, normal_rows) - # initial guess (cartesian, composite): warm-start from the fields or zero, then - # impose v_n=0 exactly on it so the iteration starts feasible. + # initial guess (cartesian, composite): warm-start from the fields or zero. + # A prescribed datum on a COLD start is imposed through the FIRST increment's + # affine lift (zeroRowsColumns x/b — the same treatment as the linear one-shot), + # NOT by snapping the zero state: v_n = ũ_n over a zero interior manufactures an + # extreme boundary-strip strain rate whose shear-thinning tangent linearises so + # badly that no line-searched step descends (measured: first step 4 orders too + # large, f(α) > f(0) down to α=2⁻¹¹). Assembling the first tangent at rest and + # letting the increment carry the datum jump gives the smooth regularised- + # viscosity response instead; every later iterate is feasible and the + # increments are homogeneous. A warm start is assumed smooth (the previous + # converged state) and takes the exact affine snap directly. + # COLLECTIVE datum-activity flag (see the linear path): datum_map is rank-local, + # and the lift branch changes both the zeroRowsColumns variant and the + # line-search collectives — branching per-rank deadlocks at np>1. + datum_active = dm.comm.tompi4py().allreduce(len(datum_map), MPI.SUM) > 0 + lift_datum = datum_active and zero_init_guess if zero_init_guess: u = dm.createGlobalVec() u.set(0.0) else: u = _gather_fields_to_global(solver) - _project_out_normal_component(u, Q, Qt, normal_rows) + if not lift_datum: + _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map) J, Jp = snes.getJacobian()[:2] pres_is = solver._subdict["pressure"][0] @@ -767,9 +806,20 @@ def rotated_residual(uvec, keep_cartesian=False): Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent if diag_scale is None: diag_scale = _velocity_diag_scale(Ahat, solver) - Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) bhat = Fhat.copy() bhat.scale(-1.0) + if lift_datum and iters == 0: + # cold-start lift: the increment FROM ZERO carries the datum jump — + # b̂[rows] = diag·x̂ and the eliminated columns' A[:,cols]·x̂ are + # subtracted from the other rows (see the linear one-shot, which this + # reproduces operation-for-operation at the rest-state tangent). + xhat = bhat.duplicate() + xhat.zeroEntries() + _set_rows_local(xhat, datum_map) + Ahat.zeroRowsColumns(normal_rows, diag=diag_scale, x=xhat, b=bhat) + xhat.destroy() + else: + Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) dhat, last_reason, ctx = _solve_rotated_iterative( solver, Ahat, bhat, Q, Qt, normal_rows, custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) @@ -782,9 +832,18 @@ def rotated_residual(uvec, keep_cartesian=False): # level). ‖u‖=0 on a cold start ⇒ this never fires prematurely (d is large). step_converged = d.norm() <= stol * (u.norm() + 1e-30) improved = False - if not step_converged: + if lift_datum and iters == 0: + # accept the lift step unconditionally (it is the smooth rest-state + # response — the standard first-Picard warmup; a line search against + # ‖F̂(0)‖, which ignores the datum mismatch, would be meaningless) and + # snap the datum rows exactly (the KSP cleanup zeroed them). + u.axpy(1.0, d) + _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map) + step_converged = False + improved = True + elif not step_converged: u, improved = _backtracking_line_search( - u, d, rnorm, rotated_residual, Q, Qt, normal_rows) + u, d, rnorm, rotated_residual, Q, Qt, normal_rows, datum_map) # every per-iteration temporary dies HERE, on all exit paths dhat.destroy() d.destroy() diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index b74ee522..d988aafd 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -668,3 +668,47 @@ def test_rotated_freeslip_prescribed_normal_datum(): err = np.abs(vn - target).max() assert err < 1e-8, f"prescribed u.n=cos(theta) not imposed: max nodal error {err:.2e}" assert vn.max() > 0.9 and vn.min() < -0.9, "prescribed normal velocity magnitude wrong" + + +def test_rotated_freeslip_nonlinear_prescribed_normal_datum(): + """The prescribed wall-normal datum u.n = cos(theta) through the NONLINEAR rotated + path (power-law annulus): the loop genuinely iterates and every iterate is kept + feasible, so the converged solution carries the datum to the same (machine) + precision as the linear path — the constraint is exact independent of where the + nonlinear iteration stops.""" + RI, RO = 0.5, 1.0 + mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + nhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable("Vdn", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable("Pdn", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + g = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + e = 0.5 * (g + g.T) + eII = sympy.sqrt(0.5 * (e[0, 0] ** 2 + e[1, 1] ** 2) + e[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation gauge) + s.add_rotated_freeslip_bc(0, "Upper", normal=nhat) + s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero flux + s.consistent_jacobian = True # Newton tangent (few iterations) + s.petsc_use_pressure_nullspace = True + s.tolerance = 1e-7 + s.solve() + + info = s._rotated_freeslip_info + assert info["nonlinear_iterations"] > 1, "nonlinear datum solve did not genuinely iterate" + assert info["converged"], "nonlinear datum solve did not converge" + vc = v.coords + rr = np.hypot(vc[:, 0], vc[:, 1]) + outer = np.abs(rr - RO) < 4.0e-2 # outer velocity nodes (incl. edge mids) + rhat = vc[outer] / rr[outer, None] + vn = np.einsum("ij,ij->i", v.data[outer], rhat) + target = vc[outer, 0] / rr[outer] # cos(theta) at the same node coords + err = np.abs(vn - target).max() + assert err < 1e-8, f"nonlinear u.n=cos(theta) not imposed: max nodal error {err:.2e}" + assert vn.max() > 0.9 and vn.min() < -0.9, "prescribed normal velocity magnitude wrong" From b40c81bc77475f966d490da027687acf60be2dd6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 09:41:05 +1000 Subject: [PATCH 02/17] refactor(rotated-bc): datum-activity reduction via PETSc Vec norm, not bare MPI The rank-consistent datum-activity flag now rides the x-hat lift vector's collective norm (PETSc-before-MPI house rule): x-hat carries the prescribed rotated-frame datum for the zeroRowsColumns lift anyway, so one object does both jobs and the hand-rolled ownership loop collapses into _set_rows_local. Removes the mpi4py allreduce introduced in the previous commit. The deeper cure is a surface submesh (the 1D-manifold gap, #202 lineage): a boundary-field datum would make this flag a field norm and retire the hand-rolled surface bookkeeping in this layer - workstream C. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 52 ++++++++++++++----------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index c552e15e..2fc1eeae 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -16,7 +16,6 @@ """ import numpy as np import sympy -from mpi4py import MPI from petsc4py import PETSc # Rank-safe printing only (mpi.pprint). Safe at module top: underworld3.mpi is a @@ -372,17 +371,25 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo datum_specs = getattr(solver, "_rotated_freeslip_datum", None) Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) - # COLLECTIVE datum-activity flag: datum_map is RANK-LOCAL (only owners of - # datum boundary nodes hold entries), so branching on it directly desyncs the - # ranks' collective sequences (the two zeroRowsColumns variants below scatter - # differently) — the np>1 deadlock class. One reduction decides for everyone. - datum_active = dm.comm.tompi4py().allreduce(len(datum_map), MPI.SUM) > 0 # rotate:  = Q A Qᵀ, b̂ = Q b Ahat = Aorig.ptap(Qt) bhat = b.duplicate() Q.mult(b, bhat) + # x̂ carries the prescribed rotated-frame datum (sgn·ũ_n at the constrained + # rows) and does DOUBLE DUTY: its collective norm is the rank-consistent + # datum-activity flag. datum_map itself is RANK-LOCAL (only owners of datum + # boundary nodes hold entries), so branching on it per-rank desyncs the ranks' + # collective sequences (the two zeroRowsColumns variants below scatter + # differently) — the np>1 deadlock class. The reduction rides a PETSc Vec norm + # rather than a bare MPI call; a true surface submesh would make the datum a + # boundary field and this flag a field norm (the 1D-manifold gap, #202 lineage). + xhat = bhat.duplicate() + xhat.zeroEntries() + _set_rows_local(xhat, datum_map) + datum_active = xhat.norm() > 0.0 + # constrain rotated normal rows (v_n = datum, datum=0 for pure free-slip): zero the # matrix rows/cols and set the RHS at those rows to the datum. The constraint # diagonal is the mean |diag(A_vv)| rather than 1.0: unit diagonals amid O(eta/h^2) @@ -395,21 +402,13 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbo # zeroRowsColumns(rows, diag, x, b) sets b[rows]=diag·x[rows] AND subtracts the # eliminated columns' contribution A[:,cols]·x from the other rows. x=0 # reproduces the pure free-slip RHS exactly, so datum=0 is bit-identical. - xhat = bhat.duplicate() - xhat.zeroEntries() - rs, re = xhat.getOwnershipRange() - xa = xhat.getArray() - for g, val in datum_map.items(): - if rs <= g < re: - xa[g - rs] = val - xhat.setArray(xa) Ahat.zeroRowsColumns(normal_rows, diag=diag, x=xhat, b=bhat) - xhat.destroy() else: # zeroRowsColumns takes GLOBAL row indices; the RHS write goes through # _zero_rows_local (ownership-relative indexing — the np>1 crash class). Ahat.zeroRowsColumns(normal_rows, diag=diag) _zero_rows_local(bhat, normal_rows) + xhat.destroy() # ITERATIVE by default (LU is almost never right): a self-contained fieldsplit- # Schur solve whose velocity block is geometric FMG on the custom prolongation @@ -730,11 +729,20 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T # viscosity response instead; every later iterate is feasible and the # increments are homogeneous. A warm start is assumed smooth (the previous # converged state) and takes the exact affine snap directly. - # COLLECTIVE datum-activity flag (see the linear path): datum_map is rank-local, - # and the lift branch changes both the zeroRowsColumns variant and the - # line-search collectives — branching per-rank deadlocks at np>1. - datum_active = dm.comm.tompi4py().allreduce(len(datum_map), MPI.SUM) > 0 + # x̂ (the prescribed rotated-frame datum as a composite vector) does double + # duty, as in the linear path: its collective PETSc norm is the rank-consistent + # datum-activity flag (datum_map is rank-local — per-rank branching on it + # desyncs the zeroRowsColumns variant and the line-search collectives, the + # np>1 deadlock class), and it is the ready-made lift vector for the + # cold-start first increment. + xhat = dm.createGlobalVec() + xhat.set(0.0) + _set_rows_local(xhat, datum_map) + datum_active = xhat.norm() > 0.0 lift_datum = datum_active and zero_init_guess + if not lift_datum: + xhat.destroy() + xhat = None if zero_init_guess: u = dm.createGlobalVec() u.set(0.0) @@ -813,11 +821,9 @@ def rotated_residual(uvec, keep_cartesian=False): # b̂[rows] = diag·x̂ and the eliminated columns' A[:,cols]·x̂ are # subtracted from the other rows (see the linear one-shot, which this # reproduces operation-for-operation at the rest-state tangent). - xhat = bhat.duplicate() - xhat.zeroEntries() - _set_rows_local(xhat, datum_map) Ahat.zeroRowsColumns(normal_rows, diag=diag_scale, x=xhat, b=bhat) xhat.destroy() + xhat = None else: Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) dhat, last_reason, ctx = _solve_rotated_iterative( @@ -879,6 +885,8 @@ def rotated_residual(uvec, keep_cartesian=False): _destroy_rotated_ksp_ctx(ctx) # KSP/PC + the owned Schur pmat if Ahat is not None: Ahat.destroy() # the reused rotated operator + if xhat is not None: + xhat.destroy() # unused lift vector (loop exited before it) removed = _finalize_rotated_solution(solver, u, Q, normal_rows, remove_rotation_gauge) return {"Q": Q, "Qt": Qt, "reaction": reaction, "U": u, From b00fcace7786750824de5be87a4d267250e79ff2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 12:01:17 +1000 Subject: [PATCH 03/17] refactor(rotated-bc): one Newton loop for linear and nonlinear rotated solves; datum via conds (#438, #403) One constraint, one path: the linear one-shot is deleted and every rotated solve runs the manual Newton/Picard loop. A linear model converges after its first increment (the cold-start affine lift IS the linear solve), so the up-front _residual_is_nonlinear probe leaves the dispatch entirely - two Jacobian assemblies saved per linear solve (test_1018: 60s -> 55s). The probe survives only as a lazy guard in the picard + pure-Newton corner, where a linear model now ignores the meaningless warmup instead of raising. - solver._rotated_use_lu becomes an increment-solver choice (direct LU per Newton increment, serial PC-free diagnostic; naive pressure pin extracted to _naive_pressure_pin with its TODO(BUG)). Verified against the iterative path to 6e-10. - One result-dict shape: reaction = F(u) always (equals A.u - b exactly for linear residuals); the linear-only A/b keys and their consumer fallback in boundary_normal_traction retire. ksp_its is always the per-increment list (test_1018 3D assertion moved to max-per-increment). - add_rotated_freeslip_bc(conds, ...) accepts the non-zero wall-normal datum (number or scalar expression, value-first; #403 item 2) and stores it in _rotated_freeslip_datum; vector values rejected. Tests and FreeSurface now use the public argument instead of poking the attribute. - FreeSurface: the penalty-fallback NotImplementedError catch in _solve_consistent is gone - consistent_constraint="strong" now runs through the SNES rotated path for nonlinear rheologies. Linear solutions change at solver-tolerance level only (agreed relaxation of bit-preservation); an occasional cheap polish increment can follow the lift when the first residual check sits at the tolerance boundary (free-surface plume wall-time unchanged). Tests: test_1018 serial 18/18; test_1066 np2; test_1070 serial + np2 (7); test_1071 serial; LU-vs-iterative and lazy-probe smokes. Underworld development team with AI support from Claude Code --- CLAUDE.md | 2 +- .../cython/petsc_generic_snes_solvers.pyx | 106 +++--- src/underworld3/systems/free_surface.py | 30 +- src/underworld3/utilities/rotated_bc.py | 324 ++++++------------ .../test_1066_rotated_datum_parallel.py | 3 +- tests/test_1018_rotated_freeslip.py | 24 +- 6 files changed, 182 insertions(+), 307 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 0f4a1272..b5f471dd 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -326,7 +326,7 @@ The PETSc-based solvers are carefully optimized and validated. **NO CHANGES with ## Boundary Conditions: Free-slip -**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is the only implemented datum) +**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is free-slip; a non-zero scalar/expression `conds` prescribes the wall-normal datum `u·n̂ = ũ_n` strongly) to impose `v·n̂ = 0`: - Enforces zero wall-normal flow to **machine precision** (Nitsche / penalty leak ~1e-3). diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 0fda2ebd..1f148103 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1031,15 +1031,13 @@ class SolverBaseClass(uw_object): ``ksp.solve()`` on the rotated operator (``utilities/rotated_bc.py``) and never touches ``self.snes`` — so the generic reader would see stale SNES state. - Handles BOTH rotated result dicts: the linear one-shot solve (``ksp_its`` an - int, one outer solve, no outer verdict) and the manual nonlinear Newton loop - (``nonlinear_iterations``, ``ksp_its`` a per-increment LIST, and an outer - ``converged`` flag from the residual/step-norm tests). ``ksp_reason`` is a KSP - code on both paths (the nonlinear dict carries the LAST increment's), so it is - rendered with the KSP reason table — the SNES table shares the integers but - names different outcomes. A malformed dict raises: this reader and the dicts in - ``rotated_bc.py`` are a contract, and a silent fallback here previously masked - a broken one.""" + The rotated result dict is the manual Newton loop's (``nonlinear_iterations``, + ``ksp_its`` a per-increment LIST, an outer ``converged`` flag from the + residual/step-norm tests, and ``ksp_reason`` the LAST increment's KSP code) — + rendered with the KSP reason table, since the SNES table shares the integers + but names different outcomes. A malformed dict raises: this reader and the + dict in ``rotated_bc.py`` are a contract, and a silent fallback here + previously masked a broken one.""" from underworld3.systems.solve_report import SolveReport, ksp_reason_string info = info or {} @@ -5379,9 +5377,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._block_constraint_bcs = [] # Rotated strong free-slip BCs: [(boundary, normal), ...]. Registered via # add_rotated_freeslip_bc; when non-empty, solve() delegates to - # underworld3.utilities.rotated_bc (per-node DOF rotation + strong v_n=0 + - # reaction = sigma_nn). Empty by default → the solve path is unchanged. + # underworld3.utilities.rotated_bc (per-node DOF rotation + strong v_n = u_n + # + reaction = sigma_nn). Empty by default → the solve path is unchanged. + # _rotated_freeslip_datum maps boundary → prescribed wall-normal datum + # (the non-zero `conds` of add_rotated_freeslip_bc; absent ⇒ u.n = 0). self._rotated_freeslip_bcs = [] + self._rotated_freeslip_datum = {} self._rotated_freeslip_info = None # Give the Lagrange-multiplier (lambda) block its own viscosity-scaled # Schur preconditioner. The constraint Schur complement S_lambda = C A^-1 C^T @@ -5554,14 +5555,19 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Parameters ---------- - conds : float or None, optional - Prescribed wall-normal velocity datum, in the canonical value-first - BC order (Style Charter, API conventions). Only the homogeneous - free-slip constraint :math:`\mathbf{u}\cdot\hat{\mathbf n}=0` is - implemented, so this must be zero (or ``None``, meaning zero); a - non-zero datum raises ``NotImplementedError`` (use - :meth:`add_nitsche_bc` or ``add_constraint_bc`` for prescribed - normal in/outflow). + conds : float, scalar sympy expression, or None, optional + Prescribed wall-normal velocity datum :math:`\tilde u_n` (the SCALAR + component along the outward normal), in the canonical value-first BC + order (Style Charter, API conventions). Zero (or ``None``) is pure + free-slip :math:`\mathbf{u}\cdot\hat{\mathbf n}=0`. A non-zero value + — a number, a sympy expression of ``mesh.X``, or a scalar field read + such as ``h_dot.sym[0]`` — is imposed strongly (machine precision) as + :math:`\mathbf{u}\cdot\hat{\mathbf n}=\tilde u_n`, evaluated at the + boundary nodes at each ``solve()``. On an enclosed boundary the + datum must be discretely flux-free (:math:`\oint \tilde u_n = 0`) + for incompressibility. A corner/edge node shared between rotated + boundaries has no single normal and stays at the free-slip pinning + (datum ignored there). Vector/matrix values are rejected. boundary : str Boundary label to constrain. normal : None or sympy 1×dim Matrix or array, optional @@ -5605,14 +5611,17 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): raise TypeError( f"add_rotated_freeslip_bc() requires a boundary label string; " f"got {type(boundary).__name__}") - # Value comparison, not sympy's structural ==: Float(0.0) != Integer(0) - # structurally, but both are zero data. is_zero is True only when sympy - # can PROVE zero, so unprovable symbolic data is rejected too (#336). - if conds is not None and sympy.sympify(conds).is_zero is not True: - raise NotImplementedError( - "add_rotated_freeslip_bc imposes u.n = 0 only; a non-zero " - "wall-normal datum is not implemented (use add_nitsche_bc or " - "add_constraint_bc for prescribed normal in/outflow)") + if conds is not None: + if isinstance(conds, (sympy.MatrixBase, list, tuple, np.ndarray)): + raise TypeError( + "add_rotated_freeslip_bc(conds=...) takes the SCALAR " + "wall-normal datum u.n; got a vector/matrix value") + # Value comparison, not sympy's structural ==: Float(0.0) != Integer(0) + # structurally, but both are zero data (the free-slip member of the + # family). is_zero is True only when sympy can PROVE zero, so a + # symbolic datum (field read, expression) is correctly kept. + if sympy.sympify(conds).is_zero is not True: + self._rotated_freeslip_datum[boundary] = conds self._rotated_freeslip_bcs.append((boundary, normal)) self.is_setup = False return @@ -5632,11 +5641,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): viscosity ⇒ ``J`` independent of ``v`` ⇒ the two assemblies are bit-identical ⇒ linear. - Used to fail-fast on the rotated-free-slip path, which is a single linear - solve (assemble ``J(0)``, ``F(0)`` once) and would otherwise SILENTLY - return one Newton linearisation from ``u=0`` for a nonlinear model. The - caller must have run the pre-solve preamble (auxiliary vector + constants) - so the assembly sees the correct coefficients. + Used as a LAZY guard in the rotated-free-slip loop's picard + pure-Newton + corner (a Picard warmup is meaningless for a linear residual and impossible + for pure Newton) — the loop itself needs no up-front probe, it + self-terminates. The caller must have run the pre-solve preamble + (auxiliary vector + constants) so the assembly sees the correct + coefficients. """ snes = self.snes dm = self.dm @@ -8338,26 +8348,18 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.dm.setAuxiliaryVec(self.mesh.lvec, None) self._update_constants() - # Two paths, keyed on whether the model is nonlinear (detected by probing - # the assembled Jacobian for solution-dependence — a symbolic test cannot - # see the JIT-substituted strain-rate term): - # * LINEAR model → the validated one-shot solve (assemble J(0),F(0); - # rotate; single self-contained fieldsplit KSP). Fast, and the linear - # solution is exact, so warm-start / Picard warmup add nothing and are - # correctly ignored here. - # * NONLINEAR model → the manual outer Newton/Picard loop that rotates - # F(u), J(u) and the v_n=0 constraint every iteration. It honours - # zero_init_guess (warm start), the Picard warmup count, and the - # consistent_jacobian tangent (Picard / Newton / continuation). - from underworld3.utilities.rotated_bc import ( - solve_rotated_freeslip, solve_rotated_freeslip_nonlinear) - if not self._residual_is_nonlinear(): - self._rotated_freeslip_info = solve_rotated_freeslip( - self, self._rotated_freeslip_bcs, verbose=verbose) - else: - self._rotated_freeslip_info = solve_rotated_freeslip_nonlinear( - self, self._rotated_freeslip_bcs, verbose=verbose, - zero_init_guess=zero_init_guess, picard=picard) + # ONE path for linear and nonlinear models: the manual outer + # Newton/Picard loop that rotates F(u), J(u) and the strong + # v_n = u_n constraint every iteration. It honours zero_init_guess + # (warm start), the Picard warmup count, and the consistent_jacobian + # tangent (Picard / Newton / continuation); a linear model converges + # after its first increment, so no up-front nonlinearity probe is + # paid (the loop self-terminates; _residual_is_nonlinear survives as + # a lazy guard inside the picard+pure-Newton corner). + from underworld3.utilities.rotated_bc import solve_rotated_freeslip + self._rotated_freeslip_info = solve_rotated_freeslip( + self, self._rotated_freeslip_bcs, verbose=verbose, + zero_init_guess=zero_init_guess, picard=picard) # This path solves via ksp.solve on the rotated operator (not self.snes), # so give it a report from the rotated result rather than leaving a stale one. self._capture_rotated_report(self._rotated_freeslip_info) diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index ebe32542..50d34a9d 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -493,10 +493,11 @@ def _build_consistent(self): A rotated per-node constraint — the same primitive as the held lid, differing only in the constraint right-hand side, so free-slip is the ``\tilde u_n = 0`` member of the same family. It enforces the datum to machine precision (no - penalty leak) and routes the solve through the rotated LINEAR path (a direct - KSP solve) rather than a Newton line search, so an isoviscous flow does not - thrash ``newtonls``. It requires the datum to be discretely flux-free, which - :meth:`_surface_mean` now guarantees by weighting the demean by arc length. + penalty leak) through the unified rotated Newton loop: an isoviscous flow + converges in one increment (the cold-start affine lift IS the linear solve), + and a nonlinear rheology iterates with the datum held exactly at every + accepted iterate. It requires the datum to be discretely flux-free, which + :meth:`_surface_mean` guarantees by weighting the demean by arc length. ``"penalty"`` A weak natural BC. It tolerates a datum that carries a small net flux, at the @@ -515,9 +516,10 @@ def _build_consistent(self): if self.consistent_constraint == "strong": # The datum is read along the rotated per-node normal — the same deform # direction the rate was measured along, so no spurious tangential-slope term - # on a bumpy or rotating surface. - self.consistent.add_rotated_freeslip_bc(0.0, self.surface, normal=self.normal) - self.consistent._rotated_freeslip_datum = {self.surface: self._un_target.sym[0]} + # on a bumpy or rotating surface. Value-first: the ũ_n field read IS the + # conds datum, re-evaluated at the boundary nodes at each solve. + self.consistent.add_rotated_freeslip_bc( + self._un_target.sym[0], self.surface, normal=self.normal) else: n_hat = (self.normal if self.normal is not None else self.mesh.boundary_normal(self.surface)) @@ -958,19 +960,7 @@ def _solve_consistent(self, increment, dt): u_tilde = self._demean(increment / dt) self._un_target.array[...] = 0.0 self._un_target.array[self._un_target_rows, 0, 0] = u_tilde - try: - self.consistent.solve(zero_init_guess=True) - except NotImplementedError as exc: - # The rotated datum is implemented on the LINEAR rotated path only; the solve - # dispatches by a measured residual probe, so this fires exactly when the - # rheology is genuinely nonlinear. Name the knob rather than leaving the - # caller with the primitive's message. - raise NotImplementedError( - "FreeSurface: the strong rotated constraint cannot prescribe u.n = ũ_n " - "for a nonlinear rheology (the rotated datum is implemented on the linear " - "path only). Construct the manager with consistent_constraint='penalty' " - "to impose the material-surface rate weakly instead." - ) from exc + self.consistent.solve(zero_init_guess=True) def _conserve_composition(self): r"""Hold :math:`\int` (conserve integrand) fixed by a uniform shift of the diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 2fc1eeae..b012e738 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1,10 +1,14 @@ """Rotated strong free-slip for the Stokes saddle (the implementation behind ``solver.add_rotated_freeslip_bc``): build a per-node rotation Q from boundary -normals, rotate the assembled saddle Â=Q A Qᵀ / b̂=Q b, impose v_n=0 on the -rotated normal rows, solve, rotate back u=Qᵀû, remove the rigid-rotation gauge, -and expose σ_nn as the constraint reaction. - -The rotated saddle is solved by a self-contained fieldsplit-Schur KSP by default: the +normals, then drive ONE manual Newton/Picard loop that rotates the residual and +tangent every iteration (Ĵ=Q J Qᵀ, F̂=Q F), imposes v_n = ũ_n strongly on the +rotated normal rows (ũ_n = 0 is pure free-slip; a datum enters cold starts via the +first increment's affine lift, feasible iterates thereafter), rotates back u=Qᵀû, +removes the rigid-rotation gauge, and exposes σ_nn as the constraint reaction. +A linear model is simply the loop converging after its first increment — there is +no separate linear path and no up-front nonlinearity probe. + +Each increment is solved by a self-contained fieldsplit-Schur KSP by default: the velocity block is geometric FMG on the custom prolongation (``set_custom_fmg``) when a hierarchy is registered (the PREFERRED route), else GAMG tuned to the native path's settings; the Schur complement is preconditioned by the native 1/mu pressure mass @@ -318,176 +322,29 @@ def _set_rows_local(vec, row_val_map): # --------------------------------------------------------------------------- # # The rotated solve # --------------------------------------------------------------------------- # -def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbose=False): - """Assemble + solve the rotated strong-free-slip Stokes saddle. Fills the - solver's velocity/pressure fields with the (rotated-back, gauge-removed) - solution. - - Called from ``SNES_Stokes_SaddlePt.solve`` after ``_build`` (so the SNES/DM - exist); when used standalone it builds the solver itself. - - Returns - ------- - dict - The solve result consumed by ``boundary_normal_traction`` / - ``dynamic_topography_field`` (as their ``solve_result`` argument), with keys: - - * ``"Q"``, ``"Qt"`` — the rotation and its transpose (PETSc Mats); - * ``"A"``, ``"b"`` — the unrotated operator and RHS, kept so the reaction - ``A·u − b`` can be reconstructed (linear path only); - * ``"U"``, ``"Uhat"`` — the Cartesian and rotated-frame solutions; - * ``"normal_rows"`` — global rows of the constrained normal components; - * ``"boundaries"`` — the boundary specs the rotation was built from; - * ``"rotation_gauge_removed"`` — whether a rigid-rotation gauge was projected out; - * ``"ksp_reason"``, ``"ksp_its"`` — outer KSP converged-reason and iteration count; - * ``"rnorm"`` — true rotated residual ‖·û − b̂‖ (feeds the solve report). - """ - if getattr(solver, "snes", None) is None: - solver._setup_pointwise_functions() - solver._setup_discretisation() - solver._setup_solver() - dm = solver.dm - snes = solver.snes - - # Assemble the operator FIRST so its parallel row layout is final before we build - # Q against it (A = exact Jacobian at 0 — linear; b = -F(0)). The Pmat is - # assembled alongside: its p-p block is the native 1/mu pressure mass (the - # DS JacobianPreconditioner term) that preconditions the Schur complement — - # petsc4py's computeJacobian(x, J) would silently pass J as its own Pmat and - # the mass block would never be assembled. - snes.setUp() - U0 = dm.getGlobalVec() - U0.set(0.0) - J, Jp = snes.getJacobian()[:2] - snes.computeJacobian(U0, J, Jp) - Aorig = J.copy() - F0 = dm.getGlobalVec() - snes.computeFunction(U0, F0) - b = F0.copy() - b.scale(-1.0) - # borrowed temporaries → return to pool - dm.restoreGlobalVec(U0) - dm.restoreGlobalVec(F0) - - datum_specs = getattr(solver, "_rotated_freeslip_datum", None) - Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) - - # rotate:  = Q A Qᵀ, b̂ = Q b - Ahat = Aorig.ptap(Qt) - bhat = b.duplicate() - Q.mult(b, bhat) - - # x̂ carries the prescribed rotated-frame datum (sgn·ũ_n at the constrained - # rows) and does DOUBLE DUTY: its collective norm is the rank-consistent - # datum-activity flag. datum_map itself is RANK-LOCAL (only owners of datum - # boundary nodes hold entries), so branching on it per-rank desyncs the ranks' - # collective sequences (the two zeroRowsColumns variants below scatter - # differently) — the np>1 deadlock class. The reduction rides a PETSc Vec norm - # rather than a bare MPI call; a true surface submesh would make the datum a - # boundary field and this flag a field norm (the 1D-manifold gap, #202 lineage). - xhat = bhat.duplicate() - xhat.zeroEntries() - _set_rows_local(xhat, datum_map) - datum_active = xhat.norm() > 0.0 - - # constrain rotated normal rows (v_n = datum, datum=0 for pure free-slip): zero the - # matrix rows/cols and set the RHS at those rows to the datum. The constraint - # diagonal is the mean |diag(A_vv)| rather than 1.0: unit diagonals amid O(eta/h^2) - # viscous entries put a spectrum outlier on every rotated boundary node, poisoning - # diagonal-based Schur approximations and MG smoothing in the boundary strip. Any - # positive diagonal is exact — the solution rows are set explicitly. - diag = _velocity_diag_scale(Ahat, solver) - if datum_active: - # v_n = ũ_n: pass the datum (rotated-frame normal component = v_n) as x so - # zeroRowsColumns(rows, diag, x, b) sets b[rows]=diag·x[rows] AND subtracts the - # eliminated columns' contribution A[:,cols]·x from the other rows. x=0 - # reproduces the pure free-slip RHS exactly, so datum=0 is bit-identical. - Ahat.zeroRowsColumns(normal_rows, diag=diag, x=xhat, b=bhat) - else: - # zeroRowsColumns takes GLOBAL row indices; the RHS write goes through - # _zero_rows_local (ownership-relative indexing — the np>1 crash class). - Ahat.zeroRowsColumns(normal_rows, diag=diag) - _zero_rows_local(bhat, normal_rows) - xhat.destroy() - - # ITERATIVE by default (LU is almost never right): a self-contained fieldsplit- - # Schur solve whose velocity block is geometric FMG on the custom prolongation - # when a hierarchy is registered (set_custom_fmg), else GAMG. Direct LU only when - # explicitly opted in via solver._rotated_use_lu. - if getattr(solver, "_rotated_use_lu", False): - # Pin one pressure DOF (datum) — row only, keeps the B^T coupling. Only the - # direct solve needs it; the iterative path fixes the pressure gauge with the - # constant-pressure null space instead. - # TODO(BUG): the datum search is a naive per-rank scan of the rank's OWN - # global-section chart, so at np>1 each rank pins a different pressure DOF - # (or none, pin=None → zeroRows([None]) fails). Parallel-unsafe; tolerated - # only because LU is opt-in via solver._rotated_use_lu. - gsec = dm.getGlobalSection() - pin = None - pS, pE = gsec.getChart() - for q in range(pS, pE): - if gsec.getFieldDof(q, _PRESSURE_FIELD) > 0 \ - and gsec.getFieldOffset(q, _PRESSURE_FIELD) >= 0: - pin = gsec.getFieldOffset(q, _PRESSURE_FIELD) - break - Ahat.zeroRows([pin], diag=1.0) - if pin is not None: - _zero_rows_local(bhat, [pin]) - ksp = PETSc.KSP().create() - ksp.setOperators(Ahat) - ksp.setType("preonly") - pc = ksp.getPC() - pc.setType("lu") - pc.setFactorSolverType("mumps") - Uhat = dm.createGlobalVec() # returned in the result dict → own it - ksp.solve(bhat, Uhat) - ksp_reason = ksp.getConvergedReason() - ksp_its = ksp.getIterationNumber() - _warn_if_ksp_diverged(ksp, kind="rotated direct-LU") - else: - Mp = _pressure_mass_schur_pmat(solver) - Uhat, ksp_reason, ctx = _solve_rotated_iterative( - solver, Ahat, bhat, Q, Qt, normal_rows, verbose=verbose, Mp=Mp) - ksp_its = ctx["ksp"].getIterationNumber() - _destroy_rotated_ksp_ctx(ctx) - - # a prescribed normal datum: the iterative path zeros the rotated normal rows of - # the solution (its exact v_n=0 cleanup); restore v_n = ũ_n before rotating back - # (the constrained matrix already drove the interior — only these rows were reset). - # Restored BEFORE the residual report below, so the reported norm describes the - # solution actually returned (the datum rows are consistent with b̂ by construction). - # (_set_rows_local is purely local surgery — datum_active keeps the branch - # rank-consistent for readability, an empty local map is simply a no-op.) - if datum_active: - _set_rows_local(Uhat, datum_map) - - # True rotated residual ‖·û − b̂‖ for the solve report (one matvec). Computed - # explicitly rather than read off the KSP: preonly/LU never computes a norm, and - # the iterative norm can be the preconditioned one. - _res = bhat.duplicate() - Ahat.mult(Uhat, _res) - _res.axpy(-1.0, bhat) - rnorm = float(_res.norm()) - _res.destroy() - - # rotate back u = Qᵀ û (U is returned in the result dict → create, don't - # borrow from the pool) - U = dm.createGlobalVec() - Qt.mult(Uhat, U) - - removed = _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) - - return {"Q": Q, "Qt": Qt, "A": Aorig, "b": b, "U": U, "Uhat": Uhat, - "normal_rows": normal_rows, "boundaries": list(boundaries), - "rotation_gauge_removed": removed, "ksp_reason": ksp_reason, - "ksp_its": ksp_its, "rnorm": rnorm} +def _naive_pressure_pin(dm): + """One owned pressure DOF (datum) for the direct-LU gauge pin — row only, so + the B^T coupling is kept. Only the direct solve needs it; the iterative path + fixes the pressure gauge with the constant-pressure null space instead. + + TODO(BUG): a naive per-rank scan of the rank's OWN global-section chart, so at + np>1 each rank pins a different pressure DOF (or none → the caller must skip). + Parallel-unsafe; tolerated only because LU is opt-in via + ``solver._rotated_use_lu`` (a serial PC-free diagnostic).""" + gsec = dm.getGlobalSection() + pS, pE = gsec.getChart() + for q in range(pS, pE): + if gsec.getFieldDof(q, _PRESSURE_FIELD) > 0 \ + and gsec.getFieldOffset(q, _PRESSURE_FIELD) >= 0: + return gsec.getFieldOffset(q, _PRESSURE_FIELD) + return None def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge): """Remove the rigid-rotation gauge (if it is a genuine null space of the constrained problem), scatter the composite global vector ``U`` into the velocity/pressure fields, and refresh the enhanced-variable caches. Shared by - the linear one-shot and the nonlinear driver. Returns whether the gauge was + every exit of the rotated solve. Returns whether the gauge was removed.""" dm = solver.dm # Remove the rigid-rotation gauge — every mode that is a genuine null space of @@ -543,7 +400,7 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) def _gather_fields_to_global(solver): """Composite global vector built from the solver's current velocity/pressure - field values (the warm-start initial guess for the nonlinear driver).""" + field values (the warm-start initial guess for the rotated Newton loop).""" dm = solver.dm U = dm.createGlobalVec() U.set(0.0) @@ -602,17 +459,25 @@ def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, return u, False -def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=True, - verbose=False, zero_init_guess=True, picard=0, - rtol=None, atol=1.0e-11, stol=1.0e-8, max_it=50): - """Nonlinear rotated strong-free-slip solve: a manual outer Newton/Picard loop - that rotates the residual F(u), the Jacobian J(u) and the strong ``v_n = ũ_n`` - constraint (``ũ_n = 0`` for pure free-slip; a prescribed wall-normal datum via - ``solver._rotated_freeslip_datum``) EVERY iteration, reusing the validated - self-contained rotated fieldsplit-Schur solve +def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, + verbose=False, zero_init_guess=True, picard=0, + rtol=None, atol=1.0e-11, stol=1.0e-8, max_it=50): + """THE rotated strong-free-slip solve (linear and nonlinear models alike): a + manual outer Newton/Picard loop that rotates the residual F(u), the Jacobian + J(u) and the strong ``v_n = ũ_n`` constraint (``ũ_n = 0`` for pure free-slip; + a prescribed wall-normal datum via ``solver._rotated_freeslip_datum``) EVERY + iteration, reusing the validated self-contained rotated fieldsplit-Schur solve (``_solve_rotated_iterative``, incl. custom geometric FMG / GAMG velocity block and the rotated coupled null space) for each Newton increment. + A LINEAR model needs no separate path: the first increment from a cold start is + assembled at rest and already solves the problem (with a datum it is exactly the + affine-lift solve), so the loop converges at the next residual check — one + linear solve, plus two cheap residual evaluations. There is therefore no + up-front nonlinearity probe; the loop self-terminates. Direct LU per increment + (a serial PC-free diagnostic arbiter, e.g. for TI anisotropy experiments) is + opt-in via ``solver._rotated_use_lu``. + Why a manual loop rather than ``snes.solve()``: the rotated operator ``Q A Qᵀ`` (a ``ptap`` result) carries no DM field information, so PETSc's DM-coupled fieldsplit + geometric-MG cannot precondition it (SUBPC_ERROR); the increment @@ -665,7 +530,8 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T * ``"Q"``, ``"Qt"`` — the rotation and its transpose (PETSc Mats); * ``"reaction"`` — the converged Cartesian residual F(u), stashed as the - constraint reaction for σ_nn recovery (this path has no ``"A"``/``"b"``); + constraint reaction for σ_nn recovery (for a linear residual this equals + the assembled ``A·u − b`` exactly); * ``"U"`` — the Cartesian solution; * ``"normal_rows"`` — global rows of the constrained normal components; * ``"boundaries"`` — the boundary specs the rotation was built from; @@ -694,13 +560,20 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T mode = getattr(solver, "consistent_jacobian", False) continuation = (mode == "continuation") if picard and not continuation and mode is True: - raise NotImplementedError( - f"rotated free-slip: a Picard->Newton warmup (picard={picard}) with the " - "consistent Newton tangent (consistent_jacobian=True) requires " - "consistent_jacobian='continuation' — with pure Newton the frozen (Picard) " - "tangent needed for the warmup is not compiled in. Use " - "consistent_jacobian='continuation' (staged Picard then Newton), or drop " - "picard to run pure Newton.") + # The ONLY place linearity must be known in advance (hence the lazy probe — + # two trial assemblies, paid only by this contradictory configuration): a + # Picard warmup is meaningless for a linear residual, so ignore it there; + # for a genuinely nonlinear model pure Newton has no frozen tangent + # compiled in to warm up with, so fail loudly. + if solver._residual_is_nonlinear(): + raise NotImplementedError( + f"rotated free-slip: a Picard->Newton warmup (picard={picard}) with the " + "consistent Newton tangent (consistent_jacobian=True) requires " + "consistent_jacobian='continuation' — with pure Newton the frozen (Picard) " + "tangent needed for the warmup is not compiled in. Use " + "consistent_jacobian='continuation' (staged Picard then Newton), or drop " + "picard to run pure Newton.") + picard = 0 switch_rtol = max(float(getattr(solver, "newton_switch_rtol", 1.0e-2)), rtol) if continuation: solver._set_newton_alpha(0.0) # start in the Picard phase @@ -715,12 +588,15 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T # pure free-slip. datum_specs = getattr(solver, "_rotated_freeslip_datum", None) Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) - custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) - nsp = _rotated_nullspace(solver, Q, normal_rows) + # Direct LU per increment: no FMG prolongation / null space to prebuild — the + # gauge is fixed by the naive pressure pin instead (see _naive_pressure_pin). + use_lu = bool(getattr(solver, "_rotated_use_lu", False)) + custom_Pl = None if use_lu else _build_rotated_custom_Pl(solver, Q, normal_rows) + nsp = None if use_lu else _rotated_nullspace(solver, Q, normal_rows) # initial guess (cartesian, composite): warm-start from the fields or zero. # A prescribed datum on a COLD start is imposed through the FIRST increment's - # affine lift (zeroRowsColumns x/b — the same treatment as the linear one-shot), + # affine lift (zeroRowsColumns x/b at the rest-state tangent), # NOT by snapping the zero state: v_n = ũ_n over a zero interior manufactures an # extreme boundary-strip strain rate whose shear-thinning tangent linearises so # badly that no line-searched step descends (measured: first step 4 orders too @@ -808,10 +684,11 @@ def rotated_residual(uvec, keep_cartesian=False): Ahat = J.ptap(Qt) else: J.ptap(Qt, result=Ahat) # same nonzero pattern → in-place refresh - if ctx is None: - Mp = _pressure_mass_schur_pmat(solver) - elif Mp is not None: - Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent + if not use_lu: + if ctx is None: + Mp = _pressure_mass_schur_pmat(solver) + elif Mp is not None: + Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent if diag_scale is None: diag_scale = _velocity_diag_scale(Ahat, solver) bhat = Fhat.copy() @@ -819,16 +696,39 @@ def rotated_residual(uvec, keep_cartesian=False): if lift_datum and iters == 0: # cold-start lift: the increment FROM ZERO carries the datum jump — # b̂[rows] = diag·x̂ and the eliminated columns' A[:,cols]·x̂ are - # subtracted from the other rows (see the linear one-shot, which this - # reproduces operation-for-operation at the rest-state tangent). + # subtracted from the other rows. For a LINEAR model this increment IS + # the whole solve; the loop converges at the next residual check. Ahat.zeroRowsColumns(normal_rows, diag=diag_scale, x=xhat, b=bhat) xhat.destroy() xhat = None else: Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) - dhat, last_reason, ctx = _solve_rotated_iterative( - solver, Ahat, bhat, Q, Qt, normal_rows, - custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) + if use_lu: + if ctx is None: + ksp_lu = PETSc.KSP().create(comm=dm.comm) + ksp_lu.setType("preonly") + pc_lu = ksp_lu.getPC() + pc_lu.setType("lu") + pc_lu.setFactorSolverType("mumps") + ctx = {"ksp": ksp_lu, "Mp": None, + "pin": _naive_pressure_pin(dm)} + pin = ctx["pin"] + if pin is not None: + Ahat.zeroRows([pin], diag=1.0) + _zero_rows_local(bhat, [pin]) + ctx["ksp"].setOperators(Ahat) # values changed in place → refactor + dhat = Ahat.createVecRight() + dhat.set(0.0) + ctx["ksp"].solve(bhat, dhat) + last_reason = ctx["ksp"].getConvergedReason() + _warn_if_ksp_diverged(ctx["ksp"], kind="rotated direct-LU") + # parity with the iterative path's exact v_n cleanup (LU is exact only + # to factorisation round-off at the decoupled constraint rows) + _zero_rows_local(dhat, normal_rows) + else: + dhat, last_reason, ctx = _solve_rotated_iterative( + solver, Ahat, bhat, Q, Qt, normal_rows, + custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) lin_its.append(ctx["ksp"].getIterationNumber()) d = dm.createGlobalVec() Qt.mult(dhat, d) @@ -900,7 +800,7 @@ def rotated_residual(uvec, keep_cartesian=False): def _build_rotated_custom_Pl(solver, Q, normal_rows): """The rotated custom-FMG prolongation list [*coarse, Q_v·P_fine] for the velocity block, or None if no hierarchy is registered. Depends only on Q and - the mesh (NOT the solution), so the nonlinear driver builds it ONCE and reuses + the mesh (NOT the solution), so the rotated Newton loop builds it ONCE and reuses it across Newton iterations (the prolongation build is the expensive part).""" if getattr(solver, "_custom_mg", None) is None: return None @@ -983,7 +883,7 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal inherit it from the operator). ``custom_Pl`` / ``nsp`` may be PREBUILT (nonlinear driver: build once, reuse each - Newton step); when None they are built here (linear one-shot). + Newton step); when None they are built here. Returns ``(Uhat, reason, ctx)``. Passing ``ctx`` back in reuses the KSP/PC across Newton iterations — the fieldsplit ISs, Schur USER pmat and FMG @@ -1214,8 +1114,8 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): """Boundary normal traction σ_nn on `boundary` from the constraint reaction of the last rotated-free-slip solve (``solve_result`` is the dict returned by - ``solve_rotated_freeslip`` / ``solve_rotated_freeslip_nonlinear`` — see their - Returns sections for the keys). Returned mean-removed (the ρg·h gauge), as + ``solve_rotated_freeslip`` — see its + Returns section for the keys). Returned mean-removed (the ρg·h gauge), as ``(xs, sigma)`` with one entry per boundary velocity node on this rank. σ_nn is recovered from the CARTESIAN nodal reaction r_c = A·u − b: the nodal load @@ -1240,20 +1140,10 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): """ dm = solver.dm dim = solver.mesh.dim - # Cartesian nodal reaction r_c = F(u) at the converged state. The nonlinear - # driver stashes it directly (the final ``computeFunction`` residual); the linear - # one-shot reconstructs it as A·u−b (with A=J(0), b=−F(0), F affine ⇒ A·u−b=F(u)). - if solve_result.get("reaction") is not None: - rc = solve_result["reaction"] # owned by the result dict — do NOT destroy - own_rc = False - else: - A = solve_result["A"] - b = solve_result["b"] - U = solve_result["U"] - rc = A.createVecLeft() - A.mult(U, rc) - rc.axpy(-1.0, b) - own_rc = True + # Cartesian nodal reaction r_c = F(u) at the converged state, stashed by the + # solve driver (the final ``computeFunction`` residual; for a linear residual + # this equals the assembled A·u−b exactly). + rc = solve_result["reaction"] # owned by the result dict — do NOT destroy rcl = dm.getLocalVec() dm.globalToLocal(rc, rcl) rca = np.asarray(rcl.getArray()) @@ -1272,8 +1162,6 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): xs.append(_point_coord(dm, dim, cvec, csec, v0, v1, q)) Rn.append(float(np.dot(nrm, rcv))) # R_i = n̂·r_c (corner-correct) dm.restoreLocalVec(rcl) - if own_rc: - rc.destroy() # reconstructed A·u−b (linear path) — free it xs = np.array(xs) Rn = np.array(Rn) # σ_nn = −(nodal reaction), de-smeared by the SHARED boundary-mass primitive and diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index a11a5815..b86a6b79 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -41,8 +41,7 @@ def test_rotated_datum_prescribed_normal_partition_independent(): blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation) - s.add_rotated_freeslip_bc(0.0, "Upper", normal=nhat) - s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) # u.n = cos(theta), mean-zero s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" s.tolerance = 1.0e-9 diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index d988aafd..02d81ada 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -98,8 +98,8 @@ def test_rotated_freeslip_spherical_shell_3d(): info = s._rotated_freeslip_info assert info["ksp_reason"] > 0, f"rotated KSP diverged: {info['ksp_reason']}" - assert info["ksp_its"] <= 25, ( - f"Schur iteration blow-out: {info['ksp_its']} outer its " + assert max(info["ksp_its"]) <= 25, ( + f"Schur iteration blow-out: {info['ksp_its']} outer its per increment " "(1/mu-mass Schur preconditioning regressed?)") assert info["rotation_gauge_removed"], "3D rotation gauge not detected/removed" @@ -611,11 +611,11 @@ def test_rotated_freeslip_dynamic_topography_field(): # --- Prescribed non-zero wall-normal datum (u.n = ũ_n) --------------------------- -# The rotated constraint imposes u.n = datum strongly: datum=0 is pure free-slip (the -# held lid); a non-zero datum is the "consistent" material-surface velocity. Both share -# the same rotated matrix and differ only in the constraint RHS, so datum=0 must be -# BIT-IDENTICAL to plain free-slip. (Set via solver._rotated_freeslip_datum until the -# add_rotated_freeslip_bc datum argument lands — underworldcode/underworld3 tracking.) +# The rotated constraint imposes u.n = datum strongly, via the value-first conds +# argument: conds=0 is pure free-slip (the held lid); a non-zero conds is the +# "consistent" material-surface velocity. Both share the same rotated matrix and +# differ only in the constraint RHS, so an explicit zero datum must reproduce plain +# free-slip exactly. def _annulus_datum_solve(mode, tag): RI, RO = 0.5, 1.0 @@ -631,11 +631,8 @@ def _annulus_datum_solve(mode, tag): blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation gauge) - s.add_rotated_freeslip_bc(0, "Upper", normal=nhat) - if mode == "zero": - s._rotated_freeslip_datum = {"Upper": 0.0} - elif mode == "cos": - s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero => ∮u.n=0 + datum = {"plain": 0, "zero": 0.0, "cos": x / r}[mode] # cos: mean-zero => ∮u.n=0 + s.add_rotated_freeslip_bc(datum, "Upper", normal=nhat) s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" s.tolerance = 1e-9 @@ -693,8 +690,7 @@ def test_rotated_freeslip_nonlinear_prescribed_normal_datum(): blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation gauge) - s.add_rotated_freeslip_bc(0, "Upper", normal=nhat) - s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero flux + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) # u.n = cos(theta), mean-zero flux s.consistent_jacobian = True # Newton tangent (few iterations) s.petsc_use_pressure_nullspace = True s.tolerance = 1e-7 From 322d8fadf6826fec184df9b1dbe3c768e5ccbd70 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 12:34:24 +1000 Subject: [PATCH 04/17] fix(rotated-bc): honest step-norm exit; operator-verified rotation null modes Two defects surfaced by the TI fault smoke (#438 acceptance): 1. The step-norm exit reported converged=True on tiny increments while the residual was still O(1) (a stiff power-law tangent at the regularisation floor produces tiny steps far from the solution). The exit is now VERIFIED against the problem's intrinsic scale ||F(0)||: a warm start at the solution passes; a stagnating crawl goes back through the line search and either progresses (bounded by max_it) or ends on the stall exit honestly. 2. _mode_satisfies_constraints admitted rigid-rotation modes that are PINNED by an essential BC on another boundary (it only checked the rotated constraint rows). The admitted mode was projected out of every increment RHS, leaving an irreducible residual component: Newton converged superlinearly to rel ~2e-5 and floored there (measured on the essential-inner + rotated-outer power-law annulus). Modes are now also verified as null vectors of the ASSEMBLED operator (||J.m|| against the velocity diagonal scale), which catches any pinning without enumerating BC types; the null space is built after the first Jacobian assembly to make that possible. Same fix protects the post-solve rotation-gauge removal. With the fix the same solve converges to rel 2e-13 in 11 iterations and correctly leaves the gauge in place. TI fault-bearing smoke (power-law background, weak-plane band, cos(theta) datum) now passes: 24 its to rel 6e-8, datum at the quadrature floor, sigma_nn recovery clean. Tests: test_1018 serial 18/18 (incl. gauge-removal cases), test_1066 np2 (reference unchanged - the linear-case imprint of the spurious projection is below the 1e-6 test tolerance), test_1070 serial (7) + test_1071. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/rotated_bc.py | 80 ++++++++++++++++++++----- 1 file changed, 65 insertions(+), 15 deletions(-) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index b012e738..428cf2b0 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -588,11 +588,14 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, # pure free-slip. datum_specs = getattr(solver, "_rotated_freeslip_datum", None) Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) - # Direct LU per increment: no FMG prolongation / null space to prebuild — the + # Direct LU per increment: no FMG prolongation / null space to build — the # gauge is fixed by the naive pressure pin instead (see _naive_pressure_pin). use_lu = bool(getattr(solver, "_rotated_use_lu", False)) custom_Pl = None if use_lu else _build_rotated_custom_Pl(solver, Q, normal_rows) - nsp = None if use_lu else _rotated_nullspace(solver, Q, normal_rows) + # The null space is built INSIDE the loop, after the first Jacobian assembly: + # _mode_satisfies_constraints verifies each candidate mode against the + # ASSEMBLED operator (‖J·m‖ ≈ 0), which an unassembled J cannot support. + nsp = None # initial guess (cartesian, composite): warm-start from the fields or zero. # A prescribed datum on a COLD start is imposed through the FIRST increment's @@ -687,6 +690,7 @@ def rotated_residual(uvec, keep_cartesian=False): if not use_lu: if ctx is None: Mp = _pressure_mass_schur_pmat(solver) + nsp = _rotated_nullspace(solver, Q, normal_rows) # J assembled above elif Mp is not None: Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent if diag_scale is None: @@ -736,7 +740,28 @@ def rotated_residual(uvec, keep_cartesian=False): # are at the solution — the exit for a warm start that is already converged # (otherwise the relative test above, with a tiny r0, chatters near machine # level). ‖u‖=0 on a cold start ⇒ this never fires prematurely (d is large). + # A tiny step alone does NOT prove convergence — a stiff tangent (power-law + # at the regularisation floor) also produces tiny increments at a large + # residual — so the exit is VERIFIED against the problem's intrinsic scale + # ‖F̂(0)‖: a warm start at the solution passes (its residual sits at the + # cold chain's converged level); a stagnating crawl does not. One extra + # residual evaluation, paid only on this rare path. On verification + # failure the tiny step goes through the NORMAL line search instead: if it + # still improves, the crawl proceeds (bounded by max_it); if not, the + # stall exit ends the loop honestly. step_converged = d.norm() <= stol * (u.norm() + 1e-30) + if step_converged: + z = dm.getGlobalVec() + z.set(0.0) + F0hat = rotated_residual(z) + fnorm0 = F0hat.norm() + F0hat.destroy() + dm.restoreGlobalVec(z) + step_converged = rnorm <= rtol * fnorm0 + atol + if not step_converged: + mpi.pprint(f"[rotated_bc] step-norm at rel |F̂| = " + f"{rnorm / (fnorm0 + 1e-300):.2e} of the rest-state " + f"residual — stagnation-or-crawl, continuing to iterate.") improved = False if lift_datum and iters == 0: # accept the lift step unconditionally (it is the smooth rest-state @@ -882,8 +907,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal complement for enclosed domains (the hand-built IS fieldsplit does not inherit it from the operator). - ``custom_Pl`` / ``nsp`` may be PREBUILT (nonlinear driver: build once, reuse each - Newton step); when None they are built here. + ``custom_Pl`` / ``nsp`` are built by the CALLER (once, reused each Newton step; + the null-space mode verification needs the assembled Jacobian, which only the + caller can guarantee); None simply means "no hierarchy" / "no null modes". Returns ``(Uhat, reason, ctx)``. Passing ``ctx`` back in reuses the KSP/PC across Newton iterations — the fieldsplit ISs, Schur USER pmat and FMG @@ -896,12 +922,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal pres_is = solver._subdict["pressure"][0] if ctx is None: - if custom_Pl is None: - custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) - - # rotated coupled null space (pressure-const ⊕ Q·rotation) on the operator - if nsp is None: - nsp = _rotated_nullspace(solver, Q, normal_rows) + # rotated coupled null space (pressure-const ⊕ Q·rotation) on the operator — + # built by the CALLER (the Newton loop, after the first Jacobian assembly: + # the mode verification needs the assembled operator). if nsp is not None: Ahat.setNullSpace(nsp) Ahat.setTransposeNullSpace(nsp) @@ -1083,10 +1106,23 @@ def _rotated_nullspace(solver, Q, normal_rows): def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): - """True iff the rigid-body mode ``tg`` satisfies all rotated v_n=0 - constraints — i.e. Q·tg is ~0 on every constrained normal row. (A closed - circular boundary admits its one rotation; a full spherical shell admits all - three; straight/partial walls pin them.) + """True iff the rigid-body mode ``tg`` is a genuine null mode of the + constrained problem: it satisfies all rotated v_n=0 constraints (Q·tg ~0 on + every constrained normal row) AND it is a null vector of the assembled + operator (‖J·tg‖ ~0 against the velocity diagonal scale). The operator test + is what catches pinning the rotated rows cannot see — an ESSENTIAL no-slip + on another boundary leaves the mode tangential to the rotated wall (first + test passes) while the eliminated operator is NOT null on it; admitting it + then projects an irreducible component out of every increment RHS, and the + Newton residual floors at that component's magnitude instead of converging + (measured: rel ~2e-5 plateau on the essential-inner + rotated-outer annulus). + (A closed circular free-slip boundary admits its one rotation; a full + spherical shell admits all three; straight/partial walls and any essential + BC pin them.) + + The caller must have ASSEMBLED the solver's Jacobian at some iterate; an + unassembled J (norm 0) skips the operator test rather than passing every + mode through a vacuous 0 ≈ 0. COLLECTIVE: every rank runs the same global-vector ops. Do NOT early-return on a per-rank ``not normal_rows`` — in parallel a rank may own no boundary node @@ -1108,7 +1144,21 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): # transient duplicates tr.destroy() trc.destroy() - return viol < tol + if viol >= tol: + return False + # operator-nullity: rigid rotation has exactly zero strain in the discrete + # space (P2 contains linear fields, affine cells integrate the form exactly), + # so a genuine null mode gives assembly round-off; a pinned mode leaves O(1) + # boundary-strip rows. + J = solver.snes.getJacobian()[0] + Jm = tg.duplicate() + J.mult(tg, Jm) + jn = Jm.norm() + Jm.destroy() + if jn == 0.0 and J.norm() == 0.0: # J never assembled → cannot verify + return True + op_viol = jn / (_velocity_diag_scale(J, solver) * (tg.norm() + 1e-30)) + return op_viol < tol def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): From d8df7a1951cc950c77b37f675c04b53f9bf124a8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 13:01:15 +1000 Subject: [PATCH 05/17] feat(free-surface): warm-start the consistent datum solve; reference-scale convergence for warm starts - FreeSurface._solve_consistent now warm-starts from the free solve's fields: the consistent solution is the free solution with a small material-boundary datum correction, and the free solve has already converged this step. Starting there keeps a power-law tangent at physical strain rates - a cold start puts it at the regularisation floor, where the Newton line search stalls at O(0.1) relative residual (measured, acceptance attempt 2). With the warm start the first acceptance step converges in 3 iterations. - The rotated loop's convergence reference is now max(r0, ||F(0)||): rtol relative to a warm start's own (small) initial residual demands ever-more absolute accuracy the better the guess (the warm-start rtol trap). The rest-state residual is the intrinsic forcing scale; cold starts unchanged (their r0 IS that scale). The step-norm verification and the non-convergence warning read the same reference. Acceptance (#438) evidence in ~/+Simulations/FreeSurface/annulus_fs_convection/b_snes_datum_acceptance/: power-law annulus FS convection with consistent_constraint="strong" through the SNES rotated path holds the material boundary at the strong-datum level (err 5.3e-3 vs penalty 3.7e-2; net flux 1.4e-4 vs 2.0e-3) at matched steps; TI fault-bearing smoke passes (24 its to rel 6e-8). Open finding recorded there: deformed-ring datum compatibility floors the interior residual at rel ~2e-3 (honestly reported unconverged; suspected arc-length-vs-FE-quadrature demean gap - the principled fix is a trace-mass-weighted demean). Tests: test_1018 18/18 serial; test_1070 serial+np2 (7); test_1071; test_1066 np2. Underworld development team with AI support from Claude Code --- src/underworld3/systems/free_surface.py | 10 +++++- src/underworld3/utilities/rotated_bc.py | 47 +++++++++++++++---------- 2 files changed, 38 insertions(+), 19 deletions(-) diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index 50d34a9d..1fbbbb4d 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -960,7 +960,15 @@ def _solve_consistent(self, increment, dt): u_tilde = self._demean(increment / dt) self._un_target.array[...] = 0.0 self._un_target.array[self._un_target_rows, 0, 0] = u_tilde - self.consistent.solve(zero_init_guess=True) + # Warm-start from the free solve: the consistent solution IS the free + # solution with the (small) material-boundary datum imposed, and the free + # solve has already converged this step. Starting there keeps a power-law + # tangent at physical strain rates — a cold start puts it at the + # regularisation floor, where the Newton line search stalls at O(0.1) + # relative residual (measured, power-law annulus acceptance run). + self.consistent.u.array[...] = self.free.u.array + self.consistent.p.array[...] = self.free.p.array + self.consistent.solve(zero_init_guess=False) def _conserve_composition(self): r"""Hold :math:`\int` (conserve integrand) fixed by a uniform shift of the diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 428cf2b0..448a7297 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -654,7 +654,24 @@ def rotated_residual(uvec, keep_cartesian=False): _zero_rows_local(Fh, normal_rows) return Fh + # Convergence reference: max(initial residual, REST-STATE residual ‖F̂(0)‖). + # A warm start's own initial residual is small, and rtol relative to it + # demands ever-more absolute accuracy the better the guess (the warm-start + # rtol trap: the FS time loop stalled at rel~2e-3 of its warm start while + # being far below tolerance on the physical scale). The rest-state residual + # is the problem's intrinsic forcing scale; a cold start's r0 IS that scale, + # so nothing changes there. + fnorm_rest = None + if not zero_init_guess: + z = dm.getGlobalVec() + z.set(0.0) + Fz = rotated_residual(z) + fnorm_rest = Fz.norm() + Fz.destroy() + dm.restoreGlobalVec(z) + r0 = None + ref = None last_reason = 0 iters = 0 converged = False @@ -664,12 +681,13 @@ def rotated_residual(uvec, keep_cartesian=False): rnorm = Fhat.norm() if r0 is None: r0 = rnorm + ref = max(r0, fnorm_rest) if fnorm_rest is not None else r0 if verbose: mpi.pprint(f"[rotated_bc] nonlinear iter {iters:2d} |F̂|={rnorm:.6e} " - f"rel={rnorm/(r0+1e-300):.3e} [{phase}]") - # residual convergence (relative to the initial residual, plus an absolute + f"rel={rnorm/(ref+1e-300):.3e} [{phase}]") + # residual convergence (relative to the reference scale, plus an absolute # floor so an already-converged warm start does not chase machine noise). - if rnorm <= rtol * r0 + atol: + if rnorm <= rtol * ref + atol: converged = True Fhat.destroy() break @@ -742,25 +760,18 @@ def rotated_residual(uvec, keep_cartesian=False): # level). ‖u‖=0 on a cold start ⇒ this never fires prematurely (d is large). # A tiny step alone does NOT prove convergence — a stiff tangent (power-law # at the regularisation floor) also produces tiny increments at a large - # residual — so the exit is VERIFIED against the problem's intrinsic scale - # ‖F̂(0)‖: a warm start at the solution passes (its residual sits at the - # cold chain's converged level); a stagnating crawl does not. One extra - # residual evaluation, paid only on this rare path. On verification + # residual — so the exit is VERIFIED against the reference scale (which + # includes the rest-state residual ‖F̂(0)‖ for warm starts): a warm start + # at the solution passes; a stagnating crawl does not. On verification # failure the tiny step goes through the NORMAL line search instead: if it # still improves, the crawl proceeds (bounded by max_it); if not, the # stall exit ends the loop honestly. step_converged = d.norm() <= stol * (u.norm() + 1e-30) if step_converged: - z = dm.getGlobalVec() - z.set(0.0) - F0hat = rotated_residual(z) - fnorm0 = F0hat.norm() - F0hat.destroy() - dm.restoreGlobalVec(z) - step_converged = rnorm <= rtol * fnorm0 + atol + step_converged = rnorm <= rtol * ref + atol if not step_converged: mpi.pprint(f"[rotated_bc] step-norm at rel |F̂| = " - f"{rnorm / (fnorm0 + 1e-300):.2e} of the rest-state " + f"{rnorm / (ref + 1e-300):.2e} of the reference " f"residual — stagnation-or-crawl, continuing to iterate.") improved = False if lift_datum and iters == 0: @@ -801,10 +812,10 @@ def rotated_residual(uvec, keep_cartesian=False): # meeting the residual / step-norm criteria. Warn — as the standard SNES path does # on divergence — so an unconverged iterate left in the fields is not silent. if not converged: - rel = (rnorm / (r0 + 1e-300)) if r0 is not None else float("nan") + rel = (rnorm / (ref + 1e-300)) if ref is not None else float("nan") mpi.pprint(f"[rotated_bc] WARNING: nonlinear rotated free-slip did NOT converge " - f"in {newton_its} iterations (rel |F̂| = {rel:.2e}); the fields hold " - f"the last (unconverged) iterate.") + f"in {newton_its} iterations (rel |F̂| = {rel:.2e} of the reference " + f"residual); the fields hold the last (unconverged) iterate.") Fc.destroy() # residual output buffer (reaction persists in the result dict) _destroy_rotated_ksp_ctx(ctx) # KSP/PC + the owned Schur pmat From 3450db89769b3cb82d8f191ecc4d8a068549a684 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 16:27:05 +1000 Subject: [PATCH 06/17] docs(constitutive): TODO(DESIGN) - verify the TI consistent tangent with the PETSc Jacobian checker Regroup note (#438 B-session): TI power-law under consistent_jacobian=True shows an early trajectory identical to the frozen tangent and Picard-paced convergence (24 its vs ~11 isotropic) - consistent with the dC/d(eps_II) director-orientation terms being dropped or mis-oriented in the JIT-lowered tangent, while the isotropic terms survive (isotropic Newton is superlinear through the same rotated machinery, which also argues the rotation algebra itself is transparent). Discriminating harness (native -snes_test_jacobian vs rotated directional-derivative probe, with an isotropic control) is staged in the #438 acceptance run directory; not yet run - joint item. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 7d37eecf..1f539369 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -2474,6 +2474,17 @@ def flux(self): class TransverseIsotropicFlowModel(ViscousFlowModel): + # TODO(DESIGN): verify this model's CONSISTENT (Newton) tangent with the PETSc + # Jacobian checker (-snes_test_jacobian, NATIVE essential-BC path — no rotated + # machinery involved, so the test isolates the constitutive tangent). Evidence + # of a possible defect (#438 B-session, 2026-07-28): with a strain-rate- + # dependent eta_0/eta_1, consistent_jacobian=True produces an early iteration + # trajectory identical to the frozen (Picard) tangent and Picard-paced + # convergence (24 its vs ~11 for the isotropic power-law) — consistent with + # the dC/d(eps_II) orientation terms (the director e⊗e structure) being + # dropped or mis-oriented in the JIT-lowered g3, while the isotropic terms + # survive. Checker harness: ti_tangent_check.py in + # the #438 acceptance run directory (see the issue thread). r""" Transversely isotropic (anisotropic) viscous flow model. From 19c1fe153cee7f2c6606200e133d0971388e2249 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 17:05:21 +1000 Subject: [PATCH 07/17] docs(constitutive): TI consistent-tangent TODO upgraded to TODO(BUG) -> #457 The PETSc Jacobian checker (bounded-curvature protocol; at 1e-12 regularisation the FD reference is invalid for power-law and every model reads the same large mismatch) convicts the TI anisotropy tangent on the NATIVE path: ||J-Jfd||/||J|| ~ 2e-3-5e-3 in developed flow vs an FD-limited 4e-6 isotropic control. The rotated machinery is exonerated. Evidence and harness recorded in #457. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 1f539369..2bc5cf33 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -2474,17 +2474,15 @@ def flux(self): class TransverseIsotropicFlowModel(ViscousFlowModel): - # TODO(DESIGN): verify this model's CONSISTENT (Newton) tangent with the PETSc - # Jacobian checker (-snes_test_jacobian, NATIVE essential-BC path — no rotated - # machinery involved, so the test isolates the constitutive tangent). Evidence - # of a possible defect (#438 B-session, 2026-07-28): with a strain-rate- - # dependent eta_0/eta_1, consistent_jacobian=True produces an early iteration - # trajectory identical to the frozen (Picard) tangent and Picard-paced - # convergence (24 its vs ~11 for the isotropic power-law) — consistent with - # the dC/d(eps_II) orientation terms (the director e⊗e structure) being - # dropped or mis-oriented in the JIT-lowered g3, while the isotropic terms - # survive. Checker harness: ti_tangent_check.py in - # the #438 acceptance run directory (see the issue thread). + # TODO(BUG): the CONSISTENT (Newton) tangent of this model is inconsistent + # with its residual when the anisotropy is active (eta_1 != eta_0): PETSc + # -snes_test_jacobian at bounded curvature reads ||J-Jfd||/||J|| ~ 2e-3-5e-3 + # in developed flow where the isotropic control is FD-limited (~4e-6) — the + # dC/d(eps_II) director-coupled terms are missing or wrong, so TI "Newton" + # runs at Picard pace. Native-path evidence (no rotated machinery); the + # isotropic limit through this class is clean at rest but shares the defect + # once eta_1 differs. See issue #457 (incl. the checker-regularisation + # methodology: at 1e-12 the FD reference itself is invalid for power-law). r""" Transversely isotropic (anisotropic) viscous flow model. From cda37bc3d05ed84ed173c1ae58eaa65bf83a9501 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 28 Jul 2026 17:05:37 +1000 Subject: [PATCH 08/17] docs+test: rotated-freeslip subsystem page, changelog entry, np2 nonlinear datum test - docs/developer/subsystems/rotated-freeslip.md: the unified rotated path - constraint semantics, datum-through-Newton design (feasible iterates, cold-start lift), convergence reference scale, operator-verified null modes, parallel rules, reaction/sigma_nn hand-off. Linked from the developer index. - Changelog entry for the #438/#403 work at the conceptual level. - tests/parallel/test_1066: nonlinear (power-law) rotated datum at np>=2 - guards the collective-sequence class through the Newton loop (lift, line-search residuals, mode verification). Passes np2 and np4. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 49 +++++++ docs/developer/index.md | 2 + docs/developer/subsystems/rotated-freeslip.md | 130 ++++++++++++++++++ .../test_1066_rotated_datum_parallel.py | 41 ++++++ 4 files changed, 222 insertions(+) create mode 100644 docs/developer/subsystems/rotated-freeslip.md diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 1cdbc5dd..a80c6a63 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,55 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026) + +**Rotated strong free-slip now takes a prescribed wall-normal velocity datum +(`add_rotated_freeslip_bc(conds, boundary, ...)` with non-zero `conds`) through +the full nonlinear Newton machinery, and the separate linear and nonlinear +solve paths have been unified into one** (#438; #403 items 2 and 4). + +The rotated constraint `u·n̂ = ũ_n` is the primitive behind both the held free-slip +lid and the free-surface material-boundary condition — they differ only in the +constraint right-hand side. Previously the non-zero datum existed only on a +linear one-shot path, so a power-law or anisotropic rheology silently fell back +to a weak penalty (which leaks worst exactly where anisotropy makes it matter). +Now every rotated solve runs a single Newton/Picard loop in which accepted +iterates carry the datum exactly; a cold start imposes it through the first +increment's affine lift at the rest-state tangent (snapping a zero state onto a +datum creates a boundary strain state a shear-thinning tangent cannot recover +from — measured, not assumed); a linear model simply converges after that first +increment, so the up-front nonlinearity probe is gone from the dispatch and +every linear rotated solve saves two Jacobian assemblies. + +Three latent solver defects were found and fixed by making the loop report +honestly along the way: + +- Branching on rank-local datum bookkeeping desynchronised the ranks' + collective sequences (an np>1 deadlock class); the datum-activity decision is + now a collective PETSc reduction. +- A rigid-rotation mode pinned by an essential condition on *another* boundary + could still enter the solver null space (only the rotated rows were checked), + silently projecting an irreducible component out of every increment — Newton + converged superlinearly and then floored, far above tolerance. Candidate + modes are now verified as null vectors of the assembled operator, which + catches any form of pinning. +- A tiny Newton step was reported as convergence even when the residual was + still large (a stiff tangent also produces tiny steps); the step-norm exit is + now verified against the problem's rest-state residual scale, which is also + the convergence reference for warm starts (rtol relative to a good warm + start's own small initial residual demands ever-more absolute accuracy). + +The free-surface manager's `consistent_constraint="strong"` therefore works for +nonlinear rheologies: the penalty fallback is removed, and the consistent solve +warm-starts from the free solve's converged fields. Acceptance: power-law +annulus free-surface convection holds the material boundary at the strong-datum +level (5e-3, versus 4e-2 for the penalty) with net surface flux 1e-4; a +transversely isotropic fault-bearing smoke test converges through the same +path. Direct LU per Newton increment remains available as a serial, +preconditioner-free diagnostic (`solver._rotated_use_lu`). + +New subsystem documentation: `subsystems/rotated-freeslip.md`. + ### Local Interpolation That Reproduces Linear Fields (July 2026) **The local scattered-point interpolator now has a linear-reproduction diff --git a/docs/developer/index.md b/docs/developer/index.md index 977e09ac..7e9617d3 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -36,6 +36,7 @@ same topic are reference or historical material subordinate to the governing doc | Coding style & API conventions | [UW3 Style Charter](UW3_STYLE_CHARTER.md) (detailed reference: [Style Guide](UW3_Style_and_Patterns_Guide.md)) | | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) | | Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) | +| Rotated free-slip & wall-normal datum | [subsystems/rotated-freeslip.md](subsystems/rotated-freeslip.md) | | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) | | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) | | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) | @@ -184,6 +185,7 @@ subsystems/mesh-shape-relaxation subsystems/discretisation subsystems/solvers subsystems/boundary-stress-and-projection-postprocessing +subsystems/rotated-freeslip subsystems/petsc-jacobian-layout subsystems/constitutive-models subsystems/constitutive-models-theory diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md new file mode 100644 index 00000000..7e30897e --- /dev/null +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -0,0 +1,130 @@ +# Rotated strong free-slip and the prescribed-normal datum + +The rotated free-slip boundary condition imposes + +$$\mathbf{u}\cdot\hat{\mathbf{n}} = \tilde u_n$$ + +strongly — as an exact per-node constraint, not a weak (Nitsche/penalty) term — +by rotating each boundary node's velocity components into a (normal, +tangential) frame and constraining the rotated normal component. `conds = 0` +(the default datum) is free-slip; a non-zero scalar `conds` prescribes the +wall-normal velocity, which is what the free-surface manager uses to keep the +surface a material boundary. + +```python +stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) # free-slip +stokes.add_rotated_freeslip_bc(h_dot.sym[0], "Upper", normal=nhat) # u·n̂ = field +``` + +`normal=None` uses the geometric facet normal; a sympy `1×dim` matrix in +`mesh.X` supplies an analytic normal (exact `X/|X|` on curved boundaries — the +preferred choice there); a constant array is also accepted. The datum must be a +*scalar* (a number, an expression of `mesh.X`, or a scalar field read); on an +enclosed boundary it must be discretely flux-free for incompressibility. A +corner or 3D-edge node shared between rotated boundaries has no single normal +and stays at the free-slip pinning (the datum is ignored there). + +Why strong rather than Nitsche/penalty: the constraint holds to machine +precision (a penalty leaks ~1e-3, and the leak grows exactly where anisotropy +makes the boundary condition matter), it is correct on curved/tilted/deformed +boundaries, and the constraint **reaction is the boundary normal traction** +σ_nn (`solver.boundary_normal_traction`, `solver.dynamic_topography`) — the +quantity the free-surface machinery consumes. Reserve Nitsche for conditions +that must morph in time (Dirichlet→Neumann ramps). + +**Implementation**: `src/underworld3/utilities/rotated_bc.py`; registration and +dispatch in `petsc_generic_snes_solvers.pyx` (`add_rotated_freeslip_bc`). + +## One solve path + +There is a single driver, `solve_rotated_freeslip`: a manual outer +Newton/Picard loop that rotates the residual and tangent every iteration +(`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows, +and solves each increment with a self-contained fieldsplit-Schur KSP (geometric +FMG on the custom prolongation when a hierarchy is registered, else GAMG; the +native 1/μ pressure mass as the Schur preconditioner). There is no separate +linear path and no up-front nonlinearity probe: a linear model converges after +its first increment and the loop self-terminates. + +The manual loop exists because the rotated operator `Q A Qᵀ` carries no DM +field information, so PETSc's DM-coupled fieldsplit cannot precondition it; +the increment is solved by an IS-built fieldsplit instead, driven from the +loop. The loop honours `zero_init_guess`, `picard`, and the solver's +`consistent_jacobian` tangent policy (frozen / Newton / continuation). + +### How the datum passes through Newton + +The constraint is affine, and the design keeps every **accepted iterate +feasible** (`u·n̂ = ũ_n` exactly, via an affine snap in the rotated frame), so +each Newton increment satisfies the *homogeneous* constraint `n̂·δ = 0` — the +datum never touches the tangent, the increment right-hand side, the FMG +prolongation, or the null-space handling. + +The one deliberate exception is a **cold start with a non-zero datum**: the +first increment carries the datum jump through the affine lift +(`zeroRowsColumns(rows, diag, x̂, b̂)`) at the rest-state tangent, and is +accepted without a line search. Snapping the zero state onto the datum instead +manufactures an extreme boundary-strip strain state whose shear-thinning +tangent produces a first step orders of magnitude too large — no line-searched +step descends (measured). A warm start is assumed smooth and takes the exact +snap directly. + +### Convergence semantics + +The loop converges on `‖F̂‖ ≤ rtol·ref + atol` with +`ref = max(‖F̂(u₀)‖, ‖F̂(0)‖)`: the **rest-state residual is the intrinsic +forcing scale**, so a good warm start (whose own initial residual is small) is +not punished with an ever-stricter absolute target. A step-norm exit (tiny +Newton step) is *verified* against the same reference before it may report +convergence — a stiff tangent also produces tiny steps far from the solution, +and an unverified tiny step goes back through the line search (a slow crawl +continues; a genuine stall ends the loop with an explicit warning). An +unconverged exit always warns; the fields hold the last iterate. + +### Null-space handling + +Rigid-rotation candidates (one mode on a closed circle, three on a spherical +shell) are admitted to the increment null space and the post-solve gauge +removal only if they pass **two** tests: tangential to every rotated +constraint row (`Q·m ≈ 0` there), **and** a null vector of the assembled +operator (`‖J·m‖ ≈ 0` against the velocity diagonal scale). The second test is +what catches pinning the first cannot see — an essential condition on another +boundary leaves the mode tangential to the rotated wall while the eliminated +operator is *not* null on it, and admitting it projects an irreducible +component out of every increment (the residual then floors far above +tolerance instead of converging). The null space is therefore built after the +first Jacobian assembly. + +### Parallel notes + +- Every branch the loop takes is decided **collectively**. The datum-activity + flag rides a PETSc `Vec` norm of the lift vector (the datum bookkeeping + itself is rank-local — only owners of datum boundary nodes hold entries — + and branching on it per-rank desynchronises the collective sequences: the + np>1 deadlock class). +- All row surgery on vectors uses ownership-relative indices + (`_zero_rows_local` / `_set_rows_local`); indexing a local slice with global + rows is the np>1 crash class. +- Direct LU per increment (`solver._rotated_use_lu = True`) is a **serial** + preconditioner-free diagnostic (the pressure-gauge pin is a naive per-rank + scan, marked `TODO(BUG)`); use it to separate "the operator/constraint is + wrong" from "the preconditioner is struggling". + +## Result and reaction + +The solve fills the solver's fields and stores a result dict +(`solver._rotated_freeslip_info`) with the rotation, the constrained rows, the +per-increment KSP iteration counts, the convergence verdict, and the +**reaction** — the converged Cartesian residual `F(u)`, which for a linear +residual equals `A·u − b` exactly. `boundary_normal_traction` projects the +nodal reaction onto the boundary normal (corner-correct) and de-smears it with +the shared boundary-mass machinery in `utilities/boundary_flux.py`; +`dynamic_topography_field` writes `h = −(σ_nn − σ̄_nn)/(Δρ g)` onto a surface +field for the free-surface integrator. + +## Tests + +`tests/test_1018_rotated_freeslip.py` (serial: essential-equivalence, FMG, +tangent policies, datum linear + nonlinear), +`tests/parallel/test_1066_rotated_datum_parallel.py` (np≥2: partition +independence of the linear datum and the nonlinear Newton datum path). diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index b86a6b79..f15e9ef8 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -56,3 +56,44 @@ def test_rotated_datum_prescribed_normal_partition_independent(): vv = float(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate()) assert abs(vv - _INT_VV_REF) / _INT_VV_REF < 1.0e-6, \ f"solve energy is partition-dependent: ∫v·v={vv:.8e} vs ref {_INT_VV_REF}" + + +def test_rotated_datum_nonlinear_parallel(): + """The NONLINEAR rotated datum path in parallel: a power-law annulus with + u.n = cos(theta) through the manual Newton loop (feasible iterates, cold-start + lift, collective datum-activity flag). Guards the np>1 collective-sequence + class: the lift branch, the line-search residual evaluations and the + nullspace-mode verification are all collective, and a rank-divergent branch + deadlocks rather than failing cleanly (hence the module timeout).""" + RI, RO = 0.5, 1.0 + mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x ** 2 + y ** 2) + nhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable("Vnl64", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable("Pnl64", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + gm = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + e = 0.5 * (gm + gm.T) + eII = sympy.sqrt(0.5 * (e[0, 0] ** 2 + e[1, 1] ** 2) + e[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) + s.consistent_jacobian = True + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-7 + s.solve() + + info = s._rotated_freeslip_info + assert info["converged"], "nonlinear parallel datum solve did not converge" + assert info["nonlinear_iterations"] > 1, "did not genuinely iterate" + # datum imposed, measured with parallel-safe boundary quadrature only + vn = (v.sym[0] * x + v.sym[1] * y) / r + num = float(uw.maths.BdIntegral(mesh=mesh, fn=(vn - x / r) ** 2, boundary="Upper").evaluate()) + den = float(uw.maths.BdIntegral(mesh=mesh, fn=(x / r) ** 2, boundary="Upper").evaluate()) + relL2 = np.sqrt(max(num, 0.0) / den) + assert relL2 < 1.0e-3, f"nonlinear u.n=cos(theta) not imposed (relL2={relL2:.3e})" From 4a557b81d8fe3a7e7f51fa5519d96d756d055e71 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 12:48:32 +1000 Subject: [PATCH 09/17] fix(rotated-bc): the custom-FMG velocity sub-KSP must converge, not just apply The rotated free-slip path drove its custom-FMG velocity block with `preonly` ("a full-MG cycle per Schur application, by design"). PCFieldSplit forms the Schur complement S = A11 - A10 A00^-1 A01 and applies A00^-1 through that same velocity KSP, so `preonly` replaces A00^-1 with a single multigrid cycle and hands the pressure Krylov a different system S~ != S, preconditioned by a 1/mu mass matrix built for S. Measured (annulus, weak plane reaching the constrained boundary, eta_1/eta_0 = 1e-3, transversely isotropic): under `preonly` the pressure residual falls 4.4e4 in ~16 iterations then STAGNATES at a floor ~3.1e-7 and burns the remaining 184 iterations of its cap for nothing, every outer iteration -- 9 outer iterations, 2.77 s. Under FGMRES it converges 1.1e8 monotonically in 17 pressure iterations -- 1 outer, 0.67 s (4.1x). The isotropic control moves the same way (5 -> 1), so this is the Schur application and not the anisotropy. SolCx box + 2-level barycentric hierarchy (the test_1018 FMG configuration): 6 outer -> 1. The Schur application is bitwise reproducible under both settings (measured), so the floor is NOT "the operator changes between applications" -- S~ is a fixed linear operator, just the wrong one. The precise origin of the floor is not isolated. Wrap the FMG in FGMRES: max_it 200 matches the GAMG fallback and the native path, rtol 0.1 x tol matches the GAMG fallback (the native path asks 0.033 x tol). Found while auditing whether the rotated path preconditions the ROTATED operator. It does: an algebraic identity probe against independently formed references (block the native Jacobian, then rotate -- the opposite order to the solver) shows the velocity sub-KSP operator, the velocity sub-PC operator and every Galerkin multigrid level operator match the rotated reference exactly, while differing from the un-rotated reference by 7e-2 to 2.3e-1. The Schur preconditioner correctly comes un-rotated from the native Pmat: Q is the identity on pressure, and the p-p block of the rotated operator is identically zero. No un-rotated matrix reaches the preconditioner. Also: - warn when the velocity sub-KSP exhausts its cap. It returns KSP_DIVERGED_ITS, which KSPCheckSolve deliberately does not escalate, so FGMRES can degrade silently where `preonly` could not fail. - report `velocity_pc`, `schur_pre`, `velocity_pc_type` (PETSc's own view of the sub-PC) and `vel_its_last` / `pres_its_last` in `_rotated_freeslip_info`. The outer count alone hides a degraded inner solve, which is how this went unnoticed. The `_last` names are deliberate: these are last-application samples, not the summed work axis solver_health uses -- wiring the rotated path into SolverInstrumentation.sub_reports() is the proper fix and is not done here. - guard the FMG configuration in test_1018 against a return to `preonly`, and assert the velocity sub-PC really is multigrid rather than restating our own bookkeeping. - fix a pre-existing failure in test_1064: `int(info["ksp_its"])` on a value the unified rotated loop made a list, so the 3D spherical partition-independence test errored at np>1 before reaching its assertions. Tests: test_1018 18 passed; test_1064 + test_1066 10 passed at np2 and np4; TI fault smoke passes (24 nonlinear its, datum relL2 4.4e-05). Underworld development team with AI support from Claude Code --- docs/developer/subsystems/rotated-freeslip.md | 48 ++++++++++++ src/underworld3/utilities/rotated_bc.py | 77 +++++++++++++++++-- .../test_1064_rotated_freeslip_parallel.py | 3 +- tests/test_1018_rotated_freeslip.py | 20 ++++- 4 files changed, 139 insertions(+), 9 deletions(-) diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 7e30897e..14737fe1 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -46,6 +46,54 @@ native 1/μ pressure mass as the Schur preconditioner). There is no separate linear path and no up-front nonlinearity probe: a linear model converges after its first increment and the loop self-terminates. +### The velocity sub-KSP must converge, not just apply + +The velocity block is preconditioned by multigrid, but the KSP *around* that +multigrid is FGMRES to `0.1 × tolerance` — never `preonly`. This is not a +tuning preference. `PCFieldSplit` forms the Schur complement +`S = A₁₁ − A₁₀ A₀₀⁻¹ A₀₁` and applies `A₀₀⁻¹` through this same velocity KSP, +so `preonly` replaces `A₀₀⁻¹` with a single multigrid cycle and hands the +pressure Krylov a *different* system `S̃ ≠ S`, preconditioned by a 1/μ mass +matrix built for `S`. + +Measured on an annulus with a weak plane reaching the constrained boundary +(`η₁/η₀ = 1e-3`, transversely isotropic): under `preonly` the pressure residual +falls by 4.4e4 in about 16 iterations and then **stagnates at a floor ≈ 3.1e-7**, +spending the remaining 184 iterations of its cap for nothing — on every outer +iteration, 9 outer iterations in total. Under FGMRES it converges by 1.1e8 +monotonically in 17 pressure iterations, and the outer solve takes 1. The +isotropic control moves the same way (5 → 1), so this is a property of the Schur +application rather than of the anisotropy. + +Note the Schur application is bitwise reproducible under both settings, so the +floor is *not* "the operator changes between applications" — `S̃` is a fixed +linear operator, just the wrong one. The precise origin of the floor (most +plausibly a range/null-space inconsistency between `S̃` and the constant-pressure +null space attached to `S`) has not been isolated. + +The native (non-rotated) Stokes path and the GAMG fallback have always wrapped +their multigrid this way; the rotated custom-FMG branch was the sole exception. +`max_it` matches both (200); `rtol` matches the GAMG fallback (0.1 × tol), where +the native path asks for 0.033 × tol. + +Two things make this easy to miss, and both are now instrumented. The outer +iteration count does not reveal it — a full Schur factorisation with a good +pressure mass still reports ~1 outer iteration while the inner solve grinds +underneath. And a velocity FGMRES that exhausts its cap returns +`KSP_DIVERGED_ITS`, which `KSPCheckSolve` deliberately does not escalate, so it +would degrade silently where `preonly` simply could not fail; the rotated path +now warns on it. + +`solver._rotated_freeslip_info` therefore reports `velocity_pc`, `schur_pre` and +`velocity_pc_type` (PETSc's own view of the sub-PC), plus `vel_its_last` and +`pres_its_last`. The last two are **last-application samples, not work**: +`KSPGetIterationNumber` reports only the most recent solve, and the velocity KSP +is applied once per Schur `MatMult`. They are named `_last` so they are not +confused with the summed counts `solver_health` uses as its work axis; wiring +the rotated path into `SolverInstrumentation.sub_reports()` is the proper fix and +has not been done. On the opt-in `solver._rotated_use_lu` path there are no +sub-KSPs at all and both fields record `None`. + The manual loop exists because the rotated operator `Q A Qᵀ` carries no DM field information, so PETSc's DM-coupled fieldsplit cannot precondition it; the increment is solved by an IS-built fieldsplit instead, driven from the diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 448a7297..856f0a44 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -644,6 +644,10 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, ctx = None diag_scale = None lin_its = [] + # velocity / pressure sub-KSP LAST-APPLICATION counts, one per Newton increment + # (iterative path only — the direct-LU path has no sub-KSPs and records None) + vel_its_last = [] + pres_its_last = [] def rotated_residual(uvec, keep_cartesian=False): snes.computeFunction(uvec, Fc) @@ -752,6 +756,8 @@ def rotated_residual(uvec, keep_cartesian=False): solver, Ahat, bhat, Q, Qt, normal_rows, custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) lin_its.append(ctx["ksp"].getIterationNumber()) + vel_its_last.append(ctx.get("vel_its_last")) + pres_its_last.append(ctx.get("pres_its_last")) d = dm.createGlobalVec() Qt.mult(dhat, d) # step-norm convergence (SNES_CONVERGED_SNORM): a tiny Newton step means we @@ -830,6 +836,12 @@ def rotated_residual(uvec, keep_cartesian=False): "rotation_gauge_removed": removed, "ksp_reason": last_reason, "nonlinear_iterations": newton_its, "converged": converged, "ksp_its": lin_its, "rnorm": rnorm, "rnorm0": r0, + "vel_its_last": vel_its_last, "pres_its_last": pres_its_last, + # keyed on use_lu, NOT on `ctx is None`: the iterative path also exits + # with ctx unset when a warm start is already converged at increment 0. + "velocity_pc": "direct-LU" if use_lu else (ctx or {}).get("velocity_pc"), + "schur_pre": "none" if use_lu else (ctx or {}).get("schur_pre"), + "velocity_pc_type": None if use_lu else (ctx or {}).get("velocity_pc_type"), "continuation_switched": continuation and phase == "newton"} @@ -979,8 +991,37 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal "fieldsplit_vel_mg_levels_ksp_converged_maxits": "true", }) else: - # full-MG cycle per Schur application, by design - cfg["fieldsplit_vel_ksp_type"] = "preonly" + # Custom-FMG velocity block, wrapped in a short FGMRES — NOT `preonly`. + # PCFieldSplit applies the Schur complement S = A11 - A10 A00^-1 A01 + # through THIS velocity KSP, so `preonly` replaces A00^-1 with a single + # multigrid cycle and the pressure Krylov is handed a different system + # S~ != S, preconditioned by a 1/mu mass built for S. + # + # Measured (annulus, weak plane reaching the constrained boundary, + # eta_1/eta_0 = 1e-3, transversely isotropic): under `preonly` the + # pressure residual falls 4.4e4 in ~16 iterations and then STAGNATES at + # a floor ~3.1e-7, burning the remaining 184 iterations of its cap for + # nothing, every outer iteration — 9 outer iterations total. Under + # FGMRES it converges 1.1e8 monotonically in 17 pressure iterations — 1 + # outer. The isotropic control moves the same way (5 -> 1), so this is + # the Schur application and not the anisotropy. + # + # NB the Schur application is bitwise reproducible under both settings + # (measured), so the floor is NOT "the operator changes between + # applications" — it is S~ being the wrong system. The precise origin of + # the floor (most likely a range/nullspace inconsistency between S~ and + # the constant-pressure null space attached to S below) is not isolated. + # + # max_it matches the GAMG fallback and the native path (200); rtol + # matches the GAMG fallback (0.1 x tol). The native path asks for + # 0.033 x tol — deliberately not copied, since the FMG cycle is a far + # stronger preconditioner than GAMG here and reaches 0.1 x tol in ~11 + # iterations. + cfg.update({ + "fieldsplit_vel_ksp_type": "fgmres", + "fieldsplit_vel_ksp_rtol": str(tol * 0.1), + "fieldsplit_vel_ksp_max_it": "200", + }) for k, v in cfg.items(): opts[pfx + k] = v try: @@ -1032,7 +1073,11 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal except Exception: pass ctx = {"ksp": ksp, "pc": pc, "Mp": Mp, "nsp": nsp, "cns": cns, - "custom_Pl": custom_Pl, "pfx": pfx} + "custom_Pl": custom_Pl, "pfx": pfx, + # which preconditioner this solve actually got — reported so a model + # (or a test) can assert on it instead of inferring it from timings + "velocity_pc": "custom-FMG" if custom_Pl is not None else "GAMG", + "schur_pre": "1/mu-mass" if Mp is not None else "selfp"} else: ksp = ctx["ksp"] nsp = ctx["nsp"] @@ -1053,12 +1098,30 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # equation → setting them to exactly 0 here makes the strong v_n=0 BC exact # independent of the iterative tolerance, without perturbing the rest. _zero_rows_local(Uhat, normal_rows) + # Sub-KSP diagnostics: the OUTER count alone hides a degraded inner solve — a + # full Schur factorisation with a good pressure mass converges in ~1 outer + # iteration while the pressure sub-solve grinds against its cap underneath. + # + # These are LAST-APPLICATION samples, not work: KSPGetIterationNumber reports + # the most recent solve, and the velocity KSP is applied once per Schur MatMult + # (i.e. once per pressure Krylov iteration). They are named `_last` so they are + # not mistaken for the summed counts the solver_health report uses as its work + # axis — wiring the rotated path into SolverInstrumentation.sub_reports() is the + # proper fix and is not done here. + vel_ksp, pres_ksp = ctx["pc"].getFieldSplitSubKSP() + ctx["vel_its_last"] = vel_ksp.getIterationNumber() + ctx["pres_its_last"] = pres_ksp.getIterationNumber() + ctx["velocity_pc_type"] = vel_ksp.getPC().getType() + # A velocity FGMRES that exhausts its cap returns KSP_DIVERGED_ITS, which + # KSPCheckSolve deliberately does NOT escalate — so without this it degrades + # silently. (`preonly` could never fail, which is why the check is new.) + _warn_if_ksp_diverged(vel_ksp, kind="rotated fieldsplit velocity sub-solve") if verbose: - kind = "custom-FMG" if ctx["custom_Pl"] is not None else "GAMG" - schur = "1/mu-mass" if ctx["Mp"] is not None else "selfp" - mpi.pprint(f"[rotated_bc] velocity block = {kind}; Schur pre = {schur}; " + mpi.pprint(f"[rotated_bc] velocity block = {ctx['velocity_pc']} " + f"(PC {ctx['velocity_pc_type']}); Schur pre = {ctx['schur_pre']}; " f"outer KSP {ksp.getConvergedReason()} in " - f"{ksp.getIterationNumber()} its") + f"{ksp.getIterationNumber()} its (last apply: vel " + f"{ctx['vel_its_last']}, pres {ctx['pres_its_last']})") return Uhat, ksp.getConvergedReason(), ctx diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 576d5340..61c630ef 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -166,7 +166,8 @@ def _spherical3d_diagnostics(): L2 = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) info = s._rotated_freeslip_info - return L2, int(info["ksp_its"]), int(info["ksp_reason"]) + # the unified rotated loop reports one KSP count per Newton increment + return L2, max(info["ksp_its"]), int(info["ksp_reason"]) def _spherical3d_topography_diagnostics(cell_size=0.25): diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 02d81ada..75a42ea9 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -279,10 +279,28 @@ def test_rotated_freeslip_geometric_fmg_velocity_block(): s.solve() # geometric MG on the velocity block converged, and matches essential - assert s._rotated_freeslip_info["ksp_reason"] > 0 + info = s._rotated_freeslip_info + assert info["ksp_reason"] > 0 rel = np.linalg.norm(v.data - vE.data) / np.linalg.norm(vE.data) assert rel < 5e-3, f"FMG rotated free-slip differs from essential by {rel:.2e}" + # The velocity block really IS multigrid (PETSc's own view of the sub-PC, not our + # bookkeeping), preconditioned by the 1/mu pressure mass. + assert info["velocity_pc"] == "custom-FMG" + assert info["schur_pre"] == "1/mu-mass" + assert info["velocity_pc_type"] == "mg", ( + f"velocity sub-PC is {info['velocity_pc_type']!r}, not multigrid — the custom " + f"FMG install did not take") + + # PCFieldSplit applies the Schur complement through the velocity sub-KSP, so that + # KSP must converge to a tolerance. A `preonly` single multigrid cycle hands the + # pressure Krylov a different system, which stagnates above its tolerance and + # leaves the outer Krylov to make up the difference. Measured on this + # configuration: FGMRES-wrapped FMG = 1 outer iteration, `preonly` = 6. + assert max(info["ksp_its"]) <= 3, ( + f"rotated FMG outer iteration blow-out: {info['ksp_its']} " + f"(inexact Schur application — is the velocity sub-KSP `preonly` again?)") + def test_rotated_freeslip_boundary_normal_traction_solcx(): """sigma_nn recovered by boundary_normal_traction on the top boundary reproduces From 6987211d43e4eb2112dd58d8d42d864dc6a442b0 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 13:58:13 +1000 Subject: [PATCH 10/17] solvers: adopt native sub-solve tolerances on the rotated path; document the Krylov axes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The rotated path was looser than the native Stokes path on both fieldsplit sub-solves, with no recorded reason: rotated (was) native velocity rtol 0.1 x tol 0.033 x tol pressure rtol 1.0 x tol 0.1 x tol House policy is that defaults err on the side of robust generality, with callers loosening guardrails at their own risk. Matching native costs ~17% wall clock with IDENTICAL outer iteration counts, measured on both velocity-block routes (custom FMG and the GAMG fallback) and both rheologies (isotropic and transversely isotropic). It also drops the TI fault smoke from 24 to 15 nonlinear iterations — tighter inner solves feed a cleaner Newton step. A too-loose inner solve does not degrade gracefully: it walks back toward the `preonly` failure this branch already fixed (velocity rtol 1e-1 -> 6 outer iterations, 1e-2 -> 4, 0.033 x tol -> 1). Cheaper remains available to anyone who measures their own configuration; it is not the default. Also documents the three axes that keep getting argued as one, in the solvers subsystem page (which listed "preconditioner selection missing" as a known gap): 1. Flexible vs non-flexible Krylov. The question is whether the PRECONDITIONER is stationary, not whether the operator is symmetric — the velocity block is SPD, so cg is admissible on symmetry grounds and still fails, because GAMG with mg_levels_ksp_converged_maxits is a non-linear application. Settled by #147. 2. `preonly` vs an iterative wrapper. Not a Krylov-taste question but a positional one: a preconditioner application may be `preonly`, an operator inverse underneath a Schur complement may not, because it defines the system the pressure Krylov solves rather than merely conditioning it. 3. Inner tolerance. The genuine judgment call, decided here by the guardrail policy. Plus how to detect a degraded sub-solve: the outer count hides it, an exhausted cap returns KSP_DIVERGED_ITS which KSPCheckSolve does not escalate, and KSPGetIterationNumber on a sub-KSP is a last-application sample rather than work. Tests: test_1018 18 passed; test_1064 + test_1066 10 passed at np2 and np4; TI fault smoke passes (15 nonlinear its, datum relL2 4.4e-05). Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 115 +++++++++++++++++++++++- src/underworld3/utilities/rotated_bc.py | 25 +++--- 2 files changed, 128 insertions(+), 12 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index e6710bbb..bc702b94 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -175,9 +175,120 @@ def validate_unknowns_sharing(multi_material_model): assert model.Unknowns.DFDt is reference_unknowns.DFDt, \ f"Model {i} $D\mathbf{{F}}/Dt$ not shared - stress history will be wrong" ``` -- ⚠️ Preconditioner selection missing +- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" below; multigrid and Schur preconditioner choice is still undocumented - Could benefit from optimization examples +## Choosing the Krylov method for a fieldsplit sub-solve + +This choice gets re-argued periodically because **three independent questions are +easily merged into one**. Two of them are settled, by different evidence and for +different reasons; only the third is a genuine judgment call. + +```{note} Guardrail policy +Defaults err on the side of robust generality. A default that is slower but +survives configurations nobody has tested yet is the right default; loosening it +is the caller's decision and the caller's risk. The three axes below are stated +in those terms. +``` + +### Axis 1 — flexible or not (`fgmres` vs `gmres` / `cg` / `fcg`) + +The question is **is the preconditioner stationary?**, not whether the operator is +symmetric. + +The velocity block of Stokes *is* SPD, so `cg` is admissible on symmetry grounds. +It still fails, because GAMG with `mg_levels_ksp_converged_maxits` performs a +variable number of smoothing iterations, making the preconditioner application +non-linear. `cg`/`fcg` residual recurrences and standard `gmres` cannot accommodate +that. FGMRES is right-preconditioned by construction and can. + +Settled empirically in [#147](https://github.com/underworldcode/underworld3/issues/147) +(spherical Kramer `case1`, free-slip Nitsche, Gadi at np=144, cellsize 1/32): +`fcg` reported an indefinite matrix / `DIVERGED_PC_FAILED`; `gmres` failed with a +residual-recursion mismatch; the same configuration converged on macOS, so this is +a robustness cliff exposed by scale, not a formulation error. The rationale is +recorded inline at `petsc_generic_snes_solvers.pyx` (velocity sub-solve block). + +**Default: `fgmres` on both sub-solves.** Anything non-flexible is only safe if you +can show your preconditioner is stationary, and ours generally is not. + +### Axis 2 — `preonly` or an iterative wrapper + +This is not a Krylov-taste question at all. It is **positional**: is this KSP a +*preconditioner application*, or an *operator inverse*? + +- **Top-level PC** (scalar and single-field vector solvers, empty + `_pc_option_prefix`): the multigrid *is* the preconditioner and an outer Krylov + cleans up after it. One cycle per application is a legitimate design. +- **Velocity block under `pc_fieldsplit_schur_fact_type=full`**: `PCFieldSplit` + forms `S = A₁₁ − A₁₀ A₀₀⁻¹ A₀₁` and applies `A₀₀⁻¹` *through this KSP*. It is no + longer preconditioning anything — it **defines the operator the pressure Krylov + iterates against**. `preonly` there does not degrade conditioning; it changes + which system is being solved, to `S̃ ≠ S`, while the 1/μ pressure mass still + preconditions `S`. + +The failure is quiet. Measured on the rotated free-slip path (annulus, transversely +isotropic, weak plane reaching the constrained boundary, η₁/η₀ = 1e-3), the pressure +residual under `preonly` falls 4.4e4 in ~16 iterations and then **stagnates at a +floor ≈ 3.1e-7**, burning the remaining 184 iterations of its cap on every outer +iteration: 9 outer iterations and 2.77 s, against 1 outer, 17 pressure iterations +and 0.67 s once the same multigrid is wrapped in FGMRES. The isotropic control moves +identically (5 → 1 outer), so this is the Schur application and not the anisotropy. + +Note the operator is *not* moving between applications — applying the Schur operator +twice to the same vector is bitwise identical under both settings. `S̃` is a fixed +linear operator, just the wrong one. (Why the floor sits where it does is not +isolated; a range/null-space inconsistency between `S̃` and the constant-pressure +null space attached to `S` is the obvious suspect.) + +**Rule: `preonly` is fine as a preconditioner application and never fine as an +inverse underneath a Schur complement.** Multigrid is an excellent preconditioner; +it cannot be the whole solve when a Schur complement is applied through it. + +### Axis 3 — how tight the inner tolerance should be + +This is the axis where judgment legitimately lives, and where the guardrail policy +decides it. Loosening the velocity sub-solve does not degrade gracefully — it walks +back toward the `preonly` failure above: + +| velocity sub-KSP rtol | outer its | +|---|---| +| 1e-1 | 6 | +| 1e-2 | 4 | +| 0.033 × tolerance | 1 | + +Current defaults, matched across the native and rotated paths: + +| sub-solve | rtol | max_it | +|---|---|---| +| velocity | `0.033 × tolerance` | 200 | +| pressure | `0.1 × tolerance` | 200 | + +The rotated path previously used `0.1 × tolerance` (velocity) and `1.0 × tolerance` +(pressure) — 3× and 10× looser than native, with no recorded reason. Adopting the +native values costs ~17% wall clock with identical outer iteration counts (measured +on both velocity-block routes, isotropic and TI) and reduced a transversely +isotropic fault smoke test from 24 to 15 nonlinear iterations. Cheaper is available +to anyone who measures their own configuration; it is not the default. + +### Detecting a degraded sub-solve + +Two properties make this family of problems hard to see, and both need instrumenting +rather than eyeballing: + +- **The outer iteration count does not reveal it.** A full Schur factorisation with a + good pressure mass still reports ~1 outer iteration while the inner solve grinds + against its cap underneath. Read the sub-KSP counts. +- **An exhausted cap does not raise.** A sub-KSP that hits `max_it` returns + `KSP_DIVERGED_ITS`, which `KSPCheckSolve` deliberately does not escalate, so it + degrades silently. `preonly` could never fail this way, so switching to an + iterative wrapper introduces a failure mode that must be checked for explicitly. + +Beware also that `KSPGetIterationNumber` on a sub-KSP reports only its **most recent** +application, and the velocity KSP is applied once per Schur `MatMult` — i.e. once per +pressure Krylov iteration. That number is a sample, not work. `SolverInstrumentation` +(`systems/solver_health.py`) is the mechanism that sums properly. + ## Critical Stability Note ```{warning} Solver Stability is Paramount @@ -188,7 +299,7 @@ def validate_unknowns_sharing(multi_material_model): ```{note} For Contributors This well-documented subsystem could benefit from: -- Preconditioner selection guidance +- Preconditioner selection guidance (Krylov sub-solve choice is now covered; multigrid and Schur preconditioner selection is not) - Performance tuning documentation - Convergence analysis examples - Scaling studies and optimization diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 856f0a44..9edc867f 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -963,10 +963,15 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal "ksp_type": "fgmres", "ksp_rtol": str(tol), "ksp_max_it": "300", "pc_type": "fieldsplit", "pc_fieldsplit_type": "schur", "pc_fieldsplit_schur_fact_type": "full", - # native-parity pressure sub-solve (pyx Stokes defaults): FGMRES at the - # solver tolerance; GASM on the 1/mu mass, jacobi if only selfp exists. + # Pressure sub-solve at native-path parity (pyx Stokes defaults): + # FGMRES at 0.1 x tol, GASM on the 1/mu mass (jacobi if only selfp + # exists). This path used to ask for tol, i.e. 10x looser than native. + # Defaults here err toward the more conservative of the two: matching + # native costs ~17% wall clock with identical outer iteration counts + # (measured, both velocity-block routes, isotropic and TI), and a + # too-loose inner solve fails by silent stagnation rather than loudly. "fieldsplit_pres_ksp_type": "fgmres", - "fieldsplit_pres_ksp_rtol": str(tol), + "fieldsplit_pres_ksp_rtol": str(tol * 0.1), "fieldsplit_pres_ksp_max_it": "200", } if Mp is not None: @@ -980,7 +985,7 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # block — this applies only when no hierarchy is registered. cfg.update({ "fieldsplit_vel_ksp_type": "fgmres", - "fieldsplit_vel_ksp_rtol": str(tol * 0.1), + "fieldsplit_vel_ksp_rtol": str(tol * 0.033), "fieldsplit_vel_ksp_max_it": "200", "fieldsplit_vel_pc_type": "gamg", "fieldsplit_vel_pc_gamg_type": "agg", @@ -1012,14 +1017,14 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # the floor (most likely a range/nullspace inconsistency between S~ and # the constant-pressure null space attached to S below) is not isolated. # - # max_it matches the GAMG fallback and the native path (200); rtol - # matches the GAMG fallback (0.1 x tol). The native path asks for - # 0.033 x tol — deliberately not copied, since the FMG cycle is a far - # stronger preconditioner than GAMG here and reaches 0.1 x tol in ~11 - # iterations. + # rtol and max_it are the native path's (0.033 x tol, 200), as is the + # GAMG fallback above. The FMG cycle reaches 0.1 x tol in ~11 iterations + # and would be cheaper there, but a sub-solve tolerance is a guardrail: + # the conservative value is the default and loosening it is the caller's + # risk to take. cfg.update({ "fieldsplit_vel_ksp_type": "fgmres", - "fieldsplit_vel_ksp_rtol": str(tol * 0.1), + "fieldsplit_vel_ksp_rtol": str(tol * 0.033), "fieldsplit_vel_ksp_max_it": "200", }) for k, v in cfg.items(): From 3ade87c8f14f7dff6f93cf29f74771765e4f5e5c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 14:57:53 +1000 Subject: [PATCH 11/17] docs(solvers): the inner-tolerance margin and the flexible outer Krylov are one design MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Provenance and framing from Louis. The Stokes configuration descends from the Citcom solver of Moresi & Solomatov (1995), whose central choice is that the inner solves are DELIBERATELY inexact. Two consequences follow and must hold together: 1. Inexact inner solves perturb the search directions, so the outer/Schur Krylov must be FLEXIBLE (fgmres or similar) — a non-flexible method's residual recurrence is invalidated by a search direction that drifts. 2. Inexact is not unbounded: an inner solve must still converge WELL BELOW the tolerance demanded of the outer solve. The 0.033 and 0.1 factors are that margin — a safety margin, not a tuned constant. So what looked like two independent axes (flexible-or-not, how-tight) is one decision: flexibility buys tolerance of inexactness, the margin bounds how inexact. Arguing them separately is why the question never settled. `preonly` is then the degenerate case — no tolerance at all, so no margin for the flexible outer Krylov to work with, which is exactly the stagnation floor measured on this branch. This also sharpens the previous commit's justification. The rotated path ran its outer KSP at rtol = tol and its pressure sub-solve ALSO at rtol = tol: the inner solve permitted to be no better than the answer it feeds. That is the invariant broken outright, not a tuning inconsistency with the native path. Recorded in the solvers subsystem page and at the point in rotated_bc.py where the margin is set. The size of the factors is noted as inherited convention with no derivation recorded beyond "well below the outer tolerance". Tests: test_1018 18 passed. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 70 +++++++++++++++++++------ src/underworld3/utilities/rotated_bc.py | 18 +++++-- 2 files changed, 67 insertions(+), 21 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index bc702b94..32eaff42 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -180,21 +180,41 @@ def validate_unknowns_sharing(multi_material_model): ## Choosing the Krylov method for a fieldsplit sub-solve -This choice gets re-argued periodically because **three independent questions are -easily merged into one**. Two of them are settled, by different evidence and for -different reasons; only the third is a genuine judgment call. +This choice gets re-argued periodically. The confusion is that it looks like three +separate questions — flexible or not, `preonly` or not, how tight — when in fact the +first and third are **two halves of one design decision**, and only the middle one is +independent. + +### The design, and where it comes from + +The Stokes configuration descends from the Citcom solver of Moresi & Solomatov +(1995). Its central choice is that the **inner solves are deliberately inexact** — +you do not solve the velocity block exactly to apply the Schur complement, because +that would be ruinous and is unnecessary. Two consequences follow, and they must +hold together: + +1. **Because the inner solves are inexact, the search directions are perturbed, so + the outer/Schur Krylov must be flexible** (`fgmres` or similar). A non-flexible + method assumes a fixed preconditioner; its residual recurrence is invalidated by + a search direction that drifts. +2. **Inexact is not unbounded.** The inner solves must still converge to *well + below* the final tolerance required of the outer solve. That margin is what the + `0.033` and `0.1` factors encode — they are a safety margin, not a tuned constant. + +Flexibility buys tolerance of inexactness; the margin bounds how inexact. Neither +works without the other, which is why arguing them separately never settles. ```{note} Guardrail policy Defaults err on the side of robust generality. A default that is slower but -survives configurations nobody has tested yet is the right default; loosening it -is the caller's decision and the caller's risk. The three axes below are stated -in those terms. +survives configurations nobody has tested yet is the right default; loosening it is +the caller's decision and the caller's risk. ``` ### Axis 1 — flexible or not (`fgmres` vs `gmres` / `cg` / `fcg`) -The question is **is the preconditioner stationary?**, not whether the operator is -symmetric. +Settled, and settled by the design above: the inner solves are inexact by +construction, so the outer and Schur Krylov methods must be flexible. The question +is **is the preconditioner stationary?**, not whether the operator is symmetric. The velocity block of Stokes *is* SPD, so `cg` is admissible on symmetry grounds. It still fails, because GAMG with `mg_levels_ksp_converged_maxits` performs a @@ -245,10 +265,23 @@ null space attached to `S` is the obvious suspect.) inverse underneath a Schur complement.** Multigrid is an excellent preconditioner; it cannot be the whole solve when a Schur complement is applied through it. -### Axis 3 — how tight the inner tolerance should be +In the terms of the design above, `preonly` is the degenerate case of axis 3: it has +no tolerance at all, so there is no margin below the outer solve for the flexible +outer Krylov to work with. Flexibility tolerates inexactness; it cannot manufacture +a margin that was never there. That is exactly the stagnation floor measured above. -This is the axis where judgment legitimately lives, and where the guardrail policy -decides it. Loosening the velocity sub-solve does not degrade gracefully — it walks +### Axis 3 — how far below the outer tolerance the inner solves must go + +Not a free parameter. The invariant from the design above is that **every inner +solve reaches well below the tolerance demanded of the outer solve**; the only +judgment is how much margin, and the guardrail policy decides that. + +The failure this catches is concrete. The rotated free-slip path ran its outer KSP +at `rtol = tolerance` and its pressure sub-solve *also* at `rtol = tolerance` — no +margin whatsoever, the inner solve asked to be no better than the answer it feeds. +That is the invariant broken outright rather than a tuning disagreement. + +Loosening the velocity sub-solve likewise does not degrade gracefully — it walks back toward the `preonly` failure above: | velocity sub-KSP rtol | outer its | @@ -265,11 +298,16 @@ Current defaults, matched across the native and rotated paths: | pressure | `0.1 × tolerance` | 200 | The rotated path previously used `0.1 × tolerance` (velocity) and `1.0 × tolerance` -(pressure) — 3× and 10× looser than native, with no recorded reason. Adopting the -native values costs ~17% wall clock with identical outer iteration counts (measured -on both velocity-block routes, isotropic and TI) and reduced a transversely -isotropic fault smoke test from 24 to 15 nonlinear iterations. Cheaper is available -to anyone who measures their own configuration; it is not the default. +(pressure). Adopting the native values costs ~17% wall clock with identical outer +iteration counts (measured on both velocity-block routes, isotropic and TI) and +reduced a transversely isotropic fault smoke test from 24 to 15 nonlinear +iterations. Cheaper is available to anyone who measures their own configuration; it +is not the default. + +The `0.033` and `0.1` factors themselves are inherited from the Citcom +configuration and have no derivation recorded here beyond "well below the outer +tolerance". They are a margin whose *existence* is principled and whose *size* is +convention. ### Detecting a degraded sub-solve diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 9edc867f..f755374d 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -965,11 +965,19 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal "pc_fieldsplit_schur_fact_type": "full", # Pressure sub-solve at native-path parity (pyx Stokes defaults): # FGMRES at 0.1 x tol, GASM on the 1/mu mass (jacobi if only selfp - # exists). This path used to ask for tol, i.e. 10x looser than native. - # Defaults here err toward the more conservative of the two: matching - # native costs ~17% wall clock with identical outer iteration counts - # (measured, both velocity-block routes, isotropic and TI), and a - # too-loose inner solve fails by silent stagnation rather than loudly. + # exists). + # + # The 0.1 is the MARGIN, and it is the point. The Citcom design this + # configuration descends from (Moresi & Solomatov 1995) makes the inner + # solves deliberately inexact, which is why the outer/Schur Krylov must + # be flexible (fgmres, above) — inexact inner solves perturb the search + # directions. But inexact is bounded: an inner solve must still converge + # WELL BELOW the tolerance demanded of the outer solve. This path used to + # ask for `tol` while the outer KSP also asks for `tol` — no margin at + # all, the inner solve permitted to be no better than the answer it + # feeds. Restoring the margin costs ~17% wall clock with identical outer + # iteration counts (measured, both velocity-block routes, isotropic and + # TI); a too-loose inner solve fails by silent stagnation, not loudly. "fieldsplit_pres_ksp_type": "fgmres", "fieldsplit_pres_ksp_rtol": str(tol * 0.1), "fieldsplit_pres_ksp_max_it": "200", From 99761e022ab24b9e3ad43b8b0adae24302a480c5 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 17:03:44 +1000 Subject: [PATCH 12/17] Give the geometric-MG option bundle one owner; pick up mesh-owned hierarchies on the rotated path (#468, #467) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three routes reach a multigrid Stokes velocity block — native (PETSc interpolation between refined DMPlex levels), custom-P on the standard solve path, and custom-P through the rotated free-slip path. They are the same preconditioner reached three ways, not alternatives: custom-P is mandatory wherever native cannot go, namely rotated BCs (the DM-coupled hierarchy cannot express a per-node rotation) and adapt() children (no DMPlex refinement relation). The option bundle was written in two places and had drifted (#468). The native path had been moved to a measured gmres+sor smoother; custom_mg._configure_pcmg never received that update, and never SET mg_levels_ksp_max_it at all, so the smoother inherited whatever last wrote that options prefix — 3 left behind by the GAMG bundle on the standard path, PETSc's own PCMG default of 2 on the rotated path. The same function smoothed differently depending on what had run before it. utilities/multigrid_options.py now owns both bundles (geometric MG and the GAMG fallback) and all three writers read from it. Each bundle DERIVES the stale keys it must clear — every key a sibling bundle sets that it does not — instead of carrying a hand-maintained delete list; the derivation reproduces all three previous hand-written lists exactly. The rotated route's svd coarse solve is a named variant of the shared bundle (coarse="svd") rather than a call-site override, because the Galerkin-coarsened rotated block inherits the rigid-rotation null space where redundant/LU hits a zero pivot (#306). The #276 single-field injection fragility deliberately stays out of the bundle: it is a routing decision, not an option value. Separately (#467), mesh.adapt() leaves a coarse tail on its refinement child so every solver on an adapted mesh gets geometric MG with no per-solver call. The rotated path never consulted it — the standard path's injection hook runs after the rotated dispatch has already returned — so an adapt child under rotated free-slip fell back to GAMG, indistinguishable from having no hierarchy at all. That is the adapt-on-top-faults workflow's own configuration. custom_mg.build_transfers is now the shared "which hierarchy does this solver get?" rule for both paths, with the same opportunistic barycentric->RBF fallback and degrade-to-GAMG on failure. Measured (same operator, RHS and coarse solve; two-level hierarchy, i.e. where the native measurement says the gmres margin is SMALLEST): route pre-fix unified custom-P standard (isotropic annulus) 5 vel its 4 vel its custom-P rotated (TI, eta1/eta0 = 1e-3) 11 vel its 5 vel its 0.683 s 0.392 s (1.74x, linear solve timed in isolation) Also here: the rotated FMG branch now drops its option keys after setup, as the KSP's own keys already were. The per-solve prefix is unique, so those keys previously accumulated in the global database once per timestep. A new test covers the failure mode that introduces — a later setFromOptions finding pc_type gone and abandoning the multigrid — across eight Newton increments. Tests: tests/test_1021_mg_option_bundle.py reads the smoother configuration back off the LIVE PCs for all three routes and asserts they agree, with the coarse difference asserted rather than tolerated (an options-database assertion would not have caught the drift, which was a key nobody wrote). Verified to fail when the pre-fix bundle is reimposed. The mesh-owned pickup is covered serially there and at np>1 in tests/parallel/test_1064, where the transfers are built cross-partition. Validation: test_1014/1015/1016/1017/1018/1020, test_0753/0835/0836/0840 pass; full "level_1 and tier_a" sweep 526 passed with 2 pre-existing failures (#470). Parallel np2 and np4: test_1017_custom_mg_parallel 6 passed, test_1064_rotated_freeslip_parallel 9 passed, test_1066_rotated_datum_parallel 2 passed. No goldens moved — GOLDEN_ANNULUS_FMG is a converged-solution quantity, so a better preconditioner reaches the same answer, and the iteration guards are upper bounds a better smoother only helps. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 43 ++++ docs/developer/subsystems/rotated-freeslip.md | 33 ++- docs/developer/subsystems/solvers.md | 79 +++++- .../cython/petsc_generic_snes_solvers.pyx | 61 +---- src/underworld3/utilities/custom_mg.py | 117 ++++++--- .../utilities/multigrid_options.py | 198 +++++++++++++++ src/underworld3/utilities/rotated_bc.py | 85 ++++--- .../test_1064_rotated_freeslip_parallel.py | 38 ++- tests/test_1021_mg_option_bundle.py | 240 ++++++++++++++++++ 9 files changed, 770 insertions(+), 124 deletions(-) create mode 100644 src/underworld3/utilities/multigrid_options.py create mode 100644 tests/test_1021_mg_option_bundle.py diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 2a8f28f2..036d63ce 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,49 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### One Owner for the Geometric-Multigrid Option Bundle (July 2026) + +**The PETSc option bundle that configures a Stokes velocity block's multigrid +now lives in exactly one module, and all three routes that reach that block read +it from there** (#468), **and rotated free-slip now picks up a mesh-owned +multigrid hierarchy instead of silently discarding it** (#467). + +Three routes reach a multigrid velocity block: native (PETSc interpolation +between refined DMPlex levels), custom-P on the standard solve path, and custom-P +through the rotated free-slip path. They are the same preconditioner reached +three ways, not alternatives — custom-P is *mandatory* wherever native cannot go, +namely rotated boundary conditions and `adapt()` children. The bundle was +written in two places and had drifted: the native path had been moved to a +`gmres`+`sor` smoother on a recorded measurement, and the custom-P routes had +not. Worse, the custom-P writer never *set* the smoother iteration count at all, +so it inherited whatever had last written that options prefix — 3 left behind by +the GAMG bundle on the standard path, PETSc's own default of 2 on the rotated +path. The same function smoothed differently depending on what had run before +it. + +Unifying the bundle recovers, on the same operator, right-hand side and coarse +solve: rotated custom-P velocity-block iterations 11 → 5 (0.68 s → 0.39 s of +linear solve, timed in isolation), standard custom-P 5 → 4, on a *two-level* +hierarchy — the depth at which the native measurement says the gmres margin is +smallest. The bundle also now derives which stale keys it must clear rather than +carrying a hand-maintained list, which is what let the iteration count go unset +in the first place. + +Separately, `mesh.adapt()` leaves a coarse tail on its refinement child so that +every solver on an adapted mesh gets geometric multigrid with no per-solver call. +The rotated path never consulted it — the standard path's injection hook runs +after the rotated dispatch has already returned — so an adapt child under rotated +free-slip fell back to algebraic multigrid, indistinguishable from having no +hierarchy at all. That is the `adapt-on-top-faults` workflow's own configuration +(a fault resolved by local refinement, with rotated free-slip chosen because it +composes with transverse isotropy). Both paths now resolve the hierarchy through +one shared rule, with the same opportunistic degrade-to-GAMG behaviour. + +The regression test reads the smoother configuration back off the **live PETSc +objects** for all three routes and asserts they agree. An options-database +assertion would not have caught the original drift, because the drift was +precisely a key nobody wrote. + ### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026) **Rotated strong free-slip now takes a prescribed wall-normal velocity datum diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 14737fe1..7c47271e 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -41,11 +41,42 @@ There is a single driver, `solve_rotated_freeslip`: a manual outer Newton/Picard loop that rotates the residual and tangent every iteration (`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows, and solves each increment with a self-contained fieldsplit-Schur KSP (geometric -FMG on the custom prolongation when a hierarchy is registered, else GAMG; the +FMG on the custom prolongation when the solver has a hierarchy, else GAMG; the native 1/μ pressure mass as the Schur preconditioner). There is no separate linear path and no up-front nonlinearity probe: a linear model converges after its first increment and the loop self-terminates. +### Where the multigrid hierarchy comes from + +Native geometric FMG cannot serve a rotated solve at all — the DM-coupled +hierarchy has no way to express the per-node rotation — so custom-P multigrid is +not an optimisation here, it is the only geometric option. The hierarchy is +resolved by `custom_mg.build_transfers`, which is the **same rule the standard +path uses**: + +1. an explicit `set_custom_fmg` registration on the solver wins; otherwise +2. a **mesh-owned** coarse tail (`mesh._custom_mg_coarse_meshes`, what + `mesh.adapt()` leaves on a refinement child) is picked up opportunistically, + with the barycentric→RBF builder fallback, so a failed build degrades to GAMG + rather than crashing the solve. + +Step 2 used to be unreachable from here. The standard path's injection hook runs +*after* the rotated dispatch has already returned, and the rotated builder +consulted only `solver._custom_mg` — so an `adapt()` child under rotated +free-slip reported `velocity_pc == "GAMG"`, indistinguishable from having no +hierarchy at all, and lost its multigrid silently (#467). That is the +`adapt-on-top-faults` workflow's own configuration: a fault resolved by local +refinement on an adapt child, with rotated free-slip chosen over Nitsche because +it composes with transverse isotropy. + +The option bundle for both the FMG and the GAMG velocity block comes from +`utilities/multigrid_options.py`, shared with the native and standard custom-P +routes, so the rotated block cannot be configured differently by accident. The +single deliberate difference is the coarse solve: `coarse="svd"`, because the +Galerkin-coarsened rotated block inherits the rigid-rotation null space where +`redundant`/LU hits a zero pivot. See +[the solvers subsystem doc](solvers.md) for the bundle itself. + ### The velocity sub-KSP must converge, not just apply The velocity block is preconditioned by multigrid, but the KSP *around* that diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index 32eaff42..3f492df1 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -175,7 +175,7 @@ def validate_unknowns_sharing(multi_material_model): assert model.Unknowns.DFDt is reference_unknowns.DFDt, \ f"Model {i} $D\mathbf{{F}}/Dt$ not shared - stress history will be wrong" ``` -- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" below; multigrid and Schur preconditioner choice is still undocumented +- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" and "The multigrid option bundle, and its one owner" below; Schur preconditioner choice is still undocumented - Could benefit from optimization examples ## Choosing the Krylov method for a fieldsplit sub-solve @@ -327,6 +327,81 @@ application, and the velocity KSP is applied once per Schur `MatMult` — i.e. o pressure Krylov iteration. That number is a sample, not work. `SolverInstrumentation` (`systems/solver_health.py`) is the mechanism that sums properly. +## The multigrid option bundle, and its one owner + +Three routes reach a multigrid velocity block, and they are **the same +preconditioner reached three ways**, not three alternatives: + +| route | when | prolongation | +|---|---|---| +| native | mesh built with `refinement >= 1`, ordinary BCs | PETSc `DMCreateInterpolation` between refined DMPlex levels | +| custom-P, standard path | `set_custom_fmg`, or an `adapt()` child's mesh-owned coarse tail | barycentric / RBF, Galerkin coarse operators | +| custom-P, rotated path | rotated free-slip, via `rotated_bc` | as above, with the fine prolongation rotated, `P̂ = Q_v·P` | + +custom-P is **mandatory** wherever native cannot go: rotated boundary conditions +(the DM-coupled hierarchy cannot express a per-node rotation) and non-nested +grids (`adapt()` children have no DMPlex refinement relation). So the routes are +not ranked — the one that matters most for adapted and curved-boundary work is +custom-P. + +The option *values* live in one module, `utilities/multigrid_options.py`. Every +writer reads a bundle from there and applies it to its own options object under +its own prefix; nobody writes a multigrid option value anywhere else. That is +structural rather than stylistic: when the bundle was written in two places, the +native path was deliberately moved to a measured `gmres`+`sor` smoother and the +custom-P path was not, and the custom-P routes ran `richardson` at an iteration +count **nobody had set** — inherited from whatever last wrote that options +prefix (3 left behind by the GAMG bundle on the standard path, PETSc's own PCMG +default of 2 on the rotated path). The same function smoothed differently +depending on what had run before it. + +### What a bundle carries + +A bundle is the settings it sets *and* the keys it must clear. The clear-list is +derived, not hand-written: it is every key any sibling bundle sets that this one +does not. These bundles share an options prefix, so switching a block from GAMG +to geometric MG leaves the GAMG-only keys behind and `setFromOptions` will +happily re-read them. Deriving the list means a key added to one bundle +automatically becomes stale for the others. + +Two consequences worth stating outright: + +- **Set the smoother iteration count, never inherit it.** A bundle that omits + `mg_levels_ksp_max_it` is not "using the default" — it is using whatever the + last writer left. +- **The measured smoother is `gmres`+`sor` at 4 iterations.** Chebyshev needs + eigenvalue estimates of the smoothed operator, which are fragile on the + indefinite / variable-viscosity velocity block. Richardson is stationary and + degrades on the non-symmetric operator the consistent-Newton tangent produces: + measured on the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested + 4-level hierarchy) the per-cycle contraction is ρ = 0.75 richardson against + 0.56 gmres at the *same* four iterations, and the gmres margin **grows with + depth** (5% at 3 levels, 25% at 4). + +### The one legitimate per-route difference + +The coarse solve. The Galerkin-coarsened **rotated** velocity block inherits the +rigid-rotation null space of the constrained problem (a closed circle: one mode; +a spherical shell: three), and `redundant`/LU hits a zero pivot there — +`SUBPC_ERROR`, outer reason −11. So the rotated route asks for +`geometric_mg_bundle(coarse="svd")`, which is null-space robust and cheap on a +small coarse level. This is a named variant of the shared bundle, not a +call-site override, so it is visible in the same place as everything else. + +The other native/custom-P asymmetry — that native FMG is unusable for +single-field solvers because `DMCreateInjection` cannot reliably be built on a +refined DMPlex (#276) — is deliberately *not* in the bundle. It is a routing +decision (which route a solver may take), not an option value, and lives with +the route choice in `_apply_preconditioner_options`. + +### Testing it + +`tests/test_1021_mg_option_bundle.py` reads the smoother configuration back off +the **live PETSc objects** after setup — not out of the options database — for +all three routes and asserts they agree, with the coarse-solve difference +asserted rather than tolerated. Options-database assertions would not have caught +the original drift, because the drift was precisely a key that was never written. + ## Critical Stability Note ```{warning} Solver Stability is Paramount @@ -337,7 +412,7 @@ pressure Krylov iteration. That number is a sample, not work. `SolverInstrumenta ```{note} For Contributors This well-documented subsystem could benefit from: -- Preconditioner selection guidance (Krylov sub-solve choice is now covered; multigrid and Schur preconditioner selection is not) +- Preconditioner selection guidance (Krylov sub-solve choice and the multigrid option bundle are now covered; Schur preconditioner selection is not) - Performance tuning documentation - Convergence analysis examples - Scaling studies and optimization diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 5ff3ce6b..21c3520b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -12,6 +12,7 @@ from underworld3.utilities._jitextension import getext, JITCallbackSet import underworld3.timing as timing from underworld3.utilities._api_tools import uw_object +from underworld3.utilities import multigrid_options from underworld3.function import expression as public_expression expression = lambda *x, **X: public_expression(*x, _unique_name_generation=True, **X) @@ -618,46 +619,14 @@ class SolverBaseClass(uw_object): ) want_fmg = False + # The option VALUES live in utilities.multigrid_options, which is the + # single owner shared with the custom-P routes (custom_mg, rotated_bc) — + # the routes are the same preconditioner reached three ways and must not + # be configured from three places (#468). Coarse solve: the native + # hierarchy is not rotated, so its coarse operator carries no inherited + # null space and redundant+LU is right here. if want_fmg: - # Geometric Full Multigrid on the refinement hierarchy. Galerkin - # (RAP) coarse operators are required because UW3 does not install - # residual/Jacobian callbacks on the coarse DMs. - opts[f"{prefix}pc_type"] = "mg" - opts[f"{prefix}pc_mg_type"] = "full" # FMG (F-cycle) - opts[f"{prefix}pc_mg_galerkin"] = "both" # RAP coarse operators - # gmres+sor, sized for a DEEP hierarchy (the only kind worth having: - # a two-level cycle is a coarse-grid correction, not a V-cycle, and is - # not worth special-casing). Chebyshev needs eigenvalue estimates of the - # smoothed operator, which are fragile on the indefinite / - # variable-viscosity velocity block and diverge. Richardson is - # stationary and degrades on the NON-SYMMETRIC operator produced by the - # consistent-Newton tangent. Measured on the Spiegelman notch (Drucker- - # Prager, eta contrast 1e26) over a nested 4-level hierarchy: contraction - # per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the SAME four - # smoother iterations -- and the gmres margin GROWS with depth (5% at 3 - # levels, 25% at 4), because deeper cycles apply the smoother on more - # coarse operators. Four iterations, not more: per unit work gmres/4 - # (rho^(1/4) = 0.87) beats gmres/8 (0.91). - opts[f"{prefix}mg_levels_ksp_type"] = "gmres" - opts[f"{prefix}mg_levels_pc_type"] = "sor" - opts[f"{prefix}mg_levels_ksp_max_it"] = 4 - # Run EXACTLY max_it smoother iterations: no residual-norm computation - # and no convergence test, so every V-cycle costs the same. A Krylov - # smoother makes the cycle non-stationary, which is why the velocity - # block is fgmres (flexible) rather than gmres -- see the fieldsplit - # defaults in the Stokes __init__. - opts[f"{prefix}mg_levels_ksp_norm_type"] = "none" - opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None - # redundant+lu, not bare lu: a bare serial LU cannot factor a - # distributed coarse matrix and fails at np>1 (DIVERGED_LINEAR_SOLVE - # after 0 iterations). redundant gathers the (small) coarse system to - # one rank and is identical to lu in serial — so it is np-safe by - # default without surprising small-np users. - opts[f"{prefix}mg_coarse_pc_type"] = "redundant" - opts[f"{prefix}mg_coarse_redundant_pc_type"] = "lu" - # Clear stale GAMG-only keys so toggling back and forth is clean. - for key in ("pc_gamg_type", "pc_gamg_repartition", "pc_gamg_agg_nsmooths"): - opts.delValue(f"{prefix}{key}") + multigrid_options.geometric_mg_bundle().apply(opts, prefix) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -668,19 +637,7 @@ class SolverBaseClass(uw_object): f"mesh with refinement >= 1 to enable geometric multigrid.", stacklevel=2, ) - opts[f"{prefix}pc_type"] = "gamg" - opts[f"{prefix}pc_gamg_type"] = "agg" - opts[f"{prefix}pc_gamg_repartition"] = True - opts[f"{prefix}pc_mg_type"] = "additive" - opts[f"{prefix}pc_gamg_agg_nsmooths"] = 2 - opts[f"{prefix}mg_levels_ksp_max_it"] = 3 - opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None - # Clear stale geometric-MG-only keys. - for key in ("pc_mg_galerkin", "mg_levels_ksp_type", - "mg_levels_pc_type", "mg_levels_ksp_norm_type", - "mg_coarse_pc_type", - "mg_coarse_redundant_pc_type"): - opts.delValue(f"{prefix}{key}") + multigrid_options.gamg_bundle().apply(opts, prefix) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index d4de0371..5d96d272 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -43,6 +43,8 @@ import numpy as np from petsc4py import PETSc +from underworld3.utilities import multigrid_options + __all__ = ["barycentric_prolongation", "rbf_prolongation", "inject_custom_mg", "CustomMGHierarchy", "set_custom_fmg", "sbr_refine", "sbr_refine_where", "nvb_refine"] @@ -593,30 +595,29 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps): +def _configure_pcmg(pc, Ps, coarse="redundant"): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. - Writes the MG bundle into the options DB under the PC's OWN options prefix + The option VALUES come from :func:`multigrid_options.geometric_mg_bundle` — + the same bundle the native (DMPlex-refinement) route applies — so custom-P + and native multigrid cannot be configured differently (#468). ``coarse`` + selects the coarse-solve variant: ``"svd"`` on the rotated path, whose + Galerkin-coarsened operator inherits the rigid-rotation null space. + + The bundle is written into the options DB under the PC's OWN options prefix (``pc.getOptionsPrefix()`` — e.g. ``Solver_N_`` for the scalar top-level PC, - ``Solver_N_fieldsplit_velocity_`` for the Stokes velocity sub-PC) and removes - any gamg keys BEFORE ``setFromOptions`` — otherwise ``setFromOptions`` re-reads - a lingering ``pc_type=gamg`` and reverts ``setType("mg")``. ``setMGInterpolation`` - persists through ``setFromOptions``; the first ``PCSetUp`` builds the coarse - operators from our P (no ``MatProductReplaceMats`` shape bug, since the PCMG is - fresh and P's size is fixed).""" + ``Solver_N_fieldsplit_velocity_`` for the Stokes velocity sub-PC) BEFORE + ``setFromOptions``, and it clears the gamg keys — otherwise ``setFromOptions`` + re-reads a lingering ``pc_type=gamg`` and reverts ``setType("mg")``. + ``setMGInterpolation`` persists through ``setFromOptions``; the first + ``PCSetUp`` builds the coarse operators from our P (no + ``MatProductReplaceMats`` shape bug, since the PCMG is fresh and P's size is + fixed).""" nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" - opts = PETSc.Options() - opts.setValue(prefix + "pc_type", "mg") - opts.setValue(prefix + "pc_mg_type", "full") - opts.setValue(prefix + "pc_mg_galerkin", "both") - opts.setValue(prefix + "mg_levels_ksp_type", "richardson") - opts.setValue(prefix + "mg_levels_pc_type", "sor") - opts.setValue(prefix + "mg_coarse_pc_type", "redundant") - opts.setValue(prefix + "mg_coarse_redundant_pc_type", "lu") - for key in ("pc_gamg_type", "pc_gamg_repartition", "pc_gamg_agg_nsmooths"): - opts.delValue(prefix + key) + multigrid_options.geometric_mg_bundle(coarse=coarse).apply( + PETSc.Options(), prefix) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -963,25 +964,51 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", solver.is_setup = False -def auto_inject_custom_mg(solver, field_id=None): - """Solve-hook entry: inject custom-P FMG from either a solver-set hierarchy - (``set_custom_fmg``) or a **mesh-owned** one (``mesh.adapt`` refinement child). +def build_transfers(solver, field_id=None): + """The custom-P prolongations this solver should drive, built and ready to + install — from either a solver-set hierarchy (``set_custom_fmg``) or a + **mesh-owned** one (a ``mesh.adapt`` refinement child). + + This is the shared resolution rule for every route that can drive custom-P + multigrid: the standard solve path via :func:`auto_inject_custom_mg`, and the + rotated free-slip path, which builds its own IS-based fieldsplit inside + ``utilities.rotated_bc`` and so never reaches the standard injection hook + (#467). Both must answer "which hierarchy does this solver get?" the same way. A refinement child carries ``mesh._custom_mg_coarse_meshes`` (the static - coarse tail). The first time a solver on such a mesh solves, we lazily build a - :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` targeting ``field_id`` - (0 for the Stokes velocity block, None for scalar/vector) and register it on - the solver — so every solver on an adapted mesh drives geometric MG with no - per-solver call. A solver-set hierarchy (if present) always wins. + coarse tail), so a :class:`CustomMGHierarchy` ``[*coarse, solver.mesh]`` + targeting ``field_id`` (0 for the Stokes velocity block, None for + scalar/vector) is built lazily on first solve — every solver on an adapted + mesh drives geometric MG with no per-solver call. A solver-set hierarchy (if + present) always wins. + + Parameters + ---------- + solver : SolverBaseClass + Built solver (``_build`` has run) whose hierarchy is wanted. + field_id : int or None + Field index for multi-field solvers (0 = Stokes velocity), None = single + field. + + Returns + ------- + (CustomMGHierarchy, list of Mat) or (None, None) + ``(None, None)`` means "no hierarchy here" — either none is registered or + an OPPORTUNISTIC mesh-owned build failed and the caller should keep its + default preconditioner. A build failure on a solver-set hierarchy raises + instead: the user asked for it explicitly. """ - # Solver-set hierarchy (set_custom_fmg): the user asked for it explicitly — - # build + install directly and let any error surface. - if solver._custom_mg is not None: - inject_custom_mg(solver) - return + cfg = getattr(solver, "_custom_mg", None) + if cfg is not None: + if not (isinstance(cfg, dict) and cfg.get("mode") == "hierarchy"): + raise NotImplementedError( + "the legacy set_custom_mg registration has no hierarchy to " + "resolve; use set_custom_fmg().") + h = cfg["hierarchy"] + return h, h.build(solver) # explicit request: errors surface # Mesh-owned hierarchy (adapt() child): OPPORTUNISTIC auto-pickup. It must never - # crash a solve, so build the transfers (which now validate the finest reduced + # crash a solve, so build the transfers (which validate the finest reduced # map against the assembled operator — see CustomMGHierarchy.build) inside a # try/except and fall back to the solver's default preconditioner on any failure. # The finest map is derived from the DM section AFTER snes.setUp() finalizes it, @@ -989,7 +1016,7 @@ def auto_inject_custom_mg(solver, field_id=None): # semi-Lagrangian advection-diffusion (which earlier had to be skipped). coarse = getattr(solver.mesh, "_custom_mg_coarse_meshes", None) if coarse is None: - return # nothing to inject + return None, None # nothing to inject builder = getattr(solver.mesh, "_custom_mg_builder", "barycentric") # Retry with the RBF builder before abandoning geometric MG. The @@ -1029,7 +1056,29 @@ def auto_inject_custom_mg(solver, field_id=None): warnings.warn( f"custom_mg: mesh-owned FMG build failed ({exc}); using the " "solver's default preconditioner.") - return + return None, None + + return h, Ps + + +def auto_inject_custom_mg(solver, field_id=None): + """Solve-hook entry on the STANDARD path: resolve this solver's custom-P + hierarchy (:func:`build_transfers`) and install it on the managed PC. + + The rotated free-slip path does not come through here — it builds its own + IS-based fieldsplit and calls :func:`build_transfers` directly. + """ + # Solver-set hierarchy (set_custom_fmg, or the deprecated set_custom_mg): + # the user asked for it explicitly — build + install directly and let any + # error surface. inject_custom_mg is the one place that still understands the + # legacy registration. + if solver._custom_mg is not None: + inject_custom_mg(solver) + return + + h, Ps = build_transfers(solver, field_id=field_id) + if h is None: + return # Dimensional guard (checkable for the monolithic operator, field_id is None): # the finest transfer must chain to the operator PCMG will Galerkin against. diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py new file mode 100644 index 00000000..2245a371 --- /dev/null +++ b/src/underworld3/utilities/multigrid_options.py @@ -0,0 +1,198 @@ +r""" +The multigrid option bundles for a solver's managed preconditioner block. + +Three routes reach a multigrid velocity block in UW3, and all three must be +configured the same way: + +============================ ============================================== +route prolongation +============================ ============================================== +native PETSc ``DMCreateInterpolation`` between refined + DMPlex levels (mesh built with ``refinement>=1``) +custom-P, standard path barycentric / RBF, Galerkin coarse operators + (:mod:`underworld3.utilities.custom_mg`) +custom-P, rotated path as above, with the FINE prolongation rotated + (:mod:`underworld3.utilities.rotated_bc`) +============================ ============================================== + +custom-P is mandatory wherever native cannot go — rotated boundary conditions +(the DM-coupled hierarchy cannot express the per-node rotation) and non-nested +grids (``adapt()`` children have no DMPlex refinement relation) — so the routes +are not alternatives, they are the same preconditioner reached three ways. + +This module is where the bundles live so that they cannot drift apart. Each +writer reads a bundle from here and applies it to its own options object under +its own prefix; nobody writes a multigrid option value anywhere else. + +A bundle carries two things: the settings it *sets*, and the **stale** keys it +must *remove*. The stale list is derived, not hand-written: it is every key any +other bundle sets that this one does not. That matters because these bundles are +written into a shared options database under a shared prefix — toggling a block +from GAMG to geometric MG leaves the GAMG-only keys behind, and ``setFromOptions`` +will happily re-read them. Hand-maintained delete lists are exactly what let the +smoother iteration count go unset and be inherited from whatever ran before +(issue #468). + +Examples +-------- +Apply the geometric bundle to a solver's managed block:: + + from underworld3.utilities import multigrid_options + multigrid_options.geometric_mg_bundle().apply(self.petsc_options, prefix) + +Read the settings without writing them (the rotated path stages its options in a +dict so it can drop them again after ``setUp``):: + + cfg.update({f"fieldsplit_vel_{k}": v + for k, v in multigrid_options.gamg_bundle().settings.items()}) +""" + +from typing import NamedTuple + +__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", + "GEOMETRIC_MG_COARSE_SOLVERS"] + + +#: Coarse-solve variants of the geometric bundle. ``"redundant"`` (redundant+LU) +#: is the default and is np-safe; ``"svd"`` is required whenever the coarse +#: operator inherits a null space — see :func:`geometric_mg_bundle`. +GEOMETRIC_MG_COARSE_SOLVERS = ("redundant", "svd") + + +class MGBundle(NamedTuple): + """A preconditioner option bundle: the keys to set, and the keys to clear. + + ``settings`` maps an option suffix (no prefix, no leading ``-``) to its + value; ``None`` means a bare PETSc flag. ``stale`` lists suffixes this + bundle does not own but a sibling bundle does, which must be removed before + ``setFromOptions`` so a previously-applied bundle cannot leak through. + """ + + settings: dict + stale: tuple + + def apply(self, opts, prefix=""): + """Write this bundle into ``opts`` (a ``PETSc.Options``) under ``prefix``. + + ``prefix`` is the managed block's own prefix — ``""`` for a top-level PC, + ``"fieldsplit_velocity_"`` for the Stokes velocity sub-block — and is + applied on top of whatever prefix ``opts`` itself carries. + """ + for key, value in self.settings.items(): + opts.setValue(prefix + key, value) + for key in self.stale: + opts.delValue(prefix + key) + + +def _geometric_mg_settings(coarse): + """The geometric-MG settings for one coarse-solve variant.""" + settings = { + "pc_type": "mg", + "pc_mg_type": "full", # FMG (F-cycle) + # Galerkin (RAP) coarse operators are REQUIRED: UW3 installs no + # residual/Jacobian callbacks on the coarse DMs, so PETSc cannot + # re-discretise the operator there. + "pc_mg_galerkin": "both", + # gmres+sor, sized for a DEEP hierarchy (the only kind worth having: a + # two-level cycle is a coarse-grid correction, not a V-cycle, and is not + # worth special-casing). Chebyshev needs eigenvalue estimates of the + # smoothed operator, which are fragile on the indefinite / + # variable-viscosity velocity block and diverge. Richardson is stationary + # and degrades on the NON-SYMMETRIC operator produced by the + # consistent-Newton tangent. Measured on the Spiegelman notch + # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy: + # contraction per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the + # SAME four smoother iterations -- and the gmres margin GROWS with depth + # (5% at 3 levels, 25% at 4), because deeper cycles apply the smoother on + # more coarse operators. Four iterations, not more: per unit work gmres/4 + # (rho^(1/4) = 0.87) beats gmres/8 (0.91). + "mg_levels_ksp_type": "gmres", + "mg_levels_pc_type": "sor", + # SET the count, never inherit it. PCMG's own default is 2 and the GAMG + # bundle below leaves 3 under the same prefix, so a bundle that omits + # this key smooths differently depending on what ran before it (#468). + "mg_levels_ksp_max_it": 4, + # Run EXACTLY max_it smoother iterations: no residual-norm computation and + # no convergence test, so every V-cycle costs the same. A Krylov smoother + # makes the cycle non-stationary, which is why the velocity block is + # fgmres (flexible) rather than gmres. + "mg_levels_ksp_norm_type": "none", + "mg_levels_ksp_converged_maxits": None, + } + if coarse == "redundant": + # redundant+lu, not bare lu: a bare serial LU cannot factor a distributed + # coarse matrix and fails at np>1 (DIVERGED_LINEAR_SOLVE after 0 + # iterations). redundant gathers the (small) coarse system to one rank and + # is identical to lu in serial — np-safe by default without surprising + # small-np users. + settings["mg_coarse_pc_type"] = "redundant" + settings["mg_coarse_redundant_pc_type"] = "lu" + elif coarse == "svd": + # The Galerkin-coarsened ROTATED velocity block inherits every + # rigid-rotation null-space mode of the constrained problem (a closed + # circle: one; a spherical shell: three), and redundant/LU hits a zero + # pivot there (SUBPC_ERROR, outer reason -11 — the #306 fix). SVD is + # null-space robust and the coarse level is small. Same choice as the + # native spherical FMG setups. + settings["mg_coarse_pc_type"] = "svd" + else: + raise ValueError( + f"coarse must be one of {GEOMETRIC_MG_COARSE_SOLVERS} (got {coarse!r})") + return settings + + +def _gamg_settings(): + """The algebraic-multigrid (GAMG) settings — the fallback whenever no + geometric hierarchy is available.""" + return { + "pc_type": "gamg", + "pc_gamg_type": "agg", + "pc_gamg_repartition": True, + "pc_mg_type": "additive", + "pc_gamg_agg_nsmooths": 2, + "mg_levels_ksp_max_it": 3, + "mg_levels_ksp_converged_maxits": None, + } + + +def _all_keys(): + """Every option key any bundle here owns — the basis for the stale lists.""" + keys = set(_gamg_settings()) + for coarse in GEOMETRIC_MG_COARSE_SOLVERS: + keys |= set(_geometric_mg_settings(coarse)) + return keys + + +def _bundle(settings): + return MGBundle(settings=settings, + stale=tuple(sorted(_all_keys() - set(settings)))) + + +def geometric_mg_bundle(coarse="redundant"): + """Geometric multigrid (FMG F-cycle) on an explicit hierarchy. + + Parameters + ---------- + coarse : {"redundant", "svd"} + The coarse-level solve. ``"redundant"`` (redundant+LU) is the default. + Use ``"svd"`` when the coarse operator inherits a null space — which the + Galerkin-coarsened *rotated* velocity block always does, because the + rigid rotations survive the constraint. + + Returns + ------- + MGBundle + Settings and stale keys; see :meth:`MGBundle.apply`. + """ + return _bundle(_geometric_mg_settings(coarse)) + + +def gamg_bundle(): + """Algebraic multigrid — the fallback when no geometric hierarchy exists. + + Returns + ------- + MGBundle + Settings and stale keys; see :meth:`MGBundle.apply`. + """ + return _bundle(_gamg_settings()) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index f755374d..b83bad4d 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -9,9 +9,12 @@ no separate linear path and no up-front nonlinearity probe. Each increment is solved by a self-contained fieldsplit-Schur KSP by default: the -velocity block is geometric FMG on the custom prolongation (``set_custom_fmg``) when a -hierarchy is registered (the PREFERRED route), else GAMG tuned to the native path's -settings; the Schur complement is preconditioned by the native 1/mu pressure mass +velocity block is geometric FMG on the custom prolongation whenever this solver has a +hierarchy — registered by ``set_custom_fmg`` or owned by an ``adapt()`` child mesh — +which is the PREFERRED route, else GAMG; the option bundle for both is the shared one +in ``underworld3.utilities.multigrid_options``, so the rotated velocity block is +configured identically to the native and standard custom-P routes. +The Schur complement is preconditioned by the native 1/mu pressure mass (the Pmat p-p block, exactly the standard path's ``schur_precondition=a11``), with a constant-pressure nullspace on the inner Schur solve for enclosed domains; direct MUMPS LU is opt-in via ``solver._rotated_use_lu``. @@ -32,6 +35,11 @@ from underworld3.utilities.boundary_flux import ( _boundary_stratum_is, _desmear, write_boundary_scalar_field) +# The multigrid option bundles are owned in ONE place and shared with the native +# and standard custom-P routes, so the rotated velocity block cannot be +# configured differently from the others (#468). +from underworld3.utilities import multigrid_options + # Monotonic counter so each rotated solve gets a UNIQUE PETSc options prefix. With a # fixed prefix, sequential rotated solves (e.g. two solvers in one script, or one # solver in a time-stepping loop) share and re-set the same global-options keys; the @@ -847,19 +855,33 @@ def rotated_residual(uvec, keep_cartesian=False): def _build_rotated_custom_Pl(solver, Q, normal_rows): """The rotated custom-FMG prolongation list [*coarse, Q_v·P_fine] for the - velocity block, or None if no hierarchy is registered. Depends only on Q and + velocity block, or None if this solver has no hierarchy. Depends only on Q and the mesh (NOT the solution), so the rotated Newton loop builds it ONCE and reuses - it across Newton iterations (the prolongation build is the expensive part).""" - if getattr(solver, "_custom_mg", None) is None: + it across Newton iterations (the prolongation build is the expensive part). + + The hierarchy is resolved by ``custom_mg.build_transfers``, which is the same + rule the standard path uses: an explicit ``set_custom_fmg`` registration wins, + otherwise a MESH-OWNED coarse tail (what ``mesh.adapt()`` leaves on a + refinement child) is picked up opportunistically. That mesh-owned pickup used + to be unreachable from here — the standard path's injection hook runs after + the rotated dispatch has already returned — so an adapt child under rotated + free-slip silently lost its multigrid and solved on GAMG (#467). This is the + ``adapt-on-top-faults`` workflow's own configuration.""" + from underworld3.utilities import custom_mg + h, Ps = custom_mg.build_transfers(solver, field_id=0) + if h is None: return None vel_is = solver._subdict["velocity"][0] vis = np.asarray(vel_is.getIndices()) g2blk = {int(g): k for k, g in enumerate(vis)} Qv = Q.createSubMatrix(vel_is, vel_is) nrows_blk = sorted({g2blk[g] for g in normal_rows if g in g2blk}) - Ps = solver._custom_mg["hierarchy"].build(solver) Pfine = Qv.matMult(Ps[-1]) Pfine.zeroRows(nrows_blk, diag=0.0) + # Remember an opportunistic mesh-owned pickup, so a later solve on this solver + # (rotated or not) reuses the hierarchy instead of re-resolving it. + if getattr(solver, "_custom_mg", None) is None: + solver._custom_mg = {"mode": "hierarchy", "hierarchy": h, "verbose": False} return list(Ps[:-1]) + [Pfine] @@ -909,8 +931,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal custom_Pl=None, nsp=None, Mp=None, ctx=None): """Solve the rotated saddle with a SELF-CONTAINED fieldsplit-Schur KSP on the rotated operator. The velocity block is geometric FMG on the CUSTOM prolongation - (PR#290, rotated) when a hierarchy is registered (``set_custom_fmg``), else GAMG - (tuned to the native path's settings). + (PR#290, rotated) whenever the solver has a hierarchy — ``set_custom_fmg`` or a + mesh-owned ``adapt()`` tail — else GAMG. Both bundles come from + ``utilities.multigrid_options``, shared with the native and standard routes. A plain rotated Mat has no DM field info, so UW3's DM-coupled fieldsplit cannot split it — we build the split from EXPLICIT velocity/pressure index sets. For the @@ -988,21 +1011,18 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal cfg["pc_fieldsplit_schur_precondition"] = "selfp" cfg["fieldsplit_pres_pc_type"] = "jacobi" if custom_Pl is None: - # GAMG fallback velocity block, tuned to native parity (pyx Stokes - # defaults). NOTE: the custom-FMG route is the preferred velocity - # block — this applies only when no hierarchy is registered. + # GAMG fallback velocity block. The bundle is the SHARED one (see + # utilities.multigrid_options), so it cannot drift from what the + # standard path applies. NOTE: the custom-FMG route is the preferred + # velocity block — this applies only when no hierarchy is available. + # No stale keys to clear: `pfx` is unique per rotated solve. cfg.update({ "fieldsplit_vel_ksp_type": "fgmres", "fieldsplit_vel_ksp_rtol": str(tol * 0.033), "fieldsplit_vel_ksp_max_it": "200", - "fieldsplit_vel_pc_type": "gamg", - "fieldsplit_vel_pc_gamg_type": "agg", - "fieldsplit_vel_pc_gamg_repartition": "true", - "fieldsplit_vel_pc_mg_type": "additive", - "fieldsplit_vel_pc_gamg_agg_nsmooths": "2", - "fieldsplit_vel_mg_levels_ksp_max_it": "3", - "fieldsplit_vel_mg_levels_ksp_converged_maxits": "true", }) + cfg.update({f"fieldsplit_vel_{key}": value for key, value + in multigrid_options.gamg_bundle().settings.items()}) else: # Custom-FMG velocity block, wrapped in a short FGMRES — NOT `preonly`. # PCFieldSplit applies the Schur complement S = A11 - A10 A00^-1 A01 @@ -1054,20 +1074,23 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal A_vv, P_vv = vel_pc.getOperators() vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) - custom_mg._configure_pcmg(vel_pc, custom_Pl) - # The Galerkin-coarsened ROTATED velocity block inherits every - # rigid-rotation nullspace mode of the constrained problem (a - # closed circle: one; a spherical shell: three) — the default - # redundant/LU coarse solve hits a zero pivot (SUBPC_ERROR, - # outer reason -11). SVD is nullspace-robust and the coarse - # level is small; same choice as the native spherical FMG setups. + # coarse="svd": the Galerkin-coarsened ROTATED velocity block + # inherits every rigid-rotation nullspace mode of the constrained + # problem (a closed circle: one; a spherical shell: three), and + # the default redundant/LU coarse solve hits a zero pivot there + # (SUBPC_ERROR, outer reason -11). Everything else in the bundle + # is identical to the native and standard custom-P routes. + custom_mg._configure_pcmg(vel_pc, custom_Pl, coarse="svd") + vel_pc.setUp() + # The bundle went into the GLOBAL options DB under the velocity + # sub-PC's prefix, which is derived from `pfx` and so is unique to + # this solve — drop it again now it is consumed, as the `finally` + # below does for the KSP's own keys. vopts = PETSc.Options() vpfx = vel_pc.getOptionsPrefix() or "" - vopts.setValue(vpfx + "mg_coarse_pc_type", "svd") - vopts.delValue(vpfx + "mg_coarse_redundant_pc_type") - vel_pc.setFromOptions() - vel_pc.setUp() - vopts.delValue(vpfx + "mg_coarse_pc_type") + for key in multigrid_options.geometric_mg_bundle( + coarse="svd").settings: + vopts.delValue(vpfx + key) # Constant-pressure nullspace on the Schur COMPLEMENT (enclosed # domains): the IS-built fieldsplit does not propagate the coupled # nullspace to the inner Schur solve, which is otherwise singular diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 61c630ef..1ebe04ed 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -236,16 +236,23 @@ def fit(mask): ) -def _annulus_fmg_diagnostics(): +def _annulus_fmg_diagnostics(mesh_owned=False): """Annulus radial free-slip whose velocity block is CUSTOM GEOMETRIC FMG on a - nested hierarchy (coarse annulus -> refine -> refine), via set_custom_fmg — no - GAMG, no direct solve. Returns (velocity L2, leakage L2 on Lower, on Upper).""" + nested hierarchy (coarse annulus -> refine -> refine) — no GAMG, no direct + solve. Returns (velocity L2, leakage L2 on Lower, on Upper). + + ``mesh_owned`` selects how the hierarchy reaches the solver: ``False`` is an + explicit ``set_custom_fmg`` registration, ``True`` is the MESH-owned coarse + tail an ``adapt()`` child carries. Both must produce the same solve (#467).""" RI, RO = 0.5, 1.0 m0 = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) dm1 = m0.dm.refine() dm2 = dm1.refine() coarse = [m0, _wrap(dm1, m0)] fine = _wrap(dm2, m0) + if mesh_owned: + fine._custom_mg_coarse_meshes = coarse + fine._custom_mg_builder = "barycentric" x, y = fine.X r = sympy.sqrt(x**2 + y**2) @@ -264,10 +271,14 @@ def _annulus_fmg_diagnostics(): s.saddle_preconditioner = 1.0 s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" - custom_mg.set_custom_fmg(s, coarse, builder="barycentric", field_id=0) + if not mesh_owned: + custom_mg.set_custom_fmg(s, coarse, builder="barycentric", field_id=0) s.solve() assert s._rotated_freeslip_info["ksp_reason"] > 0, "custom-FMG rotated solve did not converge" + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG", ( + "rotated velocity block fell back to " + f"{s._rotated_freeslip_info['velocity_pc']} instead of geometric MG") L2 = float(np.sqrt(uw.maths.Integral(fine, v.sym.dot(v.sym)).evaluate())) vr = v.sym[0] * x / r + v.sym[1] * y / r leak_lo = float(np.sqrt(uw.maths.BdIntegral( @@ -424,6 +435,25 @@ def test_rotated_freeslip_annulus_fmg_partition_independent(): f"{leak_up_ref} vs {leak_up}") +def test_rotated_freeslip_mesh_owned_fmg_pickup(): + """#467, in parallel: a coarse tail owned by the MESH — what ``adapt()`` + leaves on a refinement child — must drive geometric MG under rotated + free-slip, on every partition. The rotated path used to consult only an + explicit ``set_custom_fmg`` registration and fell back to GAMG silently. + + The transfers are built cross-partition, so this is not implied by the serial + pickup test (``tests/test_1021_mg_option_bundle.py``); the tail is attached by + hand because what is under test is whether the rotated dispatch consults it, + not how ``adapt()`` produces it.""" + L2, leak_lo, leak_up = _annulus_fmg_diagnostics(mesh_owned=True) + L2_ref, leak_lo_ref, leak_up_ref = GOLDEN_ANNULUS_FMG + assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( + f"mesh-owned FMG annulus velocity L2 differs serial vs np={uw.mpi.size}: " + f"{L2_ref} vs {L2}") + assert np.isclose(leak_lo, leak_lo_ref, rtol=1e-4, atol=0) + assert np.isclose(leak_up, leak_up_ref, rtol=1e-4, atol=0) + + def test_rotated_freeslip_spherical3d_partition_independent(): """3D spherical shell (free-slip inner+outer, all three rotation nullspace modes): the parallel solve reproduces the serial velocity L2, converges, and diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py new file mode 100644 index 00000000..53b8983d --- /dev/null +++ b/tests/test_1021_mg_option_bundle.py @@ -0,0 +1,240 @@ +"""The geometric-multigrid option bundle has ONE owner, and every route reads it. + +Three routes reach a multigrid velocity block — native (DMPlex refinement), +custom-P on the standard solve path, and custom-P through the rotated free-slip +path — and they are the same preconditioner, not alternatives: custom-P is +mandatory wherever native cannot go (rotated BCs, ``adapt()`` children). They +drifted apart because the bundle was written in two places (#468): the custom-P +route ran richardson smoothing at an iteration count nobody set, inherited from +whatever had last written that options prefix (3 left over from the GAMG bundle +on the standard path, PETSc's own default of 2 on the rotated path). + +These tests read the configuration back off the LIVE PETSc objects after setup, +not out of the options database, so they check what actually runs. + +The one legitimate per-route difference is the coarse solve: the +Galerkin-coarsened *rotated* velocity block inherits the rigid-rotation null +space, so redundant/LU hits a zero pivot there and it uses SVD (the #306 fix). +That difference is asserted, not tolerated — if it disappears, the rotated coarse +solve has silently reverted. + +Also here: the rotated path must pick up a MESH-OWNED hierarchy (the coarse tail +an ``adapt()`` child carries), which it previously ignored, silently solving on +GAMG (#467). +""" +import pytest +import sympy + +import underworld3 as uw +from underworld3.utilities import custom_mg, multigrid_options, rotated_bc + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +RES = 0.15 +R_IN, R_OUT = 0.5, 1.0 + + +# --------------------------------------------------------------------------- # +# The owner itself — no solve needed +# --------------------------------------------------------------------------- # +def test_every_bundle_sets_the_smoother_iteration_count(): + """The drift mechanism: a bundle that omits ``mg_levels_ksp_max_it`` inherits + it from whatever last wrote that prefix. Every bundle must SET it.""" + bundles = [multigrid_options.gamg_bundle()] + bundles += [multigrid_options.geometric_mg_bundle(coarse=c) + for c in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS] + for bundle in bundles: + assert "mg_levels_ksp_max_it" in bundle.settings + assert "pc_type" in bundle.settings + + +def test_bundles_clear_each_others_keys(): + """Bundles share an options prefix, so switching between them must leave no + key behind — every key a bundle does not set, a sibling does, and it must + appear in that bundle's stale list.""" + bundles = [multigrid_options.gamg_bundle()] + bundles += [multigrid_options.geometric_mg_bundle(coarse=c) + for c in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS] + owned = set().union(*(set(b.settings) for b in bundles)) + for bundle in bundles: + assert set(bundle.stale) == owned - set(bundle.settings) + + +# --------------------------------------------------------------------------- # +# The three live routes +# --------------------------------------------------------------------------- # +def _stokes(mesh, tag, rotated): + """A buoyancy-driven annulus Stokes solve: fixed inner boundary, and an outer + boundary that is either an ordinary essential BC or rotated free-slip.""" + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + rhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable(f"V{tag}", mesh, mesh.dim, degree=2, + continuous=True) + p = uw.discretisation.MeshVariable(f"P{tag}", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + blob = sympy.exp(-(((x - 0.75) ** 2 + y**2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") + if rotated: + s.add_rotated_freeslip_bc(0.0, "Upper", normal=rhat) + else: + s.add_essential_bc((sympy.oo, 0.0), "Upper") + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-8 + return s + + +def _mg_config(vel_pc): + """Read the multigrid configuration off a live velocity sub-PC. The finest + smoother is representative; level 0 is the coarse solve.""" + assert vel_pc.getType() == "mg" + nlev = vel_pc.getMGLevels() + smoother = vel_pc.getMGSmoother(nlev - 1) + coarse = vel_pc.getMGCoarseSolve() + return { + "levels": nlev, + "mg_type": str(vel_pc.getMGType()), + "smoother_ksp": smoother.getType(), + "smoother_pc": smoother.getPC().getType(), + "smoother_max_it": smoother.getTolerances()[3], + "smoother_norm": str(smoother.getNormType()), + "coarse_ksp": coarse.getType(), + "coarse_pc": coarse.getPC().getType(), + } + + +def _annulus(cell_size): + return uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=cell_size, qdegree=3) + + +def _native_config(): + """Route A: native geometric FMG on a refined DMPlex hierarchy.""" + s = _stokes(uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1), + "nat", rotated=False) + s.solve() + return _mg_config(s.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC()) + + +def _custom_standard_config(): + """Route B: custom-P FMG reached through the standard solve path.""" + s = _stokes(_annulus(RES), "std", rotated=False) + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + s.solve() + return _mg_config(s.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC()) + + +def _custom_rotated_config(monkeypatch): + """Route C: custom-P FMG reached through the rotated free-slip path. + + The rotated KSP is self-contained and is destroyed at the end of the solve, + so the configuration is captured from inside the solve — there is no live PC + left to interrogate afterwards. + """ + s = _stokes(_annulus(RES), "rot", rotated=True) + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + + real_solve = rotated_bc._solve_rotated_iterative + captured = {} + + def capture(solver, Ahat, bhat, Q, Qt, normal_rows, **kw): + result = real_solve(solver, Ahat, bhat, Q, Qt, normal_rows, **kw) + captured.setdefault("config", _mg_config( + result[2]["pc"].getFieldSplitSubKSP()[0].getPC())) + return result + + monkeypatch.setattr(rotated_bc, "_solve_rotated_iterative", capture) + s.solve() + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG" + return captured["config"] + + +def test_all_three_routes_share_one_bundle(monkeypatch): + """Native, standard custom-P and rotated custom-P must smooth identically — + same Krylov smoother, same preconditioner, same iteration count, same cycle. + Only the coarse solve legitimately differs.""" + native = _native_config() + standard = _custom_standard_config() + rotated = _custom_rotated_config(monkeypatch) + + shared = ("mg_type", "smoother_ksp", "smoother_pc", "smoother_max_it", + "smoother_norm", "coarse_ksp") + for key in shared: + assert native[key] == standard[key] == rotated[key], ( + f"route drift on {key!r}: native={native[key]!r} " + f"standard={standard[key]!r} rotated={rotated[key]!r}") + + # The measured bundle, not just self-consistency: gmres/sor at four + # iterations with no norm computation (see multigrid_options for why). + assert native["smoother_ksp"] == "gmres" + assert native["smoother_pc"] == "sor" + assert native["smoother_max_it"] == 4 + + # The ONE deliberate per-route difference (#306): the rotated coarse operator + # inherits the rigid-rotation null space, where redundant/LU hits a zero pivot. + assert native["coarse_pc"] == "redundant" + assert standard["coarse_pc"] == "redundant" + assert rotated["coarse_pc"] == "svd" + + +def test_rotated_picks_up_a_mesh_owned_hierarchy(): + """#467: a coarse tail owned by the MESH — what ``adapt()`` leaves on a + refinement child — must drive geometric MG under rotated free-slip. It used + to be ignored, and the solve fell back to GAMG silently. + + The tail is attached by hand rather than by running ``adapt()``: it is the + same attribute the child carries, and what is under test is whether the + rotated dispatch consults it. + """ + mesh = _annulus(RES) + mesh._custom_mg_coarse_meshes = [_annulus(2 * RES)] + mesh._custom_mg_builder = "barycentric" + + s = _stokes(mesh, "own", rotated=True) + s.solve() + + assert s._rotated_freeslip_info["velocity_pc"] == "custom-FMG" + assert s._rotated_freeslip_info["velocity_pc_type"] == "mg" + + +def test_rotated_fmg_survives_repeated_newton_increments(): + """The rotated path applies its bundle under a per-solve options prefix and + then drops the keys again, so the global database stays bounded under + time-stepping. That is only safe if nothing re-reads them: a later + ``setFromOptions`` on the velocity sub-PC would find ``pc_type`` gone and + silently abandon the multigrid. Several Newton increments, each re-solving + on a refreshed operator, must all keep ``pc=mg``.""" + mesh = _annulus(RES) + s = _stokes(mesh, "nl", rotated=True) + + # power-law viscosity: a genuinely nonlinear solve, several increments + x, y = mesh.X + v = s.Unknowns.u + grad = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + edot = 0.5 * (grad + grad.T) + eII = sympy.sqrt(0.5 * (edot[0, 0] ** 2 + edot[1, 1] ** 2) + + edot[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + s.consistent_jacobian = True + s.tolerance = 1.0e-7 + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + s.solve() + + info = s._rotated_freeslip_info + assert len(info["ksp_its"]) > 1, "not a multi-increment solve — test is vacuous" + assert info["velocity_pc"] == "custom-FMG" + assert info["velocity_pc_type"] == "mg" + + +def test_rotated_without_a_hierarchy_falls_back_to_gamg(): + """The negative control for the test above: no hierarchy anywhere still + reports GAMG, so ``custom-FMG`` above is a real pickup and not a label that + is always set.""" + s = _stokes(_annulus(RES), "none", rotated=True) + s.solve() + assert s._rotated_freeslip_info["velocity_pc"] == "GAMG" From 2550d9a43a603b9d725dbbf07e9094d604b53920 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 29 Jul 2026 19:02:34 +1000 Subject: [PATCH 13/17] Correct the smoother's recorded mechanism: the velocity block stays symmetric under the consistent tangent MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The comment this PR moved into utilities/multigrid_options.py carried a claim we inherited without checking, and which the Layer 3 section of design/nonlinear-solver-homotopy-warmstart.md states outright: that richardson degrades because the consistent-Newton tangent makes the velocity block NON-SYMMETRIC. It does not. For an isotropic eta(eps_II) the Newton term is 2 (eta'/eps_II) (edot(u):edot(du)) (edot(u):edot(v)) a rank-one outer product a (x) a with a = edot(u) — symmetric under swapping du and v. eta depends on grad v only THROUGH its symmetric part, so both factors project onto the same tensor and nothing is left to break the symmetry. Measured as ||A - A^T||_F / ||A||_F on the assembled velocity block (velocity_block_symmetry.py in the rotated_pmat_audit study), against a linear calibration of 5.1e-17, with a control on the size of the Newton contribution so a null result cannot be vacuous: rheology asym Newton term linear (calibration) 5.1e-17 — shear-thinning, isotropic 5.3e-17 2.1e-01 power-law, isotropic 5.3e-17 4.0e-01 Drucker-Prager, pressure-dependent yield 5.1e-17 5.1e-02 transverse isotropic, PICARD tangent 7.2e-02 — transverse isotropic, Newton tangent 8.9e-02 2.0e-01 The consistent tangent contributed 5-40% of ||A|| in every nonlinear case and the block stayed symmetric anyway. Non-symmetry appears only under transverse isotropy — and there it appears under the PICARD tangent too, where a frozen-coefficient form int edot(du):C:edot(v) is symmetric by construction. That is a defect signature, not a property of the tangent: linear TI with anisotropy fully active is symmetric at 6.3e-17. Reported on #457. The Spiegelman-notch contraction measurement (0.75 richardson vs 0.56 gmres at four iterations, margin growing with depth) is unaffected — it is a measurement, not a derivation. The mechanism that fits it AND fits the cost we measured on well-conditioned problems is operator CONDITIONING: at eta contrast 1e26 the SOR-preconditioned spectrum is spread far enough that a stationary iteration stalls where a Krylov smoother adapts its polynomial. Also records the cost, which was previously unmeasured on this route: nested annulus, LINEAR (symmetric) velocity block, eta contrast 1e6, outer KSP timed in isolation, gmres/4 against richardson/3 — iterations x1.45 / x1.60 / x1.20 at 2 / 3 / 4 levels, wall clock x0.86 / x0.77 / x0.55. gmres wins iterations everywhere and loses wall clock everywhere. It stays the default on both routes because the failure it avoids is worse than the cost it carries, and it carries that cost exactly where the problems are easy; the number is recorded so the guardrail can be removed deliberately. No behavioural change — comments, docstrings and docs/developer/subsystems/solvers.md. A TODO(DESIGN) marks the design doc's Layer 3 claim rather than rewriting it. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 73 +++++++++++++++++-- .../utilities/multigrid_options.py | 49 +++++++++++-- 2 files changed, 109 insertions(+), 13 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index 3f492df1..ba2e086b 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -371,12 +371,73 @@ Two consequences worth stating outright: last writer left. - **The measured smoother is `gmres`+`sor` at 4 iterations.** Chebyshev needs eigenvalue estimates of the smoothed operator, which are fragile on the - indefinite / variable-viscosity velocity block. Richardson is stationary and - degrades on the non-symmetric operator the consistent-Newton tangent produces: - measured on the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested - 4-level hierarchy) the per-cycle contraction is ρ = 0.75 richardson against - 0.56 gmres at the *same* four iterations, and the gmres margin **grows with - depth** (5% at 3 levels, 25% at 4). + indefinite / variable-viscosity velocity block. Against richardson, measured on + the Spiegelman notch (Drucker–Prager, η contrast 1e26, nested 4-level + hierarchy): per-cycle contraction ρ = 0.75 richardson against 0.56 gmres at the + *same* four iterations, the margin growing with depth (5% at 3 levels, 25% at 4). + +### Why richardson loses — the recorded reason is wrong + +Both this bundle's comment and the Layer 3 section of +`design/nonlinear-solver-homotopy-warmstart.md` said richardson degrades because +the consistent-Newton tangent makes the velocity block **non-symmetric**. That +mechanism does not survive measurement. + +For an isotropic η(ε̇_II) the Newton term is + +$$2\,\frac{\eta'}{\dot\varepsilon_{II}}\big(\dot\varepsilon(u):\dot\varepsilon(\delta u)\big)\big(\dot\varepsilon(u):\dot\varepsilon(v)\big)$$ + +— a rank-one outer product $a\otimes a$ with $a = \dot\varepsilon(u)$, which is +symmetric. η depends on ∇v only *through* its symmetric part, so both factors +project onto the same tensor and there is nothing left to break the symmetry. +Measured as ‖A−Aᵀ‖_F/‖A‖_F on the assembled velocity block, against a linear +calibration of 5e-17: + +| rheology | ‖A−Aᵀ‖/‖A‖ | Newton term ‖A_N−A_P‖/‖A‖ | +|---|---|---| +| linear (calibration) | 5.1e-17 | — | +| shear-thinning, isotropic | 5.3e-17 | 2.1e-01 | +| power-law, isotropic | 5.3e-17 | 4.0e-01 | +| Drucker–Prager, pressure-dependent yield | 5.1e-17 | 5.1e-02 | +| transverse isotropic, Picard tangent | **7.2e-02** | — | +| transverse isotropic, Newton tangent | **8.9e-02** | 2.0e-01 | + +The third column is the control: if the consistent tangent had contributed +nothing, the symmetry reading would be vacuous. It contributed 5–40% of ‖A‖ in +every case, and the block stayed symmetric anyway. + +Non-symmetry *does* appear under transverse isotropy — but it appears under the +**Picard** tangent too, where a frozen-coefficient form ∫ ε̇(δu):C:ε̇(v) is +symmetric by construction. That is a defect signature rather than a property of +the tangent (see #457, and the linear TI rows, which are symmetric at 6.3e-17 +even with anisotropy fully active). + +The mechanism that does fit everything measured is **operator conditioning**: at +η contrast 1e26 the SOR-preconditioned spectrum is spread far enough that a +stationary iteration stalls where a Krylov smoother adapts its polynomial. It +predicts the notch result, and it predicts the cost below on well-conditioned +problems. + +### The smoother is not free + +On a *well-conditioned* problem gmres/4 is slower than richardson/3 despite +converging in fewer iterations. Nested annulus, linear (symmetric) velocity +block, η contrast 1e6, outer KSP timed in isolation: + +| depth | iterations | wall clock | +|---|---|---| +| 2 levels | ×1.45 | **×0.86** | +| 3 levels | ×1.60 | **×0.77** | +| 4 levels | ×1.20 | **×0.55** | + +The extra sweep costs ~3%; gmres itself costs ~40% per cycle. Custom-P builds its +Galerkin coarse operators from barycentric transfers, denser than the native +nested ones, so the smoother is a larger share of the cycle on that route. + +It stays the default on both routes because the failure it avoids (a stationary +smoother stalling on a high-contrast operator) is worse than the cost it carries, +and it carries that cost exactly where the problems are easy. Removing a guardrail +is the caller's decision; the number is recorded so the decision can be made. ### The one legitimate per-route difference diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 2245a371..cd2e9d27 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -97,15 +97,32 @@ def _geometric_mg_settings(coarse): # two-level cycle is a coarse-grid correction, not a V-cycle, and is not # worth special-casing). Chebyshev needs eigenvalue estimates of the # smoothed operator, which are fragile on the indefinite / - # variable-viscosity velocity block and diverge. Richardson is stationary - # and degrades on the NON-SYMMETRIC operator produced by the - # consistent-Newton tangent. Measured on the Spiegelman notch - # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy: + # variable-viscosity velocity block and diverge. + # + # The measurement, which is what this setting rests on: Spiegelman notch + # (Drucker-Prager, eta contrast 1e26) over a nested 4-level hierarchy, # contraction per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the - # SAME four smoother iterations -- and the gmres margin GROWS with depth - # (5% at 3 levels, 25% at 4), because deeper cycles apply the smoother on - # more coarse operators. Four iterations, not more: per unit work gmres/4 + # SAME four smoother iterations, the margin growing with depth (5% at 3 + # levels, 25% at 4). Four iterations, not more: per unit work gmres/4 # (rho^(1/4) = 0.87) beats gmres/8 (0.91). + # + # TODO(DESIGN): the MECHANISM on record for that result is wrong. Both + # this comment and design/nonlinear-solver-homotopy-warmstart.md (Layer 3) + # said richardson degrades because the consistent-Newton tangent makes the + # velocity block NON-SYMMETRIC. Measured (velocity_block_symmetry.py in the + # rotated_pmat_audit study), that block stays symmetric to machine precision + # -- 5e-17, against a linear calibration of 5e-17 -- under shear-thinning, + # power-law AND Drucker-Prager with pressure-dependent yield, with a control + # confirming the Newton term contributed 5-40% of ||A||. For an isotropic + # eta(eps_II) the Newton term is a rank-one outer product a (x) a with + # a = eps_dot(u), which is symmetric. Non-symmetry does appear under + # transverse isotropy, but it appears under the PICARD tangent there too + # (7e-2), where a frozen-coefficient form is symmetric by construction -- + # i.e. a defect signature, see #457. The likely real mechanism is operator + # CONDITIONING: at eta contrast 1e26 the SOR-preconditioned spectrum is + # spread far enough that a stationary iteration stalls where a Krylov + # smoother adapts its polynomial. That also predicts the cost measured on + # well-conditioned problems (see the note on geometric_mg_bundle). "mg_levels_ksp_type": "gmres", "mg_levels_pc_type": "sor", # SET the count, never inherit it. PCMG's own default is 2 and the GAMG @@ -183,6 +200,24 @@ def geometric_mg_bundle(coarse="redundant"): ------- MGBundle Settings and stale keys; see :meth:`MGBundle.apply`. + + Notes + ----- + The ``gmres``/4 smoother is not free on a well-conditioned problem. Measured on + a nested annulus with a linear (symmetric) velocity block at :math:`\eta` + contrast 1e6, outer KSP timed in isolation, ``gmres``/4 against + ``richardson``/3: it wins on iterations at every depth (×1.45, ×1.60, ×1.20 at + 2, 3, 4 levels) and **loses on wall clock at every depth** (×0.86, ×0.77, + ×0.55) — the extra sweep costs ~3% and ``gmres`` itself ~40% per cycle. Custom-P + builds its Galerkin coarse operators from barycentric transfers, denser than + the native nested ones, so the smoother is a larger share of the cycle here. + + It is the default anyway, on both routes, because the failure it avoids is + worse than the cost it carries: a stationary smoother stalls where a Krylov one + adapts, and it does so exactly on the badly-conditioned high-contrast problems + this code exists for (the Spiegelman-notch measurement above). Removing a + guardrail is the caller's decision to take; the number is quoted so that + decision can be made. """ return _bundle(_geometric_mg_settings(coarse)) From 5a390efe12e3bee730326f23b36e21b7e39f560a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 11:12:40 +1000 Subject: [PATCH 14/17] The bundle owner must not overwrite what the user set (audit F-1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The multigrid bundle was applied wholesale on every rebuild. The only escape was the `_pc_user_override` latch, and that latch keys on `pc_type` ALONE — so of the bundle's ten keys, setting any of nine was silently discarded, and setting the tenth rescued all ten. Measured, asking for chebyshev / max_it 6 / coarse svd on the Stokes velocity block: what the user set result the three bundle keys only OVERWRITTEN (gmres/4/redundant) the three keys + fieldsplit_velocity_pc_type=mg honoured the three keys + preconditioner="fmg" OVERWRITTEN Note the inversion in the last row: explicit `preconditioner="fmg"` skips the latch entirely, because that mode "always applies". The more clearly a user stated their intent, the less control they had. And the workaround for the first row — "also set pc_type to the value it already has" — is undiscoverable. Making this PR the single owner of the bundle would have entrenched all of it, which is why the fix belongs here. `MGBundle.apply` now takes an `owned` record and leaves alone any key whose current value is not the one UW3 last wrote — in the stale-key clear too, since a user-set key is not ours to remove. Ownership is RECORDED, never inferred from the value: inference fails as soon as a second internal writer touches the same key, which is exactly how the `tolerance` and `strategy` setters defeated an earlier attempt at this (#477). Every internal writer of a bundle key therefore goes through the new `SolverBaseClass._push_managed_option`, so anything unrecorded is the user's by construction. The record is keyed by the GLOBALLY-QUALIFIED option name. The solver writes through a prefixed Options view (`Solver_N_`) while `custom_mg._configure_pcmg` reads the global database using the live PC's own full prefix; recording the unqualified name made every key look user-owned over there and the bundle silently stopped applying. That regression was caught by the defaults arm, not by any test — which is the whole argument for keeping a defaults arm. Tests: `test_user_set_bundle_keys_are_honoured` (parametrised over auto and explicit `fmg`) asserts the two keys the user set survive a rebuild AND that the key they left alone still carries the measured default — respecting one key must not abandon the bundle. `test_unset_bundle_keys_keep_the_managed_defaults` is its control: without it the first test could pass trivially by never applying the bundle, which is precisely how the first attempt failed. Validation: fmg_config_diff shows the three routes still agreeing on every row but the deliberate rotated `svd`. 68 passed across test_1014/1015/1016/1017/1018/1020 and test_0753/0835; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 green (6 / 9 / 2 passed). Found by the solver configuration audit, `docs/reviews/2026-07/solver-configuration-reachability-audit.md`. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 114 ++++++++++++------ src/underworld3/utilities/custom_mg.py | 14 ++- .../utilities/multigrid_options.py | 58 ++++++++- tests/test_1021_mg_option_bundle.py | 56 +++++++++ 4 files changed, 195 insertions(+), 47 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 21c3520b..ced17dd6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -161,11 +161,37 @@ class SolverBaseClass(uw_object): # tell "user set mg" from "we set mg" and would clobber their tuned # smoother / coarse-solver options with the framework FMG bundle. self._pc_user_override = False + # Every multigrid option value UW3 itself has written on the managed block, + # keyed by full option name. This is what lets the bundle honour a + # user-set smoother while still managing the keys the user left alone: a + # present key whose value is not the one we recorded writing is theirs. + # Ownership is RECORDED, never inferred from the value — inference fails + # the moment a second internal writer touches the same key, which is how + # the `tolerance` and `strategy` setters defeated an earlier attempt (#477). + # Internal writers therefore go through _push_managed_option(). + self._managed_pc_options = {} # Custom multigrid prolongation hierarchy (see set_custom_mg / # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + def _push_managed_option(self, key, value): + """Write a PETSc option UW3 owns, recording that we wrote it. + + Use this for any option a solver sets on its own behalf that the multigrid + bundles also write (``utilities.multigrid_options``). A plain + ``self.petsc_options[key] = value`` is indistinguishable from a user's own + write, and the bundle would then back off from a key nobody asked for. + """ + self.petsc_options[key] = value + # Key the record by the GLOBAL option name. `self.petsc_options` is a + # prefixed view (`Solver_N_`), but custom_mg._configure_pcmg reads the + # global database using the live PC's own full prefix — so an unqualified + # record makes every key look user-owned over there and the bundle + # silently stops applying. + self._managed_pc_options[self.petsc_options_prefix + key] = \ + multigrid_options.option_string(value) + @property def consistent_jacobian(self): r"""Jacobian tangent selection: ``False`` | ``True`` | ``"continuation"``. @@ -625,8 +651,18 @@ class SolverBaseClass(uw_object): # be configured from three places (#468). Coarse solve: the native # hierarchy is not rotated, so its coarse operator carries no inherited # null space and redundant+LU is right here. + # `_managed_pc_options` makes the bundle respect a key the USER set while + # still managing the ones they left alone. Before this, the bundle was + # applied wholesale on every rebuild and the only escape was the + # `_pc_user_override` latch above — which keys on `pc_type` ALONE, so a user + # who set (say) `mg_levels_ksp_max_it` had it silently discarded unless they + # also set `pc_type` to the value it already had. Explicit + # `preconditioner="fmg"` was worse: it skips the latch entirely, so the + # clearer the request the less control it carried. if want_fmg: - multigrid_options.geometric_mg_bundle().apply(opts, prefix) + multigrid_options.geometric_mg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -637,7 +673,9 @@ class SolverBaseClass(uw_object): f"mesh with refinement >= 1 to enable geometric multigrid.", stacklevel=2, ) - multigrid_options.gamg_bundle().apply(opts, prefix) + multigrid_options.gamg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "gamg" def _enforce_galerkin_for_geometric_mg(self): @@ -2944,16 +2982,16 @@ class SNES_Scalar(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) self.petsc_options["snes_rtol"] = 1.0e-4 - self.petsc_options["mg_levels_ksp_max_it"] = 3 + self._push_managed_option("mg_levels_ksp_max_it", 3) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -3857,14 +3895,14 @@ class SNES_Vector(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) self.petsc_options["snes_rtol"] = 1.0e-3 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -4869,14 +4907,14 @@ class SNES_MultiComponent(SolverBaseClass): self.petsc_options["snes_type"] = "newtonls" self.petsc_options["ksp_rtol"] = 1.0e-3 self.petsc_options["ksp_type"] = "gmres" - self.petsc_options["pc_type"] = "gamg" - self.petsc_options["pc_gamg_type"] = "agg" - self.petsc_options["pc_gamg_repartition"] = True - self.petsc_options["pc_mg_type"] = "additive" - self.petsc_options["pc_gamg_agg_nsmooths"] = 2 + self._push_managed_option("pc_type", "gamg") + self._push_managed_option("pc_gamg_type", "agg") + self._push_managed_option("pc_gamg_repartition", True) + self._push_managed_option("pc_mg_type", "additive") + self._push_managed_option("pc_gamg_agg_nsmooths", 2) self.petsc_options["snes_rtol"] = 1.0e-3 - self.petsc_options["mg_levels_ksp_max_it"] = 3 - self.petsc_options["mg_levels_ksp_converged_maxits"] = None + self._push_managed_option("mg_levels_ksp_max_it", 3) + self._push_managed_option("mg_levels_ksp_converged_maxits", None) if self.verbose == True: self.petsc_options["ksp_monitor"] = None @@ -5736,13 +5774,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options[f"fieldsplit_{v_name}_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_{v_name}_ksp_max_it"] = 200 self.petsc_options[f"fieldsplit_{v_name}_ksp_rtol"] = self._tolerance * 0.1 - self.petsc_options[f"fieldsplit_{v_name}_pc_type"] = "gamg" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_type"] = "agg" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_repartition"] = True - self.petsc_options[f"fieldsplit_{v_name}_pc_mg_type"] = "additive" - self.petsc_options[f"fieldsplit_{v_name}_pc_gamg_agg_nsmooths"] = 2 - self.petsc_options[f"fieldsplit_{v_name}_mg_levels_ksp_max_it"] = 3 - self.petsc_options[f"fieldsplit_{v_name}_mg_levels_ksp_converged_maxits"] = None + self._push_managed_option(f"fieldsplit_{v_name}_pc_type", "gamg") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_type", "agg") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_repartition", True) + self._push_managed_option(f"fieldsplit_{v_name}_pc_mg_type", "additive") + self._push_managed_option(f"fieldsplit_{v_name}_pc_gamg_agg_nsmooths", 2) + self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_max_it", 3) + self._push_managed_option(f"fieldsplit_{v_name}_mg_levels_ksp_converged_maxits", None) # Create this dict self.fields = {} @@ -6338,16 +6376,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # preconditioning and weakly-indefinite coarse operators; issue #147). self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_velocity_ksp_max_it"] = 200 - self.petsc_options[f"fieldsplit_velocity_pc_type"] = "gamg" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_type"] = "agg" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_repartition"] = True + self._push_managed_option(f"fieldsplit_velocity_pc_type", "gamg") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_type", "agg") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_repartition", True) # NOTE: "kaskade" here diverges from the "additive" set by __init__'s # velocity block — a long-standing (possibly unintentional) difference. # Do not change without benchmarking (READ-23 kept the value as-is). - self.petsc_options[f"fieldsplit_velocity_pc_mg_type"] = "kaskade" - self.petsc_options[f"fieldsplit_velocity_pc_gamg_agg_nsmooths"] = 2 - self.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_max_it"] = 3 - self.petsc_options[f"fieldsplit_velocity_mg_levels_ksp_converged_maxits"] = None + self._push_managed_option(f"fieldsplit_velocity_pc_mg_type", "kaskade") + self._push_managed_option(f"fieldsplit_velocity_pc_gamg_agg_nsmooths", 2) + self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_max_it", 3) + self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_converged_maxits", None) diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 5d96d272..50c8504f 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -595,7 +595,7 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps, coarse="redundant"): +def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -603,7 +603,10 @@ def _configure_pcmg(pc, Ps, coarse="redundant"): the same bundle the native (DMPlex-refinement) route applies — so custom-P and native multigrid cannot be configured differently (#468). ``coarse`` selects the coarse-solve variant: ``"svd"`` on the rotated path, whose - Galerkin-coarsened operator inherits the rigid-rotation null space. + Galerkin-coarsened operator inherits the rigid-rotation null space. ``owned`` + is the solver's record of the option values UW3 has written, so a key the USER + set is left alone (see :meth:`multigrid_options.MGBundle.apply`); ``None`` writes + unconditionally, which is right for the rotated path's per-solve prefix. The bundle is written into the options DB under the PC's OWN options prefix (``pc.getOptionsPrefix()`` — e.g. ``Solver_N_`` for the scalar top-level PC, @@ -617,7 +620,7 @@ def _configure_pcmg(pc, Ps, coarse="redundant"): nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" multigrid_options.geometric_mg_bundle(coarse=coarse).apply( - PETSc.Options(), prefix) + PETSc.Options(), prefix, owned=owned) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) @@ -650,7 +653,8 @@ def _install_transfers(solver, Ps, verbose=False): ksp = solver.snes.getKSP() solver.snes.setUp() ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) - _configure_pcmg(ksp.getPC(), Ps) + _configure_pcmg(ksp.getPC(), Ps, + owned=solver._managed_pc_options) if verbose: from underworld3 import mpi mpi.pprint(f"[{solver.name}] custom FMG installed: {nlev} levels, " @@ -711,7 +715,7 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): vel_ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) - _configure_pcmg(vel_pc, Ps) + _configure_pcmg(vel_pc, Ps, owned=solver._managed_pc_options) vel_pc.setUp() # 4. re-attach the coupled Stokes nullspace (operator state was touched) diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index cd2e9d27..6ecfa33e 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -49,10 +49,32 @@ from typing import NamedTuple -__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", +__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", "option_string", "GEOMETRIC_MG_COARSE_SOLVERS"] +def option_string(value): + """How PETSc stores ``value`` in its options database, so a value UW3 wrote can + be compared against what is in there now. ``None`` is a bare flag (reads back as + an empty string); booleans lower-case.""" + if value is None: + return "" + if isinstance(value, bool): + return str(value).lower() + return str(value) + + +def _user_owns(opts, name, owned): + """Was ``name`` set by someone other than UW3's own option writers? + + True when the key is present and its current value is not the one UW3 last + recorded writing. ``owned is None`` disables the check (the caller owns the + prefix outright).""" + if owned is None or not opts.hasName(name): + return False + return opts.getString(name) != owned.get(name) + + #: Coarse-solve variants of the geometric bundle. ``"redundant"`` (redundant+LU) #: is the default and is np-safe; ``"svd"`` is required whenever the coarse #: operator inherits a null space — see :func:`geometric_mg_bundle`. @@ -71,17 +93,45 @@ class MGBundle(NamedTuple): settings: dict stale: tuple - def apply(self, opts, prefix=""): + def apply(self, opts, prefix="", owned=None): """Write this bundle into ``opts`` (a ``PETSc.Options``) under ``prefix``. ``prefix`` is the managed block's own prefix — ``""`` for a top-level PC, ``"fieldsplit_velocity_"`` for the Stokes velocity sub-block — and is applied on top of whatever prefix ``opts`` itself carries. + + Parameters + ---------- + owned : dict or None + Bookkeeping of the option values UW3 itself has written, keyed by full + option name. A key that is already present but whose current value is + **not** what UW3 last wrote was set by the user, and is left alone — + both here and in the stale-key clear, because it is not ours to remove. + Pass ``None`` to write unconditionally, which is right when the caller + owns the whole prefix (the rotated path builds a fresh, per-solve + prefix no user can address). + + Ownership is *recorded*, never inferred from the value. Guessing by + value fails as soon as a second internal writer touches the same key — + which is exactly how the ``tolerance`` and ``strategy`` setters defeated + an earlier attempt (#477). Every internal writer of a bundle key goes + through this bookkeeping (``SolverBaseClass._push_managed_option``), so + anything unrecorded is the user's by construction. """ for key, value in self.settings.items(): - opts.setValue(prefix + key, value) + name = prefix + key + if _user_owns(opts, name, owned): + continue + opts.setValue(name, value) + if owned is not None: + owned[name] = option_string(value) for key in self.stale: - opts.delValue(prefix + key) + name = prefix + key + if _user_owns(opts, name, owned): + continue + opts.delValue(name) + if owned is not None: + owned.pop(name, None) def _geometric_mg_settings(coarse): diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 53b8983d..c99c09a3 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -201,6 +201,62 @@ def test_rotated_picks_up_a_mesh_owned_hierarchy(): assert s._rotated_freeslip_info["velocity_pc_type"] == "mg" +def _velocity_mg_config(s): + """(smoother ksp, smoother max_it, coarse pc) off the live velocity sub-PC.""" + pc = s.snes.getKSP().getPC() + pc.setUp() # MUST precede getFieldSplitSubKSP + vpc = pc.getFieldSplitSubKSP()[0].getPC() + assert vpc.getType() == "mg" and vpc.getMGLevels() > 1, ( + f"expected geometric MG on the velocity block, got {vpc.getType()}") + smoother = vpc.getMGSmoother(vpc.getMGLevels() - 1) + return (smoother.getType(), smoother.getTolerances()[3], + vpc.getMGCoarseSolve().getPC().getType()) + + +@pytest.mark.parametrize("preconditioner", [None, "fmg"]) +def test_user_set_bundle_keys_are_honoured(preconditioner): + """A user who sets a bundle key must get it, and must keep the managed values + for the keys they left alone — per key, not all-or-nothing. + + The bundle used to be applied wholesale on every rebuild. The only escape was + the ``_pc_user_override`` latch, which keys on ``pc_type`` ALONE, so setting + ``mg_levels_ksp_max_it`` was silently discarded unless the user also set + ``pc_type`` to the value it already had. Explicit ``preconditioner="fmg"`` was + worse — it skips that latch entirely, so the clearer the request the less + control it carried, which is why both modes are parametrised here. + """ + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, f"u{'f' if preconditioner else 'a'}", rotated=False) + if preconditioner is not None: + s.preconditioner = preconditioner + # set TWO of the bundle's keys and leave the rest to the framework + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 6 + s.petsc_options["fieldsplit_velocity_mg_coarse_pc_type"] = "svd" + s.solve() + s.solve() # a rebuild must not clobber it + + smoother, max_it, coarse = _velocity_mg_config(s) + assert max_it == 6, "user-set smoother iteration count was overwritten" + assert coarse == "svd", "user-set coarse solver was overwritten" + # ...and the key the user did NOT set still carries the measured default + assert smoother == "gmres", ( + "respecting one bundle key must not abandon the rest of the bundle") + + +def test_unset_bundle_keys_keep_the_managed_defaults(): + """The control for the test above: with nothing set, the whole managed bundle + applies. Without this, 'user-set keys are honoured' could pass trivially by + never applying the bundle at all — which is exactly how the first attempt at + this failed (the routes silently reverted to GAMG).""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "dflt", rotated=False) + s.solve() + s.solve() + assert _velocity_mg_config(s) == ("gmres", 4, "redundant") + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under From b4159f5125a0c54a695e3a6e70e7e0239bfeda76 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 11:59:36 +1000 Subject: [PATCH 15/17] Fill solver.strategy from the measured smoother regimes; rename MGBundle to MGSettings MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `solver.strategy` has been a public property accepting "default", "robust" and "fast" since it was written, validating your input and then configuring all three identically — its docstring said "(Reserved)". A property that checks a value and then ignores it is the same defect class as #477 and #478: the failure is invisible because the solve still converges. This fills it, and it turns out the machinery this PR built is exactly what was missing underneath. TWO LAYERS, which is the design point: solver.strategy the named intent public, already existed MGSettings the option values utilities/multigrid_options.py solver.petsc_options the escape hatch public, outranks both (F-1) `geometric_mg_bundle(smoother=...)` gains the variant axis, from the two regimes measured in this PR's own work — neither dominates: "robust" = gmres/4 survives an operator a stationary smoother stalls on. Spiegelman notch (eta contrast 1e26, 4 levels) per-V-cycle contraction 0.56 vs richardson's 0.75, margin growing with depth; TI rotated annulus 11 velocity its -> 5, 1.74x. "fast" = richardson/3 cheaper per cycle where the operator is benign. Linear, SYMMETRIC annulus at eta contrast 1e6: beats robust on wall clock at every depth (x1.16, x1.30, x1.82 at 2, 3, 4 levels) while taking more iterations. "default" stays "robust", so nobody's results move. "fast" is the documented opt-out for the cost the adversarial review on this PR quoted and could not otherwise offer a way out of. ALSO REMOVES DEAD CODE THE MEASUREMENT EXPOSED. The `strategy` setter wrote the whole GAMG velocity bundle plus `pc_mg_type=kaskade`, the latter carrying a comment warning against changing it without benchmarking. Every one of those writes was DEAD: `_apply_preconditioner_options` runs later, at `_build`, and is the single writer of that block — it overwrites them with the geometric bundle on a refined mesh, or the GAMG bundle ("additive") without one. Measured in both orders and with `preconditioner="fmg"` set before and after: the live PC was mg/FULL every time, so `kaskade` never once took effect. The strategy's effect on the velocity block now runs through `_mg_smoother_variant`, which selects a bundle variant rather than racing the bundle writer. `MGBundle` -> `MGSettings` (internal, no call site outside the module and custom_mg): "MG" is what the module actually holds, and "Bundle" said nothing. Tests: strategy selects a real variant; "default" reproduces the framework default exactly (the control — without it, filling the axis could silently move results); a user-set option still beats the strategy while the strategy keeps the keys the user left alone (the two layers compose in the right order). Validation: three routes still agree bar the deliberate rotated svd; 79 passed across the MG, adaptivity and solve-report suites; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 green (6 / 9 / 2). Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 43 ++++++--- .../cython/petsc_generic_snes_solvers.pyx | 85 +++++++++++----- src/underworld3/utilities/custom_mg.py | 2 +- .../utilities/multigrid_options.py | 96 ++++++++++++------- tests/test_1021_mg_option_bundle.py | 45 +++++++++ 5 files changed, 201 insertions(+), 70 deletions(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index ba2e086b..eda4c33a 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -418,26 +418,47 @@ stationary iteration stalls where a Krylov smoother adapts its polynomial. It predicts the notch result, and it predicts the cost below on well-conditioned problems. -### The smoother is not free +### Two regimes, and the named intent that selects them -On a *well-conditioned* problem gmres/4 is slower than richardson/3 despite -converging in fewer iterations. Nested annulus, linear (symmetric) velocity -block, η contrast 1e6, outer KSP timed in isolation: +Neither smoother dominates, so there are two variants and `solver.strategy` +chooses between them. This is the layering: + +| layer | what it is | where | +|---|---|---| +| `solver.strategy` | the **named intent** — `"default"` / `"robust"` / `"fast"` | public property on the solver | +| `MGSettings` | the **option values** an intent resolves to | `utilities/multigrid_options.py` | +| `solver.petsc_options` | the **escape hatch** — any key you set here outranks both | public, unchanged | + +`"robust"` (gmres/4) survives an operator a stationary smoother stalls on: the +Spiegelman notch contraction above, and 11 → 5 velocity iterations on the +transversely isotropic rotated annulus. + +`"fast"` (richardson/3) is cheaper per cycle and quicker where the operator is +benign. Nested annulus, **linear (symmetric)** velocity block, η contrast 1e6, +outer KSP timed in isolation — `"fast"` against `"robust"`: | depth | iterations | wall clock | |---|---|---| -| 2 levels | ×1.45 | **×0.86** | -| 3 levels | ×1.60 | **×0.77** | -| 4 levels | ×1.20 | **×0.55** | +| 2 levels | ×0.69 | **×1.16** | +| 3 levels | ×0.63 | **×1.30** | +| 4 levels | ×0.83 | **×1.82** | The extra sweep costs ~3%; gmres itself costs ~40% per cycle. Custom-P builds its Galerkin coarse operators from barycentric transfers, denser than the native nested ones, so the smoother is a larger share of the cycle on that route. -It stays the default on both routes because the failure it avoids (a stationary -smoother stalling on a high-contrast operator) is worse than the cost it carries, -and it carries that cost exactly where the problems are easy. Removing a guardrail -is the caller's decision; the number is recorded so the decision can be made. +`"default"` is `"robust"`: the failure it avoids is worse than the cost it carries, +and it carries that cost exactly where the problem is easy. `"fast"` is the +documented opt-out. + +```{warning} +Do not fill a strategy by writing velocity-block options from the strategy setter. +`_apply_preconditioner_options` runs later, at `_build`, and is the single writer +of that block — anything written earlier is overwritten. That is how +`pc_mg_type=kaskade` sat in the `strategy` setter for a long time carrying a +comment warning against changing it, while never once taking effect (measured: the +live PC was `mg`/FULL in every ordering). A strategy selects a bundle *variant*. +``` ### The one legitimate per-route difference diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index ced17dd6..3abc62af 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -175,6 +175,16 @@ class SolverBaseClass(uw_object): # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + @property + def _mg_smoother_variant(self): + """Which measured smoother regime this solver's strategy asks for. + + ``solver.strategy`` is the named intent ("I want speed" / "I want this to + converge"); the values live in ``utilities.multigrid_options``. Solvers with + no strategy axis get the robust default. See + :func:`multigrid_options.geometric_mg_bundle` for the measurements.""" + return "fast" if getattr(self, "_strategy", "default") == "fast" else "robust" + def _push_managed_option(self, key, value): """Write a PETSc option UW3 owns, recording that we wrote it. @@ -660,9 +670,10 @@ class SolverBaseClass(uw_object): # `preconditioner="fmg"` was worse: it skips the latch entirely, so the # clearer the request the less control it carried. if want_fmg: - multigrid_options.geometric_mg_bundle().apply( - PETSc.Options(), self.petsc_options_prefix + prefix, - owned=self._managed_pc_options) + multigrid_options.geometric_mg_bundle( + smoother=self._mg_smoother_variant).apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) self._pc_managed_value = "mg" else: if self._preconditioner == "fmg" and n_levels <= 1 and uw.mpi.rank == 0: @@ -6309,19 +6320,39 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): @property def strategy(self): """ - Solver strategy controlling preconditioner configuration. - - Currently supports: - - ``"default"``: Standard Schur complement fieldsplit with GAMG - - ``"robust"``: (Reserved) More robust but slower configuration - - ``"fast"``: (Reserved) Faster but less robust configuration - - Setting this property reconfigures the entire preconditioner stack. + What this solve should optimise for — the named intent over the + multigrid smoother's two measured regimes. + + - ``"default"``, ``"robust"``: ``gmres``/4 smoothing. Survives an operator a + stationary smoother stalls on: Spiegelman notch (:math:`\eta` contrast + 1e26, 4 levels) per-V-cycle contraction 0.56 against richardson's 0.75, the + margin growing with depth; transversely isotropic rotated annulus 11 + velocity iterations down to 5. + - ``"fast"``: ``richardson``/3 smoothing. Cheaper per cycle and quicker where + the operator is benign — on a linear, symmetric annulus at + :math:`\eta` contrast 1e6 it beats ``"robust"`` on wall clock at every + hierarchy depth tested (x1.16, x1.30, x1.82 at 2, 3, 4 levels) while taking + more iterations. It gives up the regime ``"robust"`` exists for, so it is + an opt-in. + + ``"default"`` is ``"robust"``: the failure it avoids is worse than the cost it + carries, and it carries that cost exactly where the problem is easy. + + Setting this property also resets the fieldsplit / Schur / pressure sub-solve + configuration to the framework defaults. It does **not** write the velocity + block's preconditioner directly — that is applied later, from + :mod:`underworld3.utilities.multigrid_options`, which is the single writer of + those options; this property selects which variant it applies. Values written + by hand into :attr:`petsc_options` are respected and not overwritten. Returns ------- str Current strategy name. + + See Also + -------- + preconditioner : which multigrid FAMILY to use (geometric, algebraic, auto). """ return self._strategy @@ -6332,14 +6363,19 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): f"Unknown solver strategy {value!r}: " "expected 'default', 'robust', or 'fast'." ) - # 'robust' and 'fast' are accepted names but currently configure the - # same option bundle as 'default' (their dedicated branches were empty - # placeholders and have been removed — Charter S5, READ-23). + # 'fast' and 'robust' now select a real smoother variant, via + # `_mg_smoother_variant` -> `multigrid_options.geometric_mg_bundle`. They were + # accepted-and-inert placeholders for a long time: validated on input, then + # configured identically to 'default'. A property that checks your value and + # then ignores it is the same defect class as #477 and #478 — the failure is + # invisible, because the solve still converges. # self.is_setup = False self._strategy = value - # All strategies: reset to preferred + # Common to every strategy: reset the fieldsplit / Schur / pressure + # sub-solve to the framework defaults. The strategy's effect on the VELOCITY + # BLOCK is carried by `_mg_smoother_variant`, not by writes from here. self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 @@ -6376,16 +6412,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # preconditioning and weakly-indefinite coarse operators; issue #147). self.petsc_options[f"fieldsplit_velocity_ksp_type"] = "fgmres" self.petsc_options[f"fieldsplit_velocity_ksp_max_it"] = 200 - self._push_managed_option(f"fieldsplit_velocity_pc_type", "gamg") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_type", "agg") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_repartition", True) - # NOTE: "kaskade" here diverges from the "additive" set by __init__'s - # velocity block — a long-standing (possibly unintentional) difference. - # Do not change without benchmarking (READ-23 kept the value as-is). - self._push_managed_option(f"fieldsplit_velocity_pc_mg_type", "kaskade") - self._push_managed_option(f"fieldsplit_velocity_pc_gamg_agg_nsmooths", 2) - self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_max_it", 3) - self._push_managed_option(f"fieldsplit_velocity_mg_levels_ksp_converged_maxits", None) + # The velocity BLOCK's preconditioner is not set here. `strategy` used to + # write the whole GAMG bundle plus `pc_mg_type=kaskade`, and every one of + # those writes was DEAD: `_apply_preconditioner_options` runs later (at + # `_build`) and overwrites them with the geometric bundle on a refined mesh, + # or the GAMG bundle (`additive`) without one. Measured both orders — the + # live PC was `mg`/FULL every time, so `kaskade` never once took effect + # despite the comment warning against changing it. The strategy's effect on + # the velocity block now runs through `_mg_smoother_variant`, which selects + # a bundle variant instead of racing the bundle writer. diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 50c8504f..588b0031 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -605,7 +605,7 @@ def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): selects the coarse-solve variant: ``"svd"`` on the rotated path, whose Galerkin-coarsened operator inherits the rigid-rotation null space. ``owned`` is the solver's record of the option values UW3 has written, so a key the USER - set is left alone (see :meth:`multigrid_options.MGBundle.apply`); ``None`` writes + set is left alone (see :meth:`multigrid_options.MGSettings.apply`); ``None`` writes unconditionally, which is right for the rotated path's per-solve prefix. The bundle is written into the options DB under the PC's OWN options prefix diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 6ecfa33e..4f2700a5 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -38,7 +38,9 @@ Apply the geometric bundle to a solver's managed block:: from underworld3.utilities import multigrid_options - multigrid_options.geometric_mg_bundle().apply(self.petsc_options, prefix) + multigrid_options.geometric_mg_bundle().apply( + PETSc.Options(), self.petsc_options_prefix + prefix, + owned=self._managed_pc_options) Read the settings without writing them (the rotated path stages its options in a dict so it can drop them again after ``setUp``):: @@ -49,8 +51,8 @@ from typing import NamedTuple -__all__ = ["MGBundle", "geometric_mg_bundle", "gamg_bundle", "option_string", - "GEOMETRIC_MG_COARSE_SOLVERS"] +__all__ = ["MGSettings", "geometric_mg_bundle", "gamg_bundle", "option_string", + "GEOMETRIC_MG_COARSE_SOLVERS", "GEOMETRIC_MG_SMOOTHERS"] def option_string(value): @@ -80,8 +82,14 @@ def _user_owns(opts, name, owned): #: operator inherits a null space — see :func:`geometric_mg_bundle`. GEOMETRIC_MG_COARSE_SOLVERS = ("redundant", "svd") +#: Smoother variants. These are the two measured regimes, not a taste setting: +#: ``"robust"`` survives a badly-conditioned operator where a stationary smoother +#: stalls, ``"fast"`` is cheaper per cycle where the operator is benign. Selected by +#: ``solver.strategy``; see :func:`geometric_mg_bundle` for the numbers. +GEOMETRIC_MG_SMOOTHERS = ("robust", "fast") -class MGBundle(NamedTuple): + +class MGSettings(NamedTuple): """A preconditioner option bundle: the keys to set, and the keys to clear. ``settings`` maps an option suffix (no prefix, no leading ``-``) to its @@ -134,8 +142,11 @@ def apply(self, opts, prefix="", owned=None): owned.pop(name, None) -def _geometric_mg_settings(coarse): - """The geometric-MG settings for one coarse-solve variant.""" +def _geometric_mg_settings(coarse, smoother="robust"): + """The geometric-MG settings for one coarse-solve and smoother variant. + + Both variants set the SAME keys — only values differ — so the derived stale-key + sets are variant-independent.""" settings = { "pc_type": "mg", "pc_mg_type": "full", # FMG (F-cycle) @@ -173,12 +184,7 @@ def _geometric_mg_settings(coarse): # spread far enough that a stationary iteration stalls where a Krylov # smoother adapts its polynomial. That also predicts the cost measured on # well-conditioned problems (see the note on geometric_mg_bundle). - "mg_levels_ksp_type": "gmres", "mg_levels_pc_type": "sor", - # SET the count, never inherit it. PCMG's own default is 2 and the GAMG - # bundle below leaves 3 under the same prefix, so a bundle that omits - # this key smooths differently depending on what ran before it (#468). - "mg_levels_ksp_max_it": 4, # Run EXACTLY max_it smoother iterations: no residual-norm computation and # no convergence test, so every V-cycle costs the same. A Krylov smoother # makes the cycle non-stationary, which is why the velocity block is @@ -186,6 +192,20 @@ def _geometric_mg_settings(coarse): "mg_levels_ksp_norm_type": "none", "mg_levels_ksp_converged_maxits": None, } + if smoother == "robust": + settings["mg_levels_ksp_type"] = "gmres" + settings["mg_levels_ksp_max_it"] = 4 + elif smoother == "fast": + # Cheaper per cycle and measurably faster where the operator is benign, + # at the cost of the regime "robust" exists for. See the docstring. + settings["mg_levels_ksp_type"] = "richardson" + settings["mg_levels_ksp_max_it"] = 3 + else: + raise ValueError( + f"smoother must be one of {GEOMETRIC_MG_SMOOTHERS} (got {smoother!r})") + # SET the iteration count, never inherit it. PCMG's own default is 2 and the + # GAMG bundle leaves 3 under the same prefix, so a bundle that omits this key + # smooths differently depending on what ran before it (#468). if coarse == "redundant": # redundant+lu, not bare lu: a bare serial LU cannot factor a distributed # coarse matrix and fails at np>1 (DIVERGED_LINEAR_SOLVE after 0 @@ -231,15 +251,19 @@ def _all_keys(): def _bundle(settings): - return MGBundle(settings=settings, + return MGSettings(settings=settings, stale=tuple(sorted(_all_keys() - set(settings)))) -def geometric_mg_bundle(coarse="redundant"): +def geometric_mg_bundle(coarse="redundant", smoother="robust"): """Geometric multigrid (FMG F-cycle) on an explicit hierarchy. Parameters ---------- + smoother : {"robust", "fast"} + Which of the two MEASURED regimes to configure. ``"robust"`` is + ``gmres``/4 and is the default; ``"fast"`` is ``richardson``/3. Chosen by + ``solver.strategy``, not usually here. coarse : {"redundant", "svd"} The coarse-level solve. ``"redundant"`` (redundant+LU) is the default. Use ``"svd"`` when the coarse operator inherits a null space — which the @@ -248,28 +272,34 @@ def geometric_mg_bundle(coarse="redundant"): Returns ------- - MGBundle - Settings and stale keys; see :meth:`MGBundle.apply`. + MGSettings + Settings and stale keys; see :meth:`MGSettings.apply`. Notes ----- - The ``gmres``/4 smoother is not free on a well-conditioned problem. Measured on - a nested annulus with a linear (symmetric) velocity block at :math:`\eta` - contrast 1e6, outer KSP timed in isolation, ``gmres``/4 against - ``richardson``/3: it wins on iterations at every depth (×1.45, ×1.60, ×1.20 at - 2, 3, 4 levels) and **loses on wall clock at every depth** (×0.86, ×0.77, - ×0.55) — the extra sweep costs ~3% and ``gmres`` itself ~40% per cycle. Custom-P - builds its Galerkin coarse operators from barycentric transfers, denser than - the native nested ones, so the smoother is a larger share of the cycle here. - - It is the default anyway, on both routes, because the failure it avoids is - worse than the cost it carries: a stationary smoother stalls where a Krylov one - adapts, and it does so exactly on the badly-conditioned high-contrast problems - this code exists for (the Spiegelman-notch measurement above). Removing a - guardrail is the caller's decision to take; the number is quoted so that - decision can be made. + The two smoother variants are two measured regimes, and neither dominates. + + ``"robust"`` (``gmres``/4) survives an operator a stationary smoother stalls on. + Spiegelman notch, :math:`\eta` contrast 1e26, nested 4-level hierarchy: + per-V-cycle contraction :math:`\rho` = 0.56 against richardson's 0.75 at the same + four iterations, the margin growing with depth. Transversely isotropic rotated + annulus: 11 velocity iterations down to 5, and 1.74x on the linear solve. + + ``"fast"`` (``richardson``/3) is cheaper per cycle and measurably quicker where + the operator is benign. Nested annulus, LINEAR (symmetric) velocity block at + :math:`\eta` contrast 1e6, outer KSP timed in isolation: robust wins on + iterations at every depth (x1.45, x1.60, x1.20 at 2, 3, 4 levels) and **loses on + wall clock at every depth** (x0.86, x0.77, x0.55). The extra sweep costs ~3% and + ``gmres`` itself ~40% per cycle; custom-P builds its Galerkin coarse operators + from barycentric transfers, denser than the native nested ones, so the smoother + is a larger share of the cycle on that route. + + ``"robust"`` is the default because the failure it avoids is worse than the cost + it carries, and it carries that cost exactly where the problem is easy. Choose + between them with ``solver.strategy``, which is the named-intent layer over + these values. """ - return _bundle(_geometric_mg_settings(coarse)) + return _bundle(_geometric_mg_settings(coarse, smoother)) def gamg_bundle(): @@ -277,7 +307,7 @@ def gamg_bundle(): Returns ------- - MGBundle - Settings and stale keys; see :meth:`MGBundle.apply`. + MGSettings + Settings and stale keys; see :meth:`MGSettings.apply`. """ return _bundle(_gamg_settings()) diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index c99c09a3..9934e252 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -257,6 +257,51 @@ def test_unset_bundle_keys_keep_the_managed_defaults(): assert _velocity_mg_config(s) == ("gmres", 4, "redundant") +def _strategy_mg_config(strategy): + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, f"s{strategy or 'none'}"[:6], rotated=False) + if strategy is not None: + s.strategy = strategy + s.solve() + s.solve() + return _velocity_mg_config(s) + + +def test_strategy_selects_a_real_smoother_variant(): + """``solver.strategy`` must actually change the smoother. + + ``"fast"`` and ``"robust"`` were accepted-and-inert for a long time: validated + on input, then configured identically to ``"default"``. A property that checks + your value and then ignores it is the same defect class as #477/#478 — invisible, + because the solve still converges. They now select a measured variant: + ``"robust"`` is gmres/4, ``"fast"`` is richardson/3. + """ + assert _strategy_mg_config("fast")[:2] == ("richardson", 3) + assert _strategy_mg_config("robust")[:2] == ("gmres", 4) + + +def test_default_strategy_does_not_change_behaviour(): + """The control: ``"default"`` must reproduce the framework default exactly, so + filling the strategy axis moves nobody's results.""" + assert _strategy_mg_config("default") == _strategy_mg_config(None) + + +def test_a_user_key_still_beats_the_strategy(): + """The two layers compose in the right order: an explicit option the user wrote + outranks the strategy's choice of variant.""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "sxu", rotated=False) + s.strategy = "fast" # would ask for richardson/3 + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_type"] = "chebyshev" + s.solve() + s.solve() + smoother, max_it, _ = _velocity_mg_config(s) + assert smoother == "chebyshev", "the user's smoother lost to the strategy" + assert max_it == 3, "the strategy should still own the keys the user left alone" + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under From 13bd82ceb6ec05952d544b16bde2fb0dbc6a0a24 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 12:27:12 +1000 Subject: [PATCH 16/17] solver.strategy reports the preconditioner it resolved to (#484 shape) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Asking "what am I actually running?" required knowing which nine PETSc option keys to look up, and `_pc_managed_value` / `_pc_user_override` were private. That is the gap the configuration audit filed as #484 — of twelve preconditioner fallbacks, ten left no readable state — and `solver.strategy` is the natural place to close it for this block, since the strategy is what the user set. The getter now returns the strategy NAME that also reports what it resolved to: >>> stokes.strategy 'default' — geometric multigrid (2 levels), full cycle, smoother gmresx4 + sor, coarse redundant/lu >>> stokes.strategy == "default" True `_StrategyName` subclasses `str` deliberately, so comparisons, formatting and serialisation are unchanged — nothing that used `solver.strategy` before behaves differently. Verified nothing in src/ or tests/ compares it to a string in the first place; they only set it. Three properties of the report matter more than its wording: - It NAMES a user override rather than absorbing it, so the summary cannot hide the difference between what the strategy asked for and what is running. - It REFUSES to report before it knows. The preconditioner resolves at the first solve; until then the report says so instead of presenting the constructor defaults as the answer. The first version of this did present them, which is precisely the stale-but-authoritative summary this reporting exists to prevent — hence the new `_pc_resolved` flag. - `solver.preconditioner_settings` is the same information as a dict, so a test can assert on it instead of parsing prose or inferring from timings. `multigrid_options.describe()` does the formatting, next to the values it describes, so a new bundle variant cannot be added without the report following it. Validation: 14 passed in test_1021; 49 passed across test_1014/1017/1018/1020/1058; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 15 passed each. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/solvers.md | 30 ++++++ .../cython/petsc_generic_snes_solvers.pyx | 98 ++++++++++++++++++- .../utilities/multigrid_options.py | 38 +++++++ tests/test_1021_mg_option_bundle.py | 42 ++++++++ 4 files changed, 207 insertions(+), 1 deletion(-) diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index eda4c33a..e03c9309 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -451,6 +451,36 @@ nested ones, so the smoother is a larger share of the cycle on that route. and it carries that cost exactly where the problem is easy. `"fast"` is the documented opt-out. +### Asking what you actually got + +`solver.strategy` reports it. The value is still the strategy *name* — it compares +and formats as the plain string, so nothing that used it before changes — but +displaying it shows the preconditioner that name resolved to: + +```python +>>> stokes.strategy +'default' — geometric multigrid (2 levels), full cycle, smoother gmresx4 + sor, + coarse redundant/lu +>>> stokes.strategy == "default" +True +``` + +Three properties of that report matter more than its formatting: + +- **It names a user override rather than absorbing it.** A key you set yourself is + listed as such (`; overridden by the user: mg_levels_ksp_max_it=6`), so the + summary cannot hide the difference between what the strategy asked for and what + is running. +- **It refuses to report before it knows.** The preconditioner is resolved at the + first solve; until then the report says so instead of presenting the constructor + defaults as though they were the answer. A summary that looks authoritative and + is stale is worse than no summary. +- **`solver.preconditioner_settings`** is the same information as a dict, so a test + can assert on it rather than parsing prose or inferring from timings. + +This is the general shape #484 asks for across all the managed fallbacks: report the +resolved choice as readable state, and distinguish "could not" from "chose not to". + ```{warning} Do not fill a strategy by writing velocity-block options from the strategy setter. `_apply_preconditioner_options` runs later, at `_build`, and is the single writer diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3abc62af..715c78f1 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -14,6 +14,31 @@ import underworld3.timing as timing from underworld3.utilities._api_tools import uw_object from underworld3.utilities import multigrid_options + +class _StrategyName(str): + """A strategy name that also reports what it resolved to. + + Subclasses ``str`` deliberately: ``solver.strategy == "fast"``, string + formatting and serialisation all behave exactly as before, but displaying it — + in a REPL, a notebook, or a log line — shows the preconditioner it actually + configured. Asking "what am I running?" should not require knowing which nine + PETSc option keys to look up. + + Use :attr:`SolverBaseClass.preconditioner_settings` for the machine-readable + form. + """ + + def __new__(cls, name, summary): + obj = super().__new__(cls, name) + obj._summary = summary + return obj + + def __repr__(self): + return f"{str.__repr__(self)} — {self._summary}" + + def _repr_markdown_(self): + return f"**`{str(self)}`** — {self._summary}" + from underworld3.function import expression as public_expression expression = lambda *x, **X: public_expression(*x, _unique_name_generation=True, **X) @@ -170,11 +195,59 @@ class SolverBaseClass(uw_object): # the `tolerance` and `strategy` setters defeated an earlier attempt (#477). # Internal writers therefore go through _push_managed_option(). self._managed_pc_options = {} + # Has _apply_preconditioner_options actually made the resolution decision + # yet? Until it has, the options database holds only this solver's __init__ + # defaults, which is NOT what the next solve will run — reporting them as + # resolved would be exactly the stale-but-authoritative-looking summary this + # reporting exists to prevent. + self._pc_resolved = False # Custom multigrid prolongation hierarchy (see set_custom_mg / # utilities.custom_mg). None => standard FMG/GAMG path, unchanged. self._custom_mg = None + @property + def preconditioner_settings(self): + """The option values the managed preconditioner block is configured with. + + A read-only dict of the multigrid keys UW3 currently has in the options + database for this solver's managed block, so what actually got applied can + be asserted on instead of inferred from timings. Empty for a solver with no + managed block (``_pc_option_prefix is None``), or before the first build + resolves one. + + The values are what the strategy resolved to, *including* any key you set + yourself — those are respected (see :attr:`petsc_options`). Use + :attr:`strategy` for a readable summary of the same thing. + """ + prefix = self._pc_option_prefix + if prefix is None: + return {} + keys = set(multigrid_options.gamg_bundle().settings) + for coarse in multigrid_options.GEOMETRIC_MG_COARSE_SOLVERS: + keys |= set(multigrid_options.geometric_mg_bundle(coarse=coarse).settings) + out = {} + for key in sorted(keys): + name = prefix + key + if self.petsc_options.hasName(name): + out[key] = self.petsc_options.getString(name) + return out + + @property + def _user_overridden_pc_options(self): + """The managed-block keys the USER set, as (key, value) pairs. + + A key present in the options database whose value is not the one UW3 + recorded writing is theirs — the same test the bundle writer uses to decide + what to leave alone.""" + prefix = self._pc_option_prefix + if prefix is None: + return () + qualified = self.petsc_options_prefix + return tuple( + (key, value) for key, value in self.preconditioner_settings.items() + if self._managed_pc_options.get(qualified + prefix + key) != value) + @property def _mg_smoother_variant(self): """Which measured smoother regime this solver's strategy asks for. @@ -591,6 +664,9 @@ class SolverBaseClass(uw_object): prefix = self._pc_option_prefix if prefix is None: return + # Every path from here is a resolution decision, including "the user owns + # these options, leave them alone". + self._pc_resolved = True opts = self.petsc_options @@ -6350,11 +6426,31 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): str Current strategy name. + The value returned is the strategy name (it compares and formats as the + plain string) and additionally reports the preconditioner it resolved to + when displayed:: + + >>> stokes.strategy + 'default' — geometric multigrid (3 levels), full cycle, + smoother gmresx4 + sor, coarse redundant/lu + See Also -------- preconditioner : which multigrid FAMILY to use (geometric, algebraic, auto). + preconditioner_settings : the same information as a dict, for assertions. """ - return self._strategy + settings = self.preconditioner_settings + if not settings or not self._pc_resolved: + summary = "not resolved yet — configured at the first solve" + if settings: + summary += (f" (framework defaults in place: " + f"{multigrid_options.describe(settings)})") + else: + levels = len(getattr(self.mesh, "dm_hierarchy", []) or []) or None + summary = multigrid_options.describe( + settings, levels=levels, + overridden=self._user_overridden_pc_options) + return _StrategyName(self._strategy, summary) @strategy.setter def strategy(self, value): diff --git a/src/underworld3/utilities/multigrid_options.py b/src/underworld3/utilities/multigrid_options.py index 4f2700a5..2ba90a26 100644 --- a/src/underworld3/utilities/multigrid_options.py +++ b/src/underworld3/utilities/multigrid_options.py @@ -52,6 +52,7 @@ from typing import NamedTuple __all__ = ["MGSettings", "geometric_mg_bundle", "gamg_bundle", "option_string", + "describe", "GEOMETRIC_MG_COARSE_SOLVERS", "GEOMETRIC_MG_SMOOTHERS"] @@ -142,6 +143,43 @@ def apply(self, opts, prefix="", owned=None): owned.pop(name, None) +def describe(settings, levels=None, overridden=()): + """One line of plain English for a resolved bundle. + + Parameters + ---------- + settings : dict + The bundle's settings, as applied. + levels : int or None + Hierarchy depth, when known. + overridden : sequence of (key, value) + Keys the user has taken over, which the bundle left alone. + """ + pc = settings.get("pc_type", "?") + if pc == "mg": + what = "geometric multigrid" + if levels: + what += f" ({levels} levels)" + cycle = settings.get("pc_mg_type", "?") + smoother = (f"{settings.get('mg_levels_ksp_type', '?')}" + f"x{settings.get('mg_levels_ksp_max_it', '?')}" + f" + {settings.get('mg_levels_pc_type', '?')}") + coarse = settings.get("mg_coarse_pc_type", "?") + if coarse == "redundant": + coarse += f"/{settings.get('mg_coarse_redundant_pc_type', '?')}" + text = f"{what}, {cycle} cycle, smoother {smoother}, coarse {coarse}" + elif pc == "gamg": + text = (f"algebraic multigrid (gamg/{settings.get('pc_gamg_type', '?')}), " + f"{settings.get('pc_mg_type', '?')} cycle, " + f"smoother max_it {settings.get('mg_levels_ksp_max_it', '?')}") + else: + text = f"pc_type={pc}" + if overridden: + text += ("; overridden by the user: " + + ", ".join(f"{k}={v}" for k, v in overridden)) + return text + + def _geometric_mg_settings(coarse, smoother="robust"): """The geometric-MG settings for one coarse-solve and smoother variant. diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 9934e252..7468388d 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -302,6 +302,48 @@ def test_a_user_key_still_beats_the_strategy(): assert max_it == 3, "the strategy should still own the keys the user left alone" +def test_strategy_reports_what_it_resolved_to(): + """``solver.strategy`` reports the preconditioner it configured, so "what am I + running?" does not require knowing which nine PETSc keys to look up (#484). + + It must still BE the string: comparisons, formatting and serialisation of a + strategy name are unchanged. + """ + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "rep", rotated=False) + + # before any solve: must NOT present the __init__ defaults as resolved + assert s.strategy == "default" + assert f"{s.strategy}" == "default" + assert "not resolved yet" in repr(s.strategy) + + s.solve() + assert s.strategy == "default" # still the plain string + summary = repr(s.strategy) + assert "geometric multigrid (2 levels)" in summary + assert "gmres" in summary and "sor" in summary + assert "not resolved yet" not in summary + + # and the machine-readable form agrees with it + assert s.preconditioner_settings["mg_levels_ksp_type"] == "gmres" + assert s.preconditioner_settings["mg_levels_ksp_max_it"] == "4" + + +def test_strategy_report_names_a_user_override(): + """A key the user took over must be called out, not silently folded into the + summary — otherwise the report reintroduces the ambiguity it exists to remove.""" + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "rov", rotated=False) + s.petsc_options["fieldsplit_velocity_mg_levels_ksp_max_it"] = 6 + s.solve() + summary = repr(s.strategy) + assert "overridden by the user" in summary + assert "mg_levels_ksp_max_it=6" in summary + assert s.preconditioner_settings["mg_levels_ksp_max_it"] == "6" + + def test_rotated_fmg_survives_repeated_newton_increments(): """The rotated path applies its bundle under a per-solve options prefix and then drops the keys again, so the global database stays bounded under From ebb2d3379412475058e996f3f4857bb2c25fa883 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 30 Jul 2026 16:29:33 +1000 Subject: [PATCH 17/17] Adversarial review fixes: the strategy did not reach the custom-P routes, and the strategy value broke pickle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two defects found by an adversarial pass over this PR as a whole, both introduced by its own last two commits. 1. `strategy="fast"` was honoured on the NATIVE route and silently ignored on BOTH custom-P routes. `_configure_pcmg` had no smoother argument, so the standard and rotated paths always got "robust" whatever the user asked for. Measured: native route ('richardson', 3) as asked custom-P standard ('gmres', 4) IGNORED custom-P rotated ('gmres', 4, svd) IGNORED This is precisely the drift class this PR exists to eliminate — a setting honoured on one route and dropped on the two that matter most, since custom-P is mandatory for rotated BCs and adapt() children. It survived because the only test of the strategy axis used a refined mesh, i.e. the native route. `_configure_pcmg` now takes `smoother=` and every install path passes `solver._mg_smoother_variant`. Its docstring says what happens if you forget. 2. The strategy value broke `pickle`, `copy` and `deepcopy` with a TypeError. `str.__reduce_ex__` reconstructs via `cls(value)` with ONE argument, which a two-argument `__new__` cannot accept. `__reduce__` now carries the summary. The claim in 13bd82ce's message that "serialisation is unchanged" was wrong, and UW3 has a serialisation system that would have hit it. Tests for both, each written to fail on the old code: `test_strategy_reaches_every _route` asserts richardson/3 on all three routes (and that the rotated coarse solve is still svd), and `test_strategy_value_pickles_and_copies` round-trips the value. Validation: three routes still agree bar the deliberate rotated svd; 79 passed across the MG, adaptivity and solve-report suites; `level_1 and tier_a` 526 passed with the two pre-existing #470 failures; parallel np2 and np4 17 passed each. Underworld development team with AI support from Claude Code --- .../cython/petsc_generic_snes_solvers.pyx | 10 +++- src/underworld3/utilities/custom_mg.py | 14 +++-- src/underworld3/utilities/rotated_bc.py | 7 ++- tests/test_1021_mg_option_bundle.py | 56 ++++++++++++++++++- 4 files changed, 78 insertions(+), 9 deletions(-) diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 715c78f1..a5c36ee8 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -28,11 +28,19 @@ class _StrategyName(str): form. """ - def __new__(cls, name, summary): + def __new__(cls, name, summary=""): obj = super().__new__(cls, name) obj._summary = summary return obj + def __reduce__(self): + # `str.__reduce_ex__` reconstructs via `cls(value)` with ONE argument, which + # a two-argument `__new__` cannot accept — so without this, pickling, copy + # and deepcopy of a strategy value all raise TypeError. The summary is + # derived state and is carried along rather than recomputed, because the + # solver it came from is not part of the pickle. + return (self.__class__, (str(self), self._summary)) + def __repr__(self): return f"{str.__repr__(self)} — {self._summary}" diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 588b0031..542c8791 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -595,7 +595,7 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") -def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): +def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -603,7 +603,11 @@ def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): the same bundle the native (DMPlex-refinement) route applies — so custom-P and native multigrid cannot be configured differently (#468). ``coarse`` selects the coarse-solve variant: ``"svd"`` on the rotated path, whose - Galerkin-coarsened operator inherits the rigid-rotation null space. ``owned`` + Galerkin-coarsened operator inherits the rigid-rotation null space. ``smoother`` + is the variant ``solver.strategy`` asks for — pass + ``solver._mg_smoother_variant``, or the strategy is honoured on the native route + and silently ignored here, which is the drift this module exists to prevent. + ``owned`` is the solver's record of the option values UW3 has written, so a key the USER set is left alone (see :meth:`multigrid_options.MGSettings.apply`); ``None`` writes unconditionally, which is right for the rotated path's per-solve prefix. @@ -619,7 +623,7 @@ def _configure_pcmg(pc, Ps, coarse="redundant", owned=None): fixed).""" nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" - multigrid_options.geometric_mg_bundle(coarse=coarse).apply( + multigrid_options.geometric_mg_bundle(coarse=coarse, smoother=smoother).apply( PETSc.Options(), prefix, owned=owned) pc.setType("mg") pc.setMGLevels(nlev) @@ -654,6 +658,7 @@ def _install_transfers(solver, Ps, verbose=False): solver.snes.setUp() ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) _configure_pcmg(ksp.getPC(), Ps, + smoother=solver._mg_smoother_variant, owned=solver._managed_pc_options) if verbose: from underworld3 import mpi @@ -715,7 +720,8 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): vel_ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) - _configure_pcmg(vel_pc, Ps, owned=solver._managed_pc_options) + _configure_pcmg(vel_pc, Ps, smoother=solver._mg_smoother_variant, + owned=solver._managed_pc_options) vel_pc.setUp() # 4. re-attach the coupled Stokes nullspace (operator state was touched) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index b83bad4d..948f84d2 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1080,7 +1080,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal # the default redundant/LU coarse solve hits a zero pivot there # (SUBPC_ERROR, outer reason -11). Everything else in the bundle # is identical to the native and standard custom-P routes. - custom_mg._configure_pcmg(vel_pc, custom_Pl, coarse="svd") + custom_mg._configure_pcmg( + vel_pc, custom_Pl, coarse="svd", + smoother=solver._mg_smoother_variant) vel_pc.setUp() # The bundle went into the GLOBAL options DB under the velocity # sub-PC's prefix, which is derived from `pfx` and so is unique to @@ -1089,7 +1091,8 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal vopts = PETSc.Options() vpfx = vel_pc.getOptionsPrefix() or "" for key in multigrid_options.geometric_mg_bundle( - coarse="svd").settings: + coarse="svd", + smoother=solver._mg_smoother_variant).settings: vopts.delValue(vpfx + key) # Constant-pressure nullspace on the Schur COMPLEMENT (enclosed # domains): the IS-built fieldsplit does not propagate the coupled diff --git a/tests/test_1021_mg_option_bundle.py b/tests/test_1021_mg_option_bundle.py index 7468388d..d3122b6a 100644 --- a/tests/test_1021_mg_option_bundle.py +++ b/tests/test_1021_mg_option_bundle.py @@ -128,14 +128,16 @@ def _custom_standard_config(): return _mg_config(s.snes.getKSP().getPC().getFieldSplitSubKSP()[0].getPC()) -def _custom_rotated_config(monkeypatch): +def _custom_rotated_config(monkeypatch, strategy=None): """Route C: custom-P FMG reached through the rotated free-slip path. The rotated KSP is self-contained and is destroyed at the end of the solve, so the configuration is captured from inside the solve — there is no live PC left to interrogate afterwards. """ - s = _stokes(_annulus(RES), "rot", rotated=True) + s = _stokes(_annulus(RES), f"rot{strategy or ''}"[:6], rotated=True) + if strategy is not None: + s.strategy = strategy custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) real_solve = rotated_bc._solve_rotated_iterative @@ -281,6 +283,56 @@ def test_strategy_selects_a_real_smoother_variant(): assert _strategy_mg_config("robust")[:2] == ("gmres", 4) +def test_strategy_reaches_every_route(monkeypatch): + """``strategy="fast"`` must be honoured on ALL THREE routes. + + It was honoured on the native route and silently ignored on both custom-P + routes, because ``_configure_pcmg`` had no smoother argument — the exact drift + this module exists to prevent, and on the routes that matter most (custom-P is + mandatory for rotated BCs and ``adapt()`` children). Only the native route was + covered by a test, which is why it survived. + """ + # _velocity_mg_config returns a tuple, _mg_config a dict — normalise to + # (smoother ksp, smoother max_it, coarse pc) so the three routes are comparable. + native = _strategy_mg_config("fast") + + s = _stokes(_annulus(RES), "fstd", rotated=False) + s.strategy = "fast" + custom_mg.set_custom_fmg(s, [_annulus(2 * RES)], field_id=0) + s.solve() + standard = _velocity_mg_config(s) + + rot = _custom_rotated_config(monkeypatch, strategy="fast") + rotated = (rot["smoother_ksp"], rot["smoother_max_it"], rot["coarse_pc"]) + + for route, config in (("native", native), ("custom-P standard", standard), + ("custom-P rotated", rotated)): + assert config[:2] == ("richardson", 3), ( + f"strategy='fast' not honoured on the {route} route: {config}") + assert rotated[2] == "svd", "the rotated coarse solve must still be svd" + + +def test_strategy_value_pickles_and_copies(): + """The strategy getter returns a ``str`` subclass so the report is additive. It + must therefore behave as a string everywhere, including through ``pickle`` and + ``copy`` — a two-argument ``__new__`` without ``__reduce__`` breaks both, and + UW3 has a serialisation system that would hit it.""" + import copy + import pickle + + mesh = uw.meshing.Annulus(radiusInner=R_IN, radiusOuter=R_OUT, + cellSize=2 * RES, qdegree=3, refinement=1) + s = _stokes(mesh, "pk", rotated=False) + s.solve() + value = s.strategy + + for rebuilt in (pickle.loads(pickle.dumps(value)), copy.copy(value), + copy.deepcopy(value)): + assert rebuilt == "default" + assert str(rebuilt) == "default" + assert f"{value}" == "default" and value + "!" == "default!" + + def test_default_strategy_does_not_change_behaviour(): """The control: ``"default"`` must reproduce the framework default exactly, so filling the strategy axis moves nobody's results."""