From 9ee07f9f577f95605fe4cba6577627ff5ad12d31 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 27 Jul 2026 22:58:50 +1000 Subject: [PATCH 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 8/8] 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})"