The discrete adjoint of a rotated free-slip solve (re-opens #751, which never reached development)

docs/developer/subsystems/rotated-freeslip.md

@@ -345,9 +345,111 @@ 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.** 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 — 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
+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.**
+
+### 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.
+
+`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.
+
+**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).

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -1985,9 +1985,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.
@@ -1997,12 +1997,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 "
@@ -2568,7 +2571,12 @@ class SolverBaseClass(uw_object):
         kernels rather than obtained by transposing the assembled matrix.
         The kernel route covers the volume terms; a boundary Jacobian (a
         natural or Nitsche condition with a tangent) is not swapped yet, so
-        such a solver takes the matrix route."""
+        such a solver takes the matrix route. So does a rotated constraint: its
+        adjoint is the transpose of the ROTATED operator Q K Qᵀ, which
+        rotated_bc builds from the assembled FORWARD Jacobian — handing it K^T
+        would rotate the wrong matrix."""
+        if getattr(self, "_rotated_freeslip_bcs", None):
+            return False
         return not (getattr(self, "natural_bcs", None) or [])
 
     def _install_adjoint_kernels(self):
@@ -11668,12 +11676,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.
@@ -11684,14 +11703,24 @@ 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
             The dual on the composite global vector; :meth:`dual_of` builds
             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
         -------
@@ -11748,16 +11777,41 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         x = gvec.duplicate()
         x.set(0.0)
 
-        if by_kernels:
-            # J IS K^T: assembled from the transposed kernels. Solve forwards.
-            ksp = self.snes.getKSP()
-            ksp.setOperators(J, P)
-            ksp.solve(b, x)
-            reason = int(ksp.getConvergedReason())
-            self._uninstall_adjoint_kernels(tangent)
-        else:
-            reason = _solve_transposed(self.snes.getKSP(), J, P, b, x)
-            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 K, and
+                # rotated_bc owns the whole rotate / eliminate / solve /
+                # rotate-back sequence here as it does forward.
+                #
+                # This is the MATRIX route by construction: solve_rotated_adjoint
+                # forms Q K Qᵀ from the assembled forward Jacobian and transposes
+                # the result, so it needs K and not Kᵀ. _adjoint_by_kernels()
+                # returns False whenever a rotated constraint is present, which
+                # is what guarantees J is the forward operator here.
+                from underworld3.utilities.rotated_bc import solve_rotated_adjoint
+                mu, reason = solve_rotated_adjoint(self, J, b, verbose=False)
+                mu.copy(x)
+                mu.destroy()
+            elif by_kernels:
+                # J IS K^T: assembled from the transposed kernels. Solve forwards.
+                ksp = self.snes.getKSP()
+                ksp.setOperators(J, P)
+                ksp.solve(b, x)
+                reason = int(ksp.getConvergedReason())
+            else:
+                reason = _solve_transposed(self.snes.getKSP(), J, P, b, x)
+        finally:
+            # Teardown belongs in `finally`: the rotated path refuses by design —
+            # a released rotation, a changed boundary set — and a refusal that
+            # left the solver switched to the Newton tangent would make the next
+            # FORWARD solve run Newton on a solver configured for Picard,
+            # silently and far from here.
+            if by_kernels:
+                self._uninstall_adjoint_kernels(tangent)
+            else:
+                self._restore_tangent(tangent)
 
         if target is not None:
             u_adj, p_adj = target