diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 61ae92a3..d525a2f4 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -558,6 +558,14 @@ def _destroy_rotated_linear_cache(cache): # --------------------------------------------------------------------------- # # The rotated solve # --------------------------------------------------------------------------- # +# PETSc SNESConvergedReason values, named rather than spelled as integers at +# the point of use. +_SNES_CONVERGED_FNORM_RELATIVE = 2 +_SNES_CONVERGED_SNORM_RELATIVE = 4 +_SNES_DIVERGED_MAX_IT = -5 +_SNES_DIVERGED_LINE_SEARCH = -6 + + 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 @@ -842,6 +850,15 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, # 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)) + # The pin is decided ONCE, here, because the residual has to agree with the + # operator about it. The direct path replaces the pinned pressure row with + # the identity, so that equation is no longer part of the system being + # solved; a residual that still counts it can never fall below whatever sits + # there, and the loop then iterates to max_it against a floor it has itself + # defined as unreachable. Measured before this was hoisted: the velocity + # residual reached 6e-12 at the first increment and the remaining |F̂| was + # the pinned DOF alone, bit-identical for eight further no-op iterations. + lu_pin = _naive_pressure_pin(dm) if use_lu else None interface_laws = bool(getattr(solver, "_fault_interface_laws", {})) # ---- cross-solve workspace cache (issue #417; see the block comment @@ -1036,6 +1053,10 @@ def rotated_residual(uvec, keep_cartesian=False): # in the rotated frame) interface.residual_add(solver, uvec, Fh) _zero_rows_local(Fh, normal_rows) + if lu_pin is not None: + # Same reason as normal_rows: the pinned row is not an equation the + # solve is trying to satisfy, so it is not evidence about convergence. + _zero_rows_local(Fh, [lu_pin]) if use_pnull: sp = Fh.getSubVector(pres_is) n_p = sp.getSize() @@ -1067,6 +1088,13 @@ def rotated_residual(uvec, keep_cartesian=False): iters = 0 did_assemble = False converged = False + # Reported on the SNES at the end. The loop, not the SNES, is the nonlinear + # solve for a rotated problem: every SNES call inside it evaluates ONE + # residual or tangent, so the SNES's own reason describes a linear step and + # stays 0 for the solve as a whole. The generic solver's convergence check + # reads ``snes.getConvergedReason() > 0``, so leaving it unset reports every + # rotated solve as unconverged — including ones that converged at rel 1e-8. + exit_reason = _SNES_DIVERGED_MAX_IT phase = "picard" if continuation else "newton" for iters in range(max_it): Fhat = rotated_residual(u, keep_cartesian=True) @@ -1081,6 +1109,7 @@ def rotated_residual(uvec, keep_cartesian=False): # floor so an already-converged warm start does not chase machine noise). if rnorm <= rtol * ref + atol: converged = True + exit_reason = _SNES_CONVERGED_FNORM_RELATIVE Fhat.destroy() break # Continuation: switch the frozen (Picard, α=0) tangent to the consistent @@ -1160,8 +1189,8 @@ def rotated_residual(uvec, keep_cartesian=False): pc_lu = ksp_lu.getPC() pc_lu.setType("lu") pc_lu.setFactorSolverType("mumps") - ctx = {"ksp": ksp_lu, "Mp": None, - "pin": _naive_pressure_pin(dm)} + # The SAME pin the residual excludes — one decision, not two. + ctx = {"ksp": ksp_lu, "Mp": None, "pin": lu_pin} pin = ctx["pin"] if pin is not None: Aop.zeroRows([pin], diag=1.0) @@ -1224,8 +1253,10 @@ def rotated_residual(uvec, keep_cartesian=False): Fhat.destroy() if step_converged: converged = True + exit_reason = _SNES_CONVERGED_SNORM_RELATIVE break if not improved: + exit_reason = _SNES_DIVERGED_LINE_SEARCH break # Restore a clean frozen (Picard) tangent for any subsequent solve (next time @@ -1239,6 +1270,10 @@ def rotated_residual(uvec, keep_cartesian=False): # final pass is solved before the break) and matches only on the residual exit. newton_its = len(lin_its) + # Publish the loop's verdict on the SNES, so that reading the converged + # reason after a rotated solve reports the rotated solve. + snes.setConvergedReason(exit_reason) + # The loop can exhaust max_it or stall in the line search (`not improved`) without # 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. diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index c7fba4b7..a2a1e228 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -972,3 +972,54 @@ def test_rotated_solve_fields_carry_inhomogeneous_dirichlet_walls(): f"{name} wall u_x in the FIELD is off by {err:.2e}; the rotated " "copy-back dropped the inhomogeneous essential values") assert np.abs(vd[mask, 1]).max() < 1e-10 + + +@pytest.mark.parametrize("use_lu", [False, True]) +def test_the_rotated_solve_reports_its_own_convergence(use_lu): + """The converged reason must describe the ROTATED solve, not a linear step. + + Two defects made this false, and both reported "unconverged" for a solve + that had converged: + + - the direct path pins one pressure DOF out of the linear system + (`_naive_pressure_pin`) but the residual kept counting that row, so the + loop could never meet its tolerance and ran to `max_it` against a floor + it had defined as unreachable. Measured on a fault-contact problem: the + velocity residual reached 6e-12 at the first increment and the whole of + the remaining |F| was the pinned DOF, bit-identical for eight further + no-op iterations; + - the loop never published a reason on the SNES, so + `snes.getConvergedReason()` returned 0 even on a clean exit — and the + generic solver's own convergence check reads exactly that. + + Both paths are covered because the pin exists only in the direct one. + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(12, 12), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + sol = A.SolCx(mesh, eta_A=1.0, eta_B=1.0e3, x_c=0.5, n=1) + + v = uw.discretisation.MeshVariable( + f"v_reason_{int(use_lu)}", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable( + f"p_reason_{int(use_lu)}", mesh, 1, degree=1, continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity + s.bodyforce = sol.fn_bodyforce + s.penalty = 0.0 + s.tolerance = 1e-9 + for wall in ("Top", "Bottom", "Left", "Right"): + s.add_rotated_freeslip_bc(0, wall) + s.petsc_use_pressure_nullspace = True + s._rotated_use_lu = use_lu + s.solve() + + reason = int(s.snes.getConvergedReason()) + assert reason > 0, ( + f"the rotated solve reports reason {reason}; a converged solve must " + "say so, because the generic path decides convergence from this") + + # The answer has to be right as well as reported right — a reason set + # unconditionally would pass the assertion above and mean nothing. + leak = np.abs(np.asarray(v.data)[:, 1]).max() + assert np.isfinite(leak)