From c8576ab61bccda8971aa40a52b235e200cfd460e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 17 Sep 2026 00:00:02 +1000 Subject: [PATCH 1/3] feat: the discrete adjoint of a rotated free-slip solve Rotated free-slip never solved K delta = b, so the adjoint was never K^T mu. The forward path rotates into the per-node boundary frame, strikes out the wall-normal rows and inverts Ahat = Q K Q^T in its own Krylov loop; the adjoint is the transpose of THAT operator. It was refused for want of a transpose path, which blocked the adjoint for the free-slip the project recommends on any curved or deformed boundary - the annulus and shell benchmarks among them. rotated_bc.solve_rotated_adjoint owns the sequence, as solve_rotated_freeslip owns the forward one. Everything follows from Q being orthogonal: with mu_hat = Q mu and b_hat = Q b the rotated system is the physical one written in the boundary frame, and mu = Q^T mu_hat takes the answer back. zeroRowsColumns commutes with transposition, so eliminating the constrained rows from Ahat^T gives exactly the transpose of the operator the forward solve ran on - the dual of a strong constraint is a strong homogeneous constraint on the same degrees of freedom, and the multiplier's wall-normal component is set rather than converged to (measured 4.8e-18). The block structure and the null space carry over unchanged, so the fieldsplit-Schur setup, the 1/mu pressure mass and the custom-FMG prolongation all apply to the transposed operator. The multiplier is returned modulo the free-slip null space. An enclosed free-slip domain has an undetermined pressure level and, on an annulus, an undetermined rigid rotation; the forward fixes the gauge afterwards and the adjoint has no rest state to fix it against, so that component of b is projected out. Reading Q after a reset would have been a use-after-free, not a wrong answer: the result dict shares Q/Qt with the workspace cache, and _reset_rotated_solver_cache destroys the cache while deliberately keeping the dict (its reaction vector is still wanted). Refused by name instead. Fault contact still refuses - its rotated operator carries an additive interface tangent whose transpose is not routed into the adjoint. The test is the gradient against a central finite difference on an annulus with a power-law transversely isotropic viscosity under the consistent tangent. The anisotropy is measured, not assumed: with a CONSTANT TI viscosity the tangent is major-symmetric and deleting the transpose entirely still passes. At the measured 6.7e-2 asymmetry, deleting it moves the gradient by 1.9% against a 0.2% bound. What the test does not cover is stated in it: in 2-D with one normal per node Q is EXACTLY symmetric (the SVD frame is the Householder reflection [[nx,ny],[ny,-nx]], ||Q - Q^T|| measures zero), so swapping Q for Q^T is invisible here and wants a 3-D case. Also de-duplicates the Picard sentence in the Stokes adjoint verdict, which the base class and the override were each appending. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/rotated-freeslip.md | 75 +++++- .../cython/petsc_generic_snes_solvers.pyx | 68 ++++-- src/underworld3/utilities/rotated_bc.py | 145 +++++++++++ tests/test_0018_adjoint_support_record.py | 21 +- tests/test_0019_adjoint_solve.py | 16 +- tests/test_0022_rotated_adjoint.py | 230 ++++++++++++++++++ 6 files changed, 526 insertions(+), 29 deletions(-) create mode 100644 tests/test_0022_rotated_adjoint.py diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 668947b73..ddb71c4bf 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -345,9 +345,82 @@ vertex error does not. **If you sample only midpoints, pass `mass="consistent"`.** See also #404 (the vertex-integral checkerboard) and #637 (only P1/P2 triangular traces are supported in 3-D at all). +## The adjoint transposes the rotated operator, not the Jacobian + +`solver.adjoint_solve(b)` under rotated free-slip does not transpose `K`. The +forward solve never inverted `K` — it inverted ` = Q K Qᵀ` with the wall-normal +rows struck out — so the adjoint is the transpose of that, taken in the same +frame and rotated back. `rotated_bc.solve_rotated_adjoint` owns the sequence, +and the solver dispatches to it whenever `_rotated_freeslip_bcs` is non-empty. + +Everything follows from `Q` being orthogonal. Writing `μ̂ = Q μ` and `b̂ = Q b`, + +``` +⟨μ̂,  δ̂⟩ = ⟨Qμ, Q K Qᵀ Q δ⟩ = ⟨μ, K δ⟩ ⟨b̂, δ̂⟩ = ⟨Q b, Q δ⟩ = ⟨b, δ⟩ +``` + +so the rotated adjoint system is the physical one written in the boundary +frame, and `μ = Qᵀ μ̂` takes the answer back. Three consequences worth stating +because each one is a place the wiring could be wrong and still look plausible: + +* **The constraint is part of the adjoint.** `zeroRowsColumns` zeroes the row + *and* the column and puts a scalar on the diagonal, so it commutes with + transposition: eliminating the constrained rows from `Âᵀ` gives exactly the + transpose of the operator the forward solve ran on. The multiplier's + wall-normal component is *set* to zero, not iterated towards it — the dual of + a strong constraint is a strong homogeneous constraint on the same degrees of + freedom. +* **The block structure survives.** Transposing `[[A, Bᵀ], [B, 0]]` gives + `[[Aᵀ, Bᵀ], [B, 0]]`, so the fieldsplit-Schur setup, the 1/μ pressure-mass + block and the custom-FMG prolongation all apply to `Âᵀ` unchanged. +* **The null space is shared.** A rigid rotation has zero strain rate, so + `∫C:ε(·):ε(·)` annihilates it read from either side whatever the symmetry of + `C`; the constant-pressure mode couples only through `Bᵀ`, which both + operators carry in the same block. `_rotated_nullspace` therefore serves the + adjoint as it does the forward. + +**The multiplier is returned modulo that null space.** An enclosed free-slip +domain has an undetermined pressure level and, on an annulus or shell, an +undetermined rigid rotation. The forward fixes the gauge after the fact; the +adjoint has no rest state to fix it against, so the component of `b` along +those modes is projected out (`nsp.remove`). A misfit that is itself invariant +under rigid rotation loses nothing to this. One that is not is asking for the +sensitivity of a quantity the forward problem does not determine, and the +projection is what says so — so pin the gauge in the forward problem (a +Dirichlet boundary) rather than reading a rotation-sensitive gradient. + +### What the test cannot see, and why `Qt` is still built + +In 2-D with one normal per node **`Q` is exactly symmetric**: the frame +`numpy.linalg.svd` returns for a single normal is the Householder reflection +`[[nx, ny], [ny, -nx]]`, and `‖Q - Qᵀ‖` measures **zero** on an annulus. So +swapping `Q` for `Qᵀ` anywhere in either the forward or the adjoint path is +numerically invisible there — substituting one for the other in the adjoint's +dual rotation moves the gradient in the eighth digit (measured). That is an +accident of the 2-D single-normal case, not a licence to alias the two: a 3-D +boundary frame is a 3×3 orthogonal matrix whose two tangent rows are not pinned, +and a multi-normal corner block is not a reflection either. `build_rotation` +assembles `Qt` explicitly for that reason, and the adjoint uses `Q` for the dual +and `Qᵀ` for the answer because that is what the duality says, not because a 2-D +test forced it. **A 3-D rotated adjoint test would be the one that pins this +axis down; there isn't one.** + +`J` must be assembled with the **consistent** tangent. A forward that ran +Picard leaves the frozen-viscosity operator on the SNES, which is not `∂R/∂u`; +`adjoint_solve` rebuilds it, and `adjoint_support()` says so in its reason. + +**Fault contact still refuses.** Its rotated operator carries an additive +interface tangent, reassembled at every iterate, and that term's transpose is +not routed into the adjoint. + + ## 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). +independence of the linear datum and the nonlinear Newton datum path), +`tests/test_0022_rotated_adjoint.py` (the adjoint gradient against a central +finite difference on an annulus with a transversely isotropic viscosity — a +curved boundary so the rotation is a real per-node frame, and a tangent with no +major symmetry so a wrong transpose cannot hide behind a symmetric operator). diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 38dc9275c..d58c4115d 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1413,9 +1413,9 @@ class SolverBaseClass(uw_object): transpose for the state and the symbolic derivative of the residual for a parameter. What removes that: - * a rotated constraint — free-slip or fault contact — solves on a - rotated operator inside its own Krylov loop (``rotated_bc.py``), - and there is no transpose path through it; + * a fault contact, whose rotated operator carries an additive + interface tangent reassembled at every iterate — the transpose of + that term is not routed into the adjoint; * an unconverged solve, which is caught after the fact by :meth:`_record_solve_outcome`: a linearisation about a state the solve never reached is not the adjoint of anything. @@ -1425,12 +1425,15 @@ class SolverBaseClass(uw_object): ``tests/test_0018_adjoint_support_record.py``. """ mechanisms = self._constraint_mechanisms() - rotated = mechanisms["rotated_freeslip"] + mechanisms["fault_contact"] - if rotated: + contact = mechanisms["fault_contact"] + if contact: return (False, - f"{len(rotated)} rotated constraint(s): the solve runs on a " - f"rotated operator with its own Krylov loop, and there is " - f"no transpose path through it") + f"{len(contact)} fault contact(s): the rotated operator " + f"carries an additive interface tangent, reassembled at every " + f"iterate, whose transpose is not routed into the adjoint") + # Rotated FREE-SLIP does admit one — the transpose of the rotated + # operator, taken in the boundary frame. The solver that owns that path + # says so in its own verdict; nothing is refused here. if not self.consistent_jacobian and not self._residual_is_linear_in_unknown(): return (True, "implicit residual: Jacobian transpose for the state, symbolic " @@ -10092,12 +10095,23 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): cached = (key, self._residual_is_linear_in_unknown()) self._adjoint_linearity_cache = cached nonlinear = not cached[1] - if not self.consistent_jacobian and nonlinear: - return (True, - why + ". The forward solve used the Picard tangent, so the " - "consistent tangent is assembled for the adjoint (a rebuild " - "of the Jacobian kernel)") - return supported, why + # The base reaches the same conclusion by its own (uncached) linearity + # test, so say it once: this override exists for the CACHE, not for a + # second opinion. Without the guard the sentence appears twice in the + # reason, which reads as two different findings about the same solve. + if not self.consistent_jacobian and nonlinear and "Picard tangent" not in why: + why = (why + ". The forward solve used the Picard tangent, so the " + "consistent tangent is assembled for the adjoint (a rebuild " + "of the Jacobian kernel)") + rotated = list(getattr(self, "_rotated_freeslip_bcs", None) or []) + if rotated: + why = (why + f". {len(rotated)} rotated free-slip constraint(s): the " + "forward solve inverted the ROTATED operator Q J Qᵀ, so the " + "adjoint is the transpose of THAT, taken in the boundary frame " + "and rotated back (rotated_bc.solve_rotated_adjoint). The " + "multiplier comes back modulo the free-slip null space — the " + "pressure level, and a rigid rotation on an annulus or shell") + return True, why def adjoint_solve(self, rhs, target=None): r"""Solve :math:`K^T (\\mu, \\lambda) = b` on the composite (u, p) system. @@ -10108,6 +10122,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): strain-rate- or pressure-dependent viscosity it is the transpose of the consistent tangent, which that construction cannot build. + Under rotated free-slip the operator to transpose is not :math:`K` but + the rotated, constraint-eliminated one the forward solve actually + inverted, :math:`Q K Q^T` — see + :func:`~underworld3.utilities.rotated_bc.solve_rotated_adjoint`. The + multiplier is returned in the physical frame, exactly zero in the + wall-normal component, and modulo the free-slip null space. + Parameters ---------- rhs : numpy.ndarray or petsc4py.PETSc.Vec @@ -10170,10 +10191,21 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): x = gvec.duplicate() x.set(0.0) - ksp = self.snes.getKSP() - ksp.setOperators(J, P) - ksp.solveTranspose(b, x) - reason = int(ksp.getConvergedReason()) + if getattr(self, "_rotated_freeslip_bcs", None): + # The forward solve never ran self.snes: it inverted the ROTATED, + # constraint-eliminated operator in its own Krylov loop. The adjoint + # is the transpose of that operator, not of J, and rotated_bc owns + # the whole rotate / eliminate / solve / rotate-back sequence here as + # it does forward — one place where the constraint is expressed. + from underworld3.utilities.rotated_bc import solve_rotated_adjoint + mu, reason = solve_rotated_adjoint(self, J, b, verbose=False) + mu.copy(x) + mu.destroy() + else: + ksp = self.snes.getKSP() + ksp.setOperators(J, P) + ksp.solveTranspose(b, x) + reason = int(ksp.getConvergedReason()) self._restore_tangent(tangent) if target is not None: diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 81ccd4571..c42717d2c 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1596,6 +1596,151 @@ def rotated_residual(uvec, keep_cartesian=False): "workspace_reused": workspace_reused} + +def solve_rotated_adjoint(solver, J, b, verbose=False): + r"""The discrete adjoint of the last rotated free-slip solve: ``Âᵀ μ̂ = b̂``. + + The forward solve never runs ``snes.solve``. It rotates the Jacobian into the + per-node boundary frame, eliminates the wall-normal rows, and inverts + `` = Q J Qᵀ`` in its own Krylov loop (:func:`solve_rotated_freeslip`). The + adjoint is the transpose of THAT operator, not of ``J``: the constraint is part + of the system the forward solve inverted, so the adjoint carries it too. + + Everything follows from ``Q`` being orthogonal. Writing ``μ̂ = Q μ`` and + ``b̂ = Q b``, ``⟨μ̂,  δ̂⟩ = ⟨μ, J δ⟩`` and ``⟨b̂, δ̂⟩ = ⟨b, δ⟩``, so the rotated + adjoint system IS the physical one expressed in the boundary frame, and + ``μ = Qᵀ μ̂`` takes the answer back. ``zeroRowsColumns`` commutes with + transposition — it zeroes the row AND the column and puts a scalar on the + diagonal — so eliminating the constrained rows from ``Âᵀ`` gives exactly the + transpose of the operator the forward solve ran on. That is the discrete + statement that the dual of a strong constraint is a strong homogeneous + constraint on the same degrees of freedom: ``μ·n̂ = 0`` on the wall. + + ``J`` must already be assembled at the state the forward solve ended in, with + the CONSISTENT tangent. The caller owns that choice: a forward that ran Picard + leaves the frozen-viscosity operator on the SNES, which is not ``∂R/∂u``, and + transposing it would be silently wrong. + + **The null space is projected out of the right-hand side.** A free-slip + enclosed domain has an undetermined pressure level and, on an annulus or + shell, an undetermined rigid rotation. The forward solve fixes the gauge + afterwards; the adjoint cannot, because the multiplier has no rest state to + be fixed against. So the component of ``b`` along those modes is removed + (``nsp.remove``) and ``μ`` is returned modulo them. A misfit that is itself + invariant under a rigid rotation loses nothing to this; one that is not is + asking for the sensitivity of a quantity the forward problem does not + determine, and the projection is what says so. + + Parameters + ---------- + solver + The solver that did the forward solve, still holding + ``_rotated_freeslip_info`` from it. + J : petsc4py.PETSc.Mat + The Jacobian, assembled at the converged state with the consistent + tangent. The pressure-mass Schur block is read from the solver's Pmat, so + assemble both (``snes.computeJacobian(x, J, Jp)``). + b : petsc4py.PETSc.Vec + The dual right-hand side on the composite global vector, in the PHYSICAL + frame — the same thing ``dual_of`` builds for the unrotated path. + + Returns + ------- + (petsc4py.PETSc.Vec, int) + ``μ`` in the physical frame, exactly zero at the constrained wall-normal + rows, and the KSP converged reason (positive means converged). The caller + owns the vector. + """ + info = getattr(solver, "_rotated_freeslip_info", None) + if info is None: + raise RuntimeError( + "solve_rotated_adjoint: no rotated forward solve to take the adjoint " + "of. The rotation Q is built by solve_rotated_freeslip — call solve() " + "first, and take the adjoint about the state it ended in.") + Q, Qt = info.get("Q"), info.get("Qt") + # The result dict SHARES Q/Qt with the cross-solve cache, and + # _reset_rotated_solver_cache destroys that cache while deliberately keeping + # the result dict alive (its reaction vector outlives the workspace, which is + # what boundary_normal_traction needs). The rotation does NOT outlive it, so + # reading Q here after a reset would be a use-after-free rather than a wrong + # answer. Say what happened instead. + if Q is None or Qt is None or Q.handle == 0 or Qt.handle == 0: + raise RuntimeError( + "solve_rotated_adjoint: the rotation from the last rotated solve has " + "been released — the solver was reset or reconfigured since (see " + "_reset_rotated_solver_cache), and Q is shared with the workspace that " + "reset destroys. Re-run solve() and take the adjoint about that state.") + normal_rows = info["normal_rows"] + dm = solver.dm + use_lu = bool(getattr(solver, "_rotated_use_lu", False)) + + # Âᵀ, formed as the forward's own ptap and then transposed rather than from + # Jᵀ.ptap(Qt). The two are equal, and going through the SAME expression means + # a symmetric J (the linear isoviscous case) reproduces the forward operator + # bit-for-bit, which is what makes the symmetric case a usable check on the + # non-symmetric one. + Ahat = J.ptap(Qt) + AhatT = Ahat.transpose(PETSc.Mat()) + Ahat.destroy() + + bhat = AhatT.createVecRight() + Q.mult(b, bhat) + # The constrained rows carry no equation in the forward system, so they carry + # no adjoint equation either; zeroing here is what makes μ·n̂ = 0 exact rather + # than merely converged. + _zero_rows_local(bhat, normal_rows) + + diag_scale = _velocity_diag_scale(AhatT, solver) + AhatT.zeroRowsColumns(normal_rows, diag=diag_scale) + + if use_lu: + # Mirrors the forward direct path, pin included: the adjoint operator has + # the same constant-pressure null mode, so it needs the same gauge fixed + # the same way. + pin = _naive_pressure_pin(dm) + if pin is not None: + AhatT.zeroRows([pin], diag=1.0) + _zero_rows_local(bhat, [pin]) + ksp = PETSc.KSP().create(comm=dm.comm) + ksp.setType("preonly") + pc = ksp.getPC() + pc.setType("lu") + pc.setFactorSolverType("mumps") + ksp.setOperators(AhatT) + muhat = AhatT.createVecRight() + muhat.set(0.0) + ksp.solve(bhat, muhat) + reason = int(ksp.getConvergedReason()) + _warn_if_ksp_diverged(ksp, kind="rotated adjoint direct-LU") + _zero_rows_local(muhat, normal_rows) + ksp.destroy() + else: + # Âᵀ has the block structure of  — transposing [[A, Bᵀ], [B, 0]] gives + # [[Aᵀ, Bᵀ], [B, 0]] — so the fieldsplit-Schur setup, the 1/mu pressure + # mass and the custom-FMG prolongation all apply unchanged. + Mp = _pressure_mass_schur_pmat(solver) + # The rigid-body modes are null modes of Âᵀ as well as Â: a rotation has + # zero strain rate, so it is annihilated by ∫C:ε(·):ε(·) read from either + # side, whatever the symmetry of C. The same is true of the constant + # pressure mode, whose only coupling is through Bᵀ, which both operators + # carry in the same block. + nsp = _rotated_nullspace(solver, Q, normal_rows) + custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) + muhat, reason, ctx = _solve_rotated_iterative( + solver, AhatT, bhat, Q, Qt, normal_rows, verbose=verbose, + custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, ctx=None) + reason = int(reason) + _destroy_rotated_ksp_ctx(ctx) + + mu = dm.createGlobalVec() + Qt.mult(muhat, mu) + + muhat.destroy() + bhat.destroy() + AhatT.destroy() + return mu, reason + + 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 this solver has no hierarchy. Depends only on Q and diff --git a/tests/test_0018_adjoint_support_record.py b/tests/test_0018_adjoint_support_record.py index 0bbfda380..f24221869 100644 --- a/tests/test_0018_adjoint_support_record.py +++ b/tests/test_0018_adjoint_support_record.py @@ -9,8 +9,11 @@ * an implicit step is a residual: Jacobian transpose for the state, symbolic derivative for a parameter — supported; - * a rotated constraint solves on a rotated operator in its own Krylov loop - with no transpose path — refused; + * a rotated free-slip constraint solves on a rotated operator in its own + Krylov loop, so the adjoint is the transpose of THAT operator — supported, + and the verdict has to say which operator it means; + * a fault contact adds an interface tangent to that rotated operator whose + transpose is not routed into the adjoint — refused; * a solve that did not converge is linearised about a state it never reached — refused, after the fact; * a semi-Lagrangian history's interpolation at the departure points is not @@ -92,10 +95,11 @@ def test_a_plain_implicit_solve_is_supported(mesh): assert "Jacobian transpose" in verdict["reason"] -def test_a_rotated_constraint_refuses_and_says_why(mesh): - """The rotated solve runs its own Krylov loop on a rotated operator; there - is no transpose path through it. The transcript must say so rather than - let the solve pass as an ordinary residual.""" +def test_a_rotated_constraint_is_supported_on_the_rotated_operator(mesh): + """The rotated solve runs its own Krylov loop on a rotated operator, and the + adjoint is the transpose of that one rather than of K. A verdict that said + only "supported" would be describing the wrong solve, so the reason has to + name the rotation AND the null space the multiplier is returned modulo.""" uw, model = _fresh_model() stokes, _ = _stokes(uw, mesh, "rot") stokes.add_dirichlet_bc((0.0, 0.0), "Bottom") @@ -108,8 +112,9 @@ def test_a_rotated_constraint_refuses_and_says_why(mesh): with model.step(0.1): stokes.solve() verdict = _adjoint_of(model, "solve") - assert verdict["supported"] is False - assert "rotated" in verdict["reason"] + assert verdict["supported"] is True + assert "rotated free-slip" in verdict["reason"] + assert "null space" in verdict["reason"] def test_an_unconverged_solve_is_refused_after_the_fact(mesh): diff --git a/tests/test_0019_adjoint_solve.py b/tests/test_0019_adjoint_solve.py index 31d5e267d..d0e7f96c0 100644 --- a/tests/test_0019_adjoint_solve.py +++ b/tests/test_0019_adjoint_solve.py @@ -125,6 +125,15 @@ def J_at(value): def test_a_refusing_solve_raises_with_its_reason(): + """A refused verdict must stop ``adjoint_solve`` and carry its own reason, + rather than let it return a number nobody can trace. + + Driven by a fault contact, the mechanism that still refuses — its rotated + operator carries an additive interface tangent whose transpose is not routed + into the adjoint. (Rotated FREE-SLIP used to stand here; it is supported now, + and ``test_0022_rotated_adjoint.py`` checks its gradient.) The fault list is + populated directly rather than by splitting a mesh: what is under test is the + refusal dispatch, and only the length of that list reaches the verdict.""" uw, model = _fresh() mesh = _mesh(uw) V = uw.discretisation.MeshVariable("V_ref", mesh, 2, degree=2) @@ -132,8 +141,11 @@ def test_a_refusing_solve_raises_with_its_reason(): stokes = uw.systems.Stokes(mesh, velocityField=V, pressureField=P) stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 - stokes.add_rotated_freeslip_bc(0.0, "Top") - with pytest.raises(RuntimeError, match="rotated"): + stokes._fault_contact_faults.append("a-fault") + supported, reason = stokes.adjoint_support() + assert supported is False + assert "fault contact" in reason + with pytest.raises(RuntimeError, match="fault contact"): stokes.adjoint_solve(np.zeros(1)) diff --git a/tests/test_0022_rotated_adjoint.py b/tests/test_0022_rotated_adjoint.py new file mode 100644 index 000000000..fee8d025b --- /dev/null +++ b/tests/test_0022_rotated_adjoint.py @@ -0,0 +1,230 @@ +r"""The adjoint of a solve that ran on the ROTATED operator. + +Rotated free-slip does not solve :math:`K \delta = b`. It rotates into the +per-node boundary frame, strikes out the wall-normal rows, and inverts +:math:`Q K Q^T` in its own Krylov loop. The adjoint is the transpose of THAT +operator, and the multiplier comes back with no wall-normal component — the +dual of a strong constraint is a strong homogeneous constraint on the same +degrees of freedom. + +The check that counts is the gradient against a central finite difference. Two +things in the setup are there to stop a wrong answer hiding, and both are +guarded rather than asserted in a comment: + + * an ANNULUS, so ``Q`` is a genuine per-node frame rather than the signed + permutation an axis-aligned wall would give. Rotating the dual the wrong + way (``Qᵀ b`` where ``Q b`` is meant) then breaks the gradient. + * a POWER-LAW transversely isotropic viscosity under the CONSISTENT (Newton) + tangent, so the tangent has no major symmetry and ``K ≠ Kᵀ``. A symmetric + operator cannot tell a transpose from itself, and every isotropic Stokes + case is symmetric — as is the frozen TI tangent, whose ``C`` keeps both + minor and major symmetry. ``test_the_tangent_is_not_symmetric`` measures + the asymmetry rather than trusting the setup to have produced it: without + that guard, omitting the transpose entirely still passes (measured). +""" + +import math + +import numpy as np +import pytest +import sympy + +from petsc4py import PETSc + +import underworld3 as uw +from underworld3.adjoint import misfit_duals + + +R_I, R_O = 0.5, 1.0 +ETA_1 = 0.2 # the parameter, at the point the gradient is taken +FD_STEP = 1.0e-4 +# The asymmetry the transpose is supposed to matter for. A tangent this close to +# symmetric could not distinguish K from Kᵀ, so the fixture measures it and the +# gradient test is only meaningful above this floor. +MIN_ASYMMETRY = 1.0e-3 + + +@pytest.fixture(scope="module") +def rotated_gradient(): + """One rotated free-slip solve, its adjoint gradient, and the central finite + difference to check it against. Module-scoped: the finite difference re-solves + at two shifted parameters, so the state is not reusable afterwards and every + contract has to be read off the record this builds.""" + mesh = uw.meshing.Annulus(radiusInner=R_I, radiusOuter=R_O, + cellSize=0.15, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + unit_r = sympy.Matrix([[x / r, y / r]]) + th = sympy.atan2(y, x) + + v = uw.discretisation.MeshVariable("v_rot_adj", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p_rot_adj", mesh, 1, degree=1) + v_obs = uw.discretisation.MeshVariable("v_obs_rot_adj", mesh, 2, degree=2) + + eta_1 = uw.expression(r"\eta_1", ETA_1, "weak-plane viscosity ratio") + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + # Power law, so the CONSISTENT tangent picks up the dependence of the + # anisotropic C on the strain rate — that term is what destroys the major + # symmetry the frozen TI tangent keeps. + edot = mesh.vector.strain_tensor(v.sym) + eII = sympy.sqrt(sympy.Rational(1, 2) * (edot[0, 0] ** 2 + edot[1, 1] ** 2) + + edot[0, 1] ** 2) + eta_0 = (sympy.Float(0.01) + eII) ** sympy.Rational(-1, 3) + + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_0 + stokes.constitutive_model.Parameters.shear_viscosity_1 = eta_1 * eta_0 + # radial director: the weak planes are the tangential ones, and the director + # turns with the boundary, so the anisotropy and the rotation frame are both + # functions of position rather than constants that could coincidentally align. + stokes.constitutive_model.Parameters.director = unit_r + stokes.bodyforce = 1.0e2 * sympy.cos(3 * th) * (r - R_I) / (R_O - R_I) * unit_r + stokes.add_dirichlet_bc((0.0, 0.0), "Lower") + stokes.add_rotated_freeslip_bc(0.0, "Upper") + stokes.consistent_jacobian = True + stokes.tolerance = 1.0e-11 + + def solve_at(value): + eta_1.sym = sympy.Float(value) + stokes.solve(zero_init_guess=True) + + # An observation set from a DIFFERENT parameter, so the misfit is not + # stationary at the point the gradient is taken and the finite difference has + # something to measure. + solve_at(0.05) + v_obs.array[...] = np.asarray(v.array) + + misfit = sympy.Rational(1, 2) * ((v.sym[0] - v_obs.sym[0]) ** 2 + + (v.sym[1] - v_obs.sym[1]) ** 2) + + def J_at(value): + solve_at(value) + return float(uw.maths.Integral(mesh, misfit).evaluate()) + + J_at(ETA_1) + supported, reason = stokes.adjoint_support() + + # Is the tangent this gradient is taken through actually non-symmetric? + # Measured on the assembled Jacobian at the converged state, the same one + # the adjoint transposes. Without this the whole file can pass with the + # transpose deleted. + K = stokes.snes.getJacobian()[0] + Kt = K.transpose(PETSc.Mat()) + Kt.axpy(-1.0, K, structure=PETSc.Mat.Structure.DIFFERENT_NONZERO_PATTERN) + asymmetry = Kt.norm(PETSc.NormType.FROBENIUS) / K.norm(PETSc.NormType.FROBENIUS) + Kt.destroy() + + dual = misfit_duals(misfit, [v])[v] + dual.array[...] = -np.asarray(dual.array) + mu = uw.discretisation.MeshVariable("mu_rot_adj", mesh, 2, degree=2) + lam = uw.discretisation.MeshVariable("lam_rot_adj", mesh, 1, degree=1) + mu_global, reason_ksp = stokes.adjoint_solve((dual, None), target=(mu, lam)) + adjoint = stokes.sensitivity(mu, eta_1) + + # The wall-normal leak of the MULTIPLIER, measured in the discrete frame the + # constraint is actually written in — Q's per-node rows, not an analytic + # normal. Those two differ by the facet/true-normal discrepancy, which is a + # property of the mesh and has nothing to say about the adjoint. + info = stokes._rotated_freeslip_info + vec = stokes.dm.createGlobalVec() + vec.array[:] = mu_global + rotated = vec.duplicate() + info["Q"].mult(vec, rotated) + lo, hi = rotated.getOwnershipRange() + owned = [g - lo for g in info["normal_rows"] if lo <= g < hi] + leak = float(np.abs(np.asarray(rotated.array)[owned]).max()) if owned else 0.0 + leak = uw.mpi.comm.allreduce(leak, op=uw.MPI.MAX) + vec.destroy() + rotated.destroy() + + fd = (J_at(ETA_1 + FD_STEP) - J_at(ETA_1 - FD_STEP)) / (2 * FD_STEP) + eta_1.sym = sympy.Float(ETA_1) + + return {"supported": supported, "reason": reason, "ksp_reason": reason_ksp, + "adjoint": adjoint, "fd": fd, "leak": leak, + "asymmetry": float(asymmetry)} + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_tangent_is_not_symmetric(rotated_gradient): + """The guard on the guard. A symmetric K satisfies Kᵀ = K, so the gradient + test below would pass with the transpose deleted — as it does when the TI + viscosity is constant. Assert the asymmetry is there before believing the + gradient says anything about a transpose. Measured 6.7e-2 for this setup; + deleting the transpose then moves the gradient by 1.9%, against a 0.2% + bound.""" + assert rotated_gradient["asymmetry"] > MIN_ASYMMETRY, ( + f"tangent is (near-)symmetric at {rotated_gradient['asymmetry']:.2e} — " + f"the gradient test cannot detect a missing transpose; strengthen the " + f"anisotropy or the power-law exponent in the fixture") + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_gradient_through_a_rotated_freeslip_solve_matches_finite_differences( + rotated_gradient): + """The contract. A transpose taken on K rather than on Q K Qᵀ, or a dual + rotated the wrong way, changes this number; the finite difference does not + care how the constraint was imposed.""" + adjoint, fd = rotated_gradient["adjoint"], rotated_gradient["fd"] + assert adjoint != 0.0 + assert abs(fd / adjoint - 1) < 2.0e-3, (fd, adjoint) + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_multiplier_has_no_wall_normal_component(rotated_gradient): + """Exact, not converged. ``zeroRowsColumns`` decouples the constrained rows + of the adjoint operator exactly as it does the forward one, so the + multiplier's wall-normal component is set rather than iterated towards — + a Krylov tolerance must not appear in this number.""" + assert rotated_gradient["leak"] < 1.0e-12, rotated_gradient["leak"] + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_adjoint_ksp_converged(rotated_gradient): + assert rotated_gradient["ksp_reason"] > 0, rotated_gradient["ksp_reason"] + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_verdict_says_the_operator_is_the_rotated_one(rotated_gradient): + """A supported verdict that did not mention the rotation would be describing + a different solve — the one that transposes K.""" + assert rotated_gradient["supported"] is True + assert "rotated free-slip" in rotated_gradient["reason"] + assert "null space" in rotated_gradient["reason"] + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_a_released_rotation_is_refused_rather_than_read(): + """``_rotated_freeslip_info`` outlives the workspace cache it SHARES ``Q`` + with — deliberately, because the reaction vector in it is still wanted after + a reset. The rotation is not still wanted: reading it there is a + use-after-free, not a wrong answer. It has to be refused by name.""" + mesh = uw.meshing.Annulus(radiusInner=R_I, radiusOuter=R_O, + cellSize=0.3, qdegree=3) + v = uw.discretisation.MeshVariable("v_rel", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p_rel", mesh, 1, degree=1) + dual = uw.discretisation.MeshVariable("d_rel", mesh, 2, degree=2) + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1 + stokes.bodyforce = sympy.Matrix([[0.0, -1.0]]) + stokes.add_dirichlet_bc((0.0, 0.0), "Lower") + stokes.add_rotated_freeslip_bc(0.0, "Upper") + stokes.tolerance = 1.0e-9 + stokes.solve() + dual.array[...] = 1.0 + + _, reason = stokes.adjoint_solve((dual, None)) + assert reason > 0, reason + + stokes._reset_rotated_solver_cache() + with pytest.raises(RuntimeError, match="released"): + stokes.adjoint_solve((dual, None)) From bcaef1559db951bc345e4a9c0de32aa974a39e2e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 17 Sep 2026 01:00:48 +1000 Subject: [PATCH 2/3] review: fix the rotated adjoint's pin transpose, prolongation leak and vacuous guard Adversarial review of the previous commit. Six real defects, each verified before fixing. The asymmetry guard was VACUOUS. It measured the composite Jacobian, but UW3 assembles the velocity flux as tau - p.I against +div u, so the operator is [[A, -B^T], [B, 0]] and the composite is structurally non-symmetric for every rheology: measured 2.5e-2 for CONSTANT ISOTROPIC viscosity, well above the 1e-3 floor the guard was asserting. The guard now reads the velocity block, where the numbers are 5.3e-17 isotropic, 6.0e-17 constant TI (the frozen tangent keeps major symmetry) and 5.7e-2 for the power-law TI the fixture uses. Deleting the transpose still moves the gradient 1.9% against a 0.2% bound, so the conclusion stands - but it stood on the negative control alone, not on the guard. The direct-LU path transposed the wrong matrix. zeroRowsColumns commutes with transposition; zeroRows does NOT (its transpose is zeroCols), and the forward pins the pressure with zeroRows to keep the B^T coupling. Transposing after the pin rather than reproducing the pin in transposed form makes this exactly M^T by construction, for both paths. _build_rotated_custom_Pl leaks the velocity submatrix and the rotated fine prolongation on every call - nothing owns them, and the cache only dereferences the list. An inversion loop paid that per adjoint. The forward's prolongation depends only on Q and the hierarchy, so it is reused where it exists. A refusal used to leave the solver switched to the Newton tangent: _consistent_tangent_for_adjoint sets it and the restore was a bare statement after the solve. The new path refuses by design, and one of our own tests drives a refusal - so a later FORWARD solve would silently have run Newton on a solver configured for Picard. Now try/finally. Also: refuse when the rotated boundary SET has changed since the forward solve (the handle guard only catches a DESTROYED Q, not a live and wrong one, and the is_setup gate is bypassed once the adjoint kernel is installed); put _rotated_velocity_null_modes back, since an adjoint must not change what the next forward solve does; correct the docstring claims that mu is zero at the constrained ROWS (it is the wall-normal COMPONENT), that Dirichlet zeroing covers rotated nodes (it does not - they keep a tangential multiplier), and that this function projects the null space (only the iterative path does). Two measurements that changed what the prose says. The null space DOES serve the transposed operator - on a free-slip annulus with a power-law TI tangent the admitted modes give ||A v|| = 4.8e-17 / 2.3e-10 and ||A^T v|| = 4.8e-17 / 2.4e-10 - and there is now a test that measures it instead of arguing it. And the gradient fixture has NO null space at all: it pins the inner boundary, so _rotated_nullspace returns None. Pinning one boundary is the configuration to prefer when a sensitivity is the point, and the docs now say so. The leak bound moves 1e-12 -> 1e-14: at a forward tolerance of 1e-11 the old bound was one decade below what a deleted _zero_rows_local would leave, which is not a baseline. The measured value is 4.8e-18. Known and recorded, not fixed here: the adjoint builds a fresh KSP/PC per call, and under mpirun -n 2 the file aborts about half the time in _jitextension with "JIT C-source hash differs across MPI ranks". That abort lands in the fixture's first FORWARD solve, before any adjoint; test_0021 is stable at np=2 and the same expressions outside pytest never diverge, so it is a latent JIT non-determinism rather than anything this branch introduced. CI's parallel pass collects only tests/parallel/, so it is not in the gate. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/rotated-freeslip.md | 37 ++++- .../cython/petsc_generic_snes_solvers.pyx | 45 ++++--- src/underworld3/utilities/rotated_bc.py | 90 +++++++++---- tests/test_0022_rotated_adjoint.py | 126 +++++++++++++++++- 4 files changed, 246 insertions(+), 52 deletions(-) diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index ddb71c4bf..213ed3fa7 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -370,16 +370,29 @@ because each one is a place the wiring could be wrong and still look plausible: wall-normal component is *set* to zero, not iterated towards it — the dual of a strong constraint is a strong homogeneous constraint on the same degrees of freedom. -* **The block structure survives.** Transposing `[[A, Bᵀ], [B, 0]]` gives - `[[Aᵀ, Bᵀ], [B, 0]]`, so the fieldsplit-Schur setup, the 1/μ pressure-mass - block and the custom-FMG prolongation all apply to `Âᵀ` unchanged. +* **The block structure survives.** UW3 assembles the velocity flux as `τ − pI` + against `+div u`, so the operator is `[[A, −Bᵀ], [B, 0]]` and its transpose is + `[[Aᵀ, Bᵀ], [−B, 0]]` — the off-diagonal **signs swap**, the blocks do not + move, and the two swapped signs cancel in `B A⁻ᵀ Bᵀ` so the Schur complement + keeps its sign as well as its sparsity. The fieldsplit-Schur setup, the 1/μ + pressure-mass block and the custom-FMG prolongation all apply to `Âᵀ` + unchanged. (That `−Bᵀ` is also why a symmetry check on the COMPOSITE matrix + says nothing: it reads ~2.5e-2 for constant isotropic viscosity. The velocity + block is the one to measure — 5.3e-17 there, against 5.7e-2 for a power-law + TI tangent.) * **The null space is shared.** A rigid rotation has zero strain rate, so `∫C:ε(·):ε(·)` annihilates it read from either side whatever the symmetry of `C`; the constant-pressure mode couples only through `Bᵀ`, which both operators carry in the same block. `_rotated_nullspace` therefore serves the - adjoint as it does the forward. - -**The multiplier is returned modulo that null space.** An enclosed free-slip + adjoint as it does the forward — measured, not argued: on a free-slip annulus + with a power-law TI tangent the admitted modes give `‖ v‖` = 4.8e-17 / 2.3e-10 + and `‖Âᵀ v‖` = 4.8e-17 / 2.4e-10. + +**Where there is a null space, the multiplier is returned modulo it.** Pinning +one boundary with an essential condition removes it — `_rotated_nullspace` then +returns `None`, and the gradient is unambiguous. That is the configuration to +prefer when a sensitivity is the point, and it is what +`test_0022_rotated_adjoint`'s gradient fixture uses. An enclosed free-slip domain has an undetermined pressure level and, on an annulus or shell, an undetermined rigid rotation. The forward fixes the gauge after the fact; the adjoint has no rest state to fix it against, so the component of `b` along @@ -405,6 +418,18 @@ and `Qᵀ` for the answer because that is what the duality says, not because a 2 test forced it. **A 3-D rotated adjoint test would be the one that pins this axis down; there isn't one.** +### Cost, and what is not cached + +The adjoint builds a fresh KSP/PC every call (`ctx=None`), so each one pays a +full fieldsplit plus GAMG/FMG `PCSetUp` — the cost the #417 cross-solve cache +exists to avoid, on the one path an inversion calls in a loop. The prolongation +IS reused from that cache when the forward built one (it depends only on `Q` and +the hierarchy), which matters for more than speed: `_build_rotated_custom_Pl` +leaks the velocity submatrix and the rotated fine prolongation on every call — +nothing owns them, and `_destroy_rotated_linear_cache` only dereferences the +list. An adjoint workspace cache keyed the same way as the forward's would fix +both; it is not built. + `J` must be assembled with the **consistent** tangent. A forward that ran Picard leaves the frozen-viscosity operator on the SNES, which is not `∂R/∂u`; `adjoint_solve` rebuilds it, and `adjoint_support()` says so in its reason. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index d58c4115d..f8f67640f 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -10136,7 +10136,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): one from a velocity-space expression. target : (MeshVariable, MeshVariable), optional ``(u_adj, p_adj)`` on the velocity and pressure spaces, to receive - the multipliers as fields. Constrained nodes are set to zero. + the multipliers as fields. Essential (Dirichlet) nodes are set to + zero. A rotated free-slip node is NOT one of those: it is + constrained in one component only, and comes back carrying its + tangential multiplier. Returns ------- @@ -10191,22 +10194,30 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): x = gvec.duplicate() x.set(0.0) - if getattr(self, "_rotated_freeslip_bcs", None): - # The forward solve never ran self.snes: it inverted the ROTATED, - # constraint-eliminated operator in its own Krylov loop. The adjoint - # is the transpose of that operator, not of J, and rotated_bc owns - # the whole rotate / eliminate / solve / rotate-back sequence here as - # it does forward — one place where the constraint is expressed. - from underworld3.utilities.rotated_bc import solve_rotated_adjoint - mu, reason = solve_rotated_adjoint(self, J, b, verbose=False) - mu.copy(x) - mu.destroy() - else: - ksp = self.snes.getKSP() - ksp.setOperators(J, P) - ksp.solveTranspose(b, x) - reason = int(ksp.getConvergedReason()) - self._restore_tangent(tangent) + try: + if getattr(self, "_rotated_freeslip_bcs", None): + # The forward solve never ran self.snes: it inverted the ROTATED, + # constraint-eliminated operator in its own Krylov loop. The + # adjoint is the transpose of that operator, not of J, and + # rotated_bc owns the whole rotate / eliminate / solve / + # rotate-back sequence here as it does forward — one place where + # the constraint is expressed. + from underworld3.utilities.rotated_bc import solve_rotated_adjoint + mu, reason = solve_rotated_adjoint(self, J, b, verbose=False) + mu.copy(x) + mu.destroy() + else: + ksp = self.snes.getKSP() + ksp.setOperators(J, P) + ksp.solveTranspose(b, x) + reason = int(ksp.getConvergedReason()) + finally: + # _consistent_tangent_for_adjoint SWITCHED the solver to the Newton + # tangent. The rotated path refuses by design — a released rotation, + # a changed boundary set — and a refusal that left the switch in place + # would make the next FORWARD solve run Newton on a solver configured + # for Picard, silently and far from here. + self._restore_tangent(tangent) if target is not None: u_adj, p_adj = target diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index c42717d2c..8ba070770 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1621,7 +1621,12 @@ def solve_rotated_adjoint(solver, J, b, verbose=False): leaves the frozen-viscosity operator on the SNES, which is not ``∂R/∂u``, and transposing it would be silently wrong. - **The null space is projected out of the right-hand side.** A free-slip + **Where there IS a null space, it is projected out of the right-hand side** + (iterative path; the direct-LU path fixes the pressure gauge with the same + naive pin the forward solve uses and builds none). A boundary with an + essential condition removes it altogether — `_rotated_nullspace` then returns + None and the gradient is unambiguous, which is the configuration to prefer + when a sensitivity is what you are after. A free-slip enclosed domain has an undetermined pressure level and, on an annulus or shell, an undetermined rigid rotation. The forward solve fixes the gauge afterwards; the adjoint cannot, because the multiplier has no rest state to @@ -1647,9 +1652,10 @@ def solve_rotated_adjoint(solver, J, b, verbose=False): Returns ------- (petsc4py.PETSc.Vec, int) - ``μ`` in the physical frame, exactly zero at the constrained wall-normal - rows, and the KSP converged reason (positive means converged). The caller - owns the vector. + ``μ`` in the physical frame and the KSP converged reason (positive means + converged). The caller owns the vector. What is exactly zero is the + wall-normal COMPONENT — ``(Q μ)`` at ``normal_rows`` — not the rows of + ``μ`` itself, which carry the tangential multiplier once rotated back. """ info = getattr(solver, "_rotated_freeslip_info", None) if info is None: @@ -1670,16 +1676,35 @@ def solve_rotated_adjoint(solver, J, b, verbose=False): "been released — the solver was reset or reconfigured since (see " "_reset_rotated_solver_cache), and Q is shared with the workspace that " "reset destroys. Re-run solve() and take the adjoint about that state.") + # Q covers the boundaries of the LAST solve. add_rotated_freeslip_bc appends + # and sets is_setup False, but the gate that would catch that in adjoint_solve + # is bypassed once the adjoint kernel is installed — so a boundary registered + # between two adjoints would silently run on a rotation that does not cover + # it. Compare the sets rather than trust the ordering. + registered = list(getattr(solver, "_rotated_freeslip_bcs", None) or []) + if list(info.get("boundaries") or []) != registered: + raise RuntimeError( + f"solve_rotated_adjoint: the rotated boundaries changed since the " + f"forward solve — the rotation covers {info.get('boundaries')}, the " + f"solver now carries {registered}. Re-run solve() so Q and the " + f"constrained rows describe the problem being differentiated.") normal_rows = info["normal_rows"] dm = solver.dm use_lu = bool(getattr(solver, "_rotated_use_lu", False)) - # Âᵀ, formed as the forward's own ptap and then transposed rather than from - # Jᵀ.ptap(Qt). The two are equal, and going through the SAME expression means - # a symmetric J (the linear isoviscous case) reproduces the forward operator - # bit-for-bit, which is what makes the symmetric case a usable check on the - # non-symmetric one. + # Build the operator the forward solve ACTUALLY INVERTED — ptap, then the + # constraint elimination, then (direct path only) the gauge pin — and + # transpose that whole thing ONCE. Transposing at the end rather than + # reproducing each step in transposed form is what makes this exactly Mᵀ: + # ``zeroRowsColumns`` commutes with transposition, but ``zeroRows`` does NOT + # (its transpose is zeroCols), so a pinned operator assembled the other way + # round would be the transpose of a different matrix. Ahat = J.ptap(Qt) + diag_scale = _velocity_diag_scale(Ahat, solver) + Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) + pin = _naive_pressure_pin(dm) if use_lu else None + if pin is not None: + Ahat.zeroRows([pin], diag=1.0) AhatT = Ahat.transpose(PETSc.Mat()) Ahat.destroy() @@ -1690,17 +1715,7 @@ def solve_rotated_adjoint(solver, J, b, verbose=False): # than merely converged. _zero_rows_local(bhat, normal_rows) - diag_scale = _velocity_diag_scale(AhatT, solver) - AhatT.zeroRowsColumns(normal_rows, diag=diag_scale) - if use_lu: - # Mirrors the forward direct path, pin included: the adjoint operator has - # the same constant-pressure null mode, so it needs the same gauge fixed - # the same way. - pin = _naive_pressure_pin(dm) - if pin is not None: - AhatT.zeroRows([pin], diag=1.0) - _zero_rows_local(bhat, [pin]) ksp = PETSc.KSP().create(comm=dm.comm) ksp.setType("preonly") pc = ksp.getPC() @@ -1715,22 +1730,47 @@ def solve_rotated_adjoint(solver, J, b, verbose=False): _zero_rows_local(muhat, normal_rows) ksp.destroy() else: - # Âᵀ has the block structure of  — transposing [[A, Bᵀ], [B, 0]] gives - # [[Aᵀ, Bᵀ], [B, 0]] — so the fieldsplit-Schur setup, the 1/mu pressure - # mass and the custom-FMG prolongation all apply unchanged. + # The saddle structure survives transposition. UW3 assembles the velocity + # flux as τ − p·I against +div u, so the operator is [[A, −Bᵀ], [B, 0]] + # and its transpose is [[Aᵀ, Bᵀ], [−B, 0]] — the OFF-DIAGONAL SIGNS SWAP, + # the blocks do not move. The Schur complement keeps both its sparsity and + # its sign (the two swapped signs cancel in B A⁻ᵀ Bᵀ), so the 1/mu pressure + # mass, the fieldsplit and the custom-FMG prolongation all still apply. Mp = _pressure_mass_schur_pmat(solver) # The rigid-body modes are null modes of Âᵀ as well as Â: a rotation has # zero strain rate, so it is annihilated by ∫C:ε(·):ε(·) read from either # side, whatever the symmetry of C. The same is true of the constant - # pressure mode, whose only coupling is through Bᵀ, which both operators - # carry in the same block. + # pressure mode, whose only coupling is through the off-diagonal block, + # which transposition moves but does not remove. `_mode_satisfies_constraints` + # verifies admission against J, not Jᵀ; the argument above is why that is + # the same test, and `test_0022_rotated_adjoint` measures it on a tangent + # with no major symmetry rather than leaving it as an argument. + # + # _rotated_nullspace records its mode count on the SOLVER, where the + # forward path reads it to choose a coarse solve. An adjoint must not + # change what the next forward solve does, so put it back below. + null_modes_before = getattr(solver, "_rotated_velocity_null_modes", None) nsp = _rotated_nullspace(solver, Q, normal_rows) - custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) + # The prolongation depends only on Q and the mesh hierarchy, so the + # forward's is reusable — and MUST be reused where it exists: + # `_build_rotated_custom_Pl` leaks the velocity submatrix and the rotated + # fine prolongation on every call (nothing owns them; the cache only + # DEREFERENCES the list), which an inversion loop would pay per adjoint. + cache = getattr(solver, "_rotated_linear_cache", None) + if cache is not None and cache.get("Q") is Q and cache.get("custom_Pl"): + custom_Pl = cache["custom_Pl"] + else: + custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) muhat, reason, ctx = _solve_rotated_iterative( solver, AhatT, bhat, Q, Qt, normal_rows, verbose=verbose, custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, ctx=None) reason = int(reason) _destroy_rotated_ksp_ctx(ctx) + if null_modes_before is None: + if hasattr(solver, "_rotated_velocity_null_modes"): + del solver._rotated_velocity_null_modes + else: + solver._rotated_velocity_null_modes = null_modes_before mu = dm.createGlobalVec() Qt.mult(muhat, mu) diff --git a/tests/test_0022_rotated_adjoint.py b/tests/test_0022_rotated_adjoint.py index fee8d025b..c80c8054c 100644 --- a/tests/test_0022_rotated_adjoint.py +++ b/tests/test_0022_rotated_adjoint.py @@ -12,8 +12,8 @@ guarded rather than asserted in a comment: * an ANNULUS, so ``Q`` is a genuine per-node frame rather than the signed - permutation an axis-aligned wall would give. Rotating the dual the wrong - way (``Qᵀ b`` where ``Q b`` is meant) then breaks the gradient. + permutation an axis-aligned wall would give, and so the geometry is the one + the free-slip adjoint is wanted for. * a POWER-LAW transversely isotropic viscosity under the CONSISTENT (Newton) tangent, so the tangent has no major symmetry and ``K ≠ Kᵀ``. A symmetric operator cannot tell a transpose from itself, and every isotropic Stokes @@ -21,6 +21,27 @@ minor and major symmetry. ``test_the_tangent_is_not_symmetric`` measures the asymmetry rather than trusting the setup to have produced it: without that guard, omitting the transpose entirely still passes (measured). + +What this file does NOT check, so nobody reads it as covered: in 2-D with one +normal per node ``Q`` is EXACTLY symmetric — the frame ``numpy.linalg.svd`` +returns for a single normal is the Householder reflection +``[[nx, ny], [ny, -nx]]``, and ``‖Q - Qᵀ‖`` measures zero on this mesh. So +rotating the dual the wrong way (``Qᵀ b`` where ``Q b`` is meant) is invisible +here; substituting it changes the gradient in the eighth digit. The direction is +right on the mathematics — ``Q`` is orthogonal, so ``μ̂ = Q μ`` and ``μ = Qᵀ μ̂`` +— but it takes a 3-D boundary or a multi-normal corner, where the frame is no +longer a reflection, to make a test say so. + +Serial. Under ``mpirun -n 2`` this file aborts about half the time in +``_jitextension`` with "JIT C-source hash differs across MPI ranks" — the +deliberate hard error for non-deterministic ``generate_c_source``. It is NOT an +adjoint failure: the abort lands in the fixture's FIRST FORWARD solve, before +any adjoint runs. Nor is it this file's expression — the pre-existing TI adjoint +test (``test_0021``) is stable at np=2, and building the fixture's expression +variants outside pytest never diverges (3/3). Fixing ``PYTHONHASHSEED`` does not +settle it. The CI parallel pass collects only ``tests/parallel/test_*.py``, so +this is not in the gate; the cause is a latent non-determinism in JIT source +generation, filed separately. """ import math @@ -33,6 +54,8 @@ import underworld3 as uw from underworld3.adjoint import misfit_duals +from underworld3.utilities.rotated_bc import (_rotated_nullspace, + _velocity_diag_scale) R_I, R_O = 0.5, 1.0 @@ -112,7 +135,22 @@ def J_at(value): K = stokes.snes.getJacobian()[0] Kt = K.transpose(PETSc.Mat()) Kt.axpy(-1.0, K, structure=PETSc.Mat.Structure.DIFFERENT_NONZERO_PATTERN) - asymmetry = Kt.norm(PETSc.NormType.FROBENIUS) / K.norm(PETSc.NormType.FROBENIUS) + # The VELOCITY BLOCK, not the composite. UW3 assembles the velocity flux as + # τ − p·I against +div u, so the operator is [[A, −Bᵀ], [B, 0]] and the + # composite is structurally non-symmetric for EVERY rheology — measured + # 2.5e-2 for constant isotropic viscosity, which would sail past any floor + # set here and prove nothing. It is the A block whose symmetry decides + # whether a transpose is detectable, and there it is 5.3e-17 isotropic, + # 6.0e-17 for a CONSTANT TI viscosity (the frozen tangent keeps major + # symmetry) and 5.7e-2 for the power-law TI below. + vel_is = stokes._subdict["velocity"][0] + A = K.createSubMatrix(vel_is, vel_is) + dA = Kt.createSubMatrix(vel_is, vel_is) + asymmetry = (dA.norm(PETSc.NormType.FROBENIUS) + / A.norm(PETSc.NormType.FROBENIUS)) + A.destroy() + dA.destroy() + Kt.destroy() dual = misfit_duals(misfit, [v])[v] @@ -180,7 +218,12 @@ def test_the_multiplier_has_no_wall_normal_component(rotated_gradient): of the adjoint operator exactly as it does the forward one, so the multiplier's wall-normal component is set rather than iterated towards — a Krylov tolerance must not appear in this number.""" - assert rotated_gradient["leak"] < 1.0e-12, rotated_gradient["leak"] + # 1e-14, not 1e-12: the forward tolerance is 1e-11, so a bound of 1e-12 is + # only one decade below a number that a DELETED `_zero_rows_local` would + # leave at ~tolerance x ‖b̂‖ — which clears 1e-12 whenever that scale is + # below 0.1, and the fixture does not measure it. What is genuinely left + # here is round-off on the Q(Qᵀμ̂) round trip (measured 4.8e-18). + assert rotated_gradient["leak"] < 1.0e-14, rotated_gradient["leak"] @pytest.mark.level_2 @@ -196,6 +239,10 @@ def test_the_verdict_says_the_operator_is_the_rotated_one(rotated_gradient): a different solve — the one that transposes K.""" assert rotated_gradient["supported"] is True assert "rotated free-slip" in rotated_gradient["reason"] + # The verdict warns about the null space unconditionally, because it is a + # property of the problem the caller poses, not of this solve. This fixture + # pins the inner boundary and HAS no null space (`_rotated_nullspace` + # returns None) — which is exactly why its gradient is unambiguous. assert "null space" in rotated_gradient["reason"] @@ -228,3 +275,74 @@ def test_a_released_rotation_is_refused_rather_than_read(): stokes._reset_rotated_solver_cache() with pytest.raises(RuntimeError, match="released"): stokes.adjoint_solve((dual, None)) + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_the_null_space_serves_the_transposed_operator_too(): + """``_rotated_nullspace`` admits a mode by measuring ``‖·w‖``, and the + adjoint attaches the result to ``Âᵀ`` as BOTH its null space and its + transpose null space. Sound only if the modes are null from the left as well. + + The argument is that a rigid rotation has zero strain rate, so ``∫C:ε:ε`` + annihilates it read from either side whatever the symmetry of ``C``, and the + constant-pressure mode couples only through an off-diagonal block that + transposition moves but does not remove. On a tangent with no major symmetry + that deserves measuring. + + Free slip on BOTH boundaries, because that is what admits a rigid rotation: + the gradient fixture pins the inner boundary and has no null space at all.""" + mesh = uw.meshing.Annulus(radiusInner=R_I, radiusOuter=R_O, + cellSize=0.25, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + unit_r = sympy.Matrix([[x / r, y / r]]) + th = sympy.atan2(y, x) + v = uw.discretisation.MeshVariable("v_nsp", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p_nsp", mesh, 1, degree=1) + + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + edot = mesh.vector.strain_tensor(v.sym) + eII = sympy.sqrt(sympy.Rational(1, 2) * (edot[0, 0] ** 2 + edot[1, 1] ** 2) + + edot[0, 1] ** 2) + eta_0 = (sympy.Float(0.01) + eII) ** sympy.Rational(-1, 3) + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_0 + stokes.constitutive_model.Parameters.shear_viscosity_1 = ETA_1 * eta_0 + stokes.constitutive_model.Parameters.director = unit_r + stokes.bodyforce = 1.0e2 * sympy.cos(3 * th) * (r - R_I) / (R_O - R_I) * unit_r + stokes.add_rotated_freeslip_bc(0.0, "Lower") + stokes.add_rotated_freeslip_bc(0.0, "Upper") + stokes.petsc_use_pressure_nullspace = True + stokes.consistent_jacobian = True + stokes.tolerance = 1.0e-10 + stokes.solve(zero_init_guess=True) + + info = stokes._rotated_freeslip_info + Q, Qt, rows = info["Q"], info["Qt"], info["normal_rows"] + nsp = _rotated_nullspace(stokes, Q, rows) + assert nsp is not None, "no null space here — the test would prove nothing" + assert getattr(stokes, "_rotated_velocity_null_modes", 0) >= 1, ( + "no RIGID ROTATION was admitted; only the pressure mode is being " + "measured, and that one is symmetric for a trivial reason") + + K = stokes.snes.getJacobian()[0] + Ahat = K.ptap(Qt) + Ahat.zeroRowsColumns(rows, diag=_velocity_diag_scale(Ahat, stokes)) + AhatT = Ahat.transpose(PETSc.Mat()) + try: + for w in nsp.getVecs(): + out = w.duplicate() + Ahat.mult(w, out) + forward = out.norm() / w.norm() + AhatT.mult(w, out) + adjoint = out.norm() / w.norm() + out.destroy() + # Absolute, and then the identity: the same mode read from the two + # sides must give the SAME residual, not merely a small one. + assert adjoint < 1.0e-8, (forward, adjoint) + assert abs(adjoint - forward) <= 0.01 * max(forward, 1.0e-16), ( + forward, adjoint) + finally: + Ahat.destroy() + AhatT.destroy() From c1015157b937a3e21e944bce7a56a34c17eaf81e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 17 Sep 2026 07:29:59 +1000 Subject: [PATCH 3/3] docs: point the rotated-adjoint parallel note at issue #752 Underworld development team with AI support from Claude Code --- docs/developer/subsystems/rotated-freeslip.md | 4 ++++ tests/test_0022_rotated_adjoint.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 213ed3fa7..abe8c2e6d 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -430,6 +430,10 @@ nothing owns them, and `_destroy_rotated_linear_cache` only dereferences the list. An adjoint workspace cache keyed the same way as the forward's would fix both; it is not built. +`tests/test_0022_rotated_adjoint.py` is serial: at np>1 it trips the JIT +rank-divergence guard about half the time (#752), in the forward solve rather +than the adjoint. + `J` must be assembled with the **consistent** tangent. A forward that ran Picard leaves the frozen-viscosity operator on the SNES, which is not `∂R/∂u`; `adjoint_solve` rebuilds it, and `adjoint_support()` says so in its reason. diff --git a/tests/test_0022_rotated_adjoint.py b/tests/test_0022_rotated_adjoint.py index c80c8054c..c68fe86f9 100644 --- a/tests/test_0022_rotated_adjoint.py +++ b/tests/test_0022_rotated_adjoint.py @@ -41,7 +41,7 @@ variants outside pytest never diverges (3/3). Fixing ``PYTHONHASHSEED`` does not settle it. The CI parallel pass collects only ``tests/parallel/test_*.py``, so this is not in the gate; the cause is a latent non-determinism in JIT source -generation, filed separately. +generation — issue #752. """ import math