The discrete adjoint of a rotated free-slip solve

[lmoresi]

Rotated free-slip never solved K δ = b, so its adjoint was never Kᵀ μ. The forward path rotates into the per-node boundary frame, strikes out the wall-normal rows and inverts Âhat = Q K Qᵀ 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 this 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. SNES_Stokes.adjoint_solve dispatches to it whenever _rotated_freeslip_bcs is non-empty. Fault contact still refuses: its rotated operator carries an additive interface tangent whose transpose is not routed into the adjoint.

The algebra

Everything follows from Q being orthogonal. With μ̂ = Q μ and b̂ = Q b, ⟨μ̂, Âhat δ̂⟩ = ⟨μ, K δ⟩ and ⟨b̂, δ̂⟩ = ⟨b, δ⟩, so the rotated adjoint system is the physical one written in the boundary frame, and μ = Qᵀ μ̂ takes the answer back.

• zeroRowsColumns commutes with transposition (row AND column zeroed, scalar diagonal), so eliminating the constrained rows from Âhatᵀ 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 DOFs: the multiplier's wall-normal component is set, not converged to — measured 4.8e-18.
• The saddle structure survives. UW3 assembles τ − 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ᵀ. Fieldsplit-Schur, the 1/μ pressure mass and the custom-FMG prolongation all apply unchanged.
• The null space serves both sides. Measured, not argued: on a free-slip annulus with a power-law TI tangent the admitted modes give ‖Âhat v‖ = 4.8e-17 / 2.3e-10 against ‖Âhatᵀ v‖ = 4.8e-17 / 2.4e-10.

The test, and what it took to make it discriminate

The contract 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 keeps major symmetry and deleting the transpose entirely still passes — we verified that by deleting it. At the measured 5.7e-2 velocity-block asymmetry, deleting it moves the gradient 1.9% against a 0.2% bound.

The guard on that guard was itself wrong first time round: it read the composite Jacobian, which is structurally non-symmetric for every rheology (2.5e-2 for constant isotropic viscosity) and so proved nothing. It now reads the velocity block — 5.3e-17 isotropic, 6.0e-17 constant TI, 5.7e-2 power-law TI.

Stated in the test rather than left implied: in 2-D with one normal per node Q is EXACTLY symmetric — the SVD frame for a single normal is the Householder reflection [[nx,ny],[ny,−nx]], and ‖Q − Qᵀ‖ measures zero. So swapping Q for Qᵀ is invisible here; it takes a 3-D boundary or a multi-normal corner to pin that axis down, and there is no such test.

Adversarial review before the PR

Six real defects found and fixed in the second commit, each verified first:

Asymmetry guard vacuous	read the composite, not the velocity block
Direct-LU transposed the wrong matrix	zeroRows does not commute with transposition (its transpose is zeroCols); now builds the forward operator and transposes once
_build_rotated_custom_Pl leaks	velocity submatrix + rotated fine prolongation, per call; the forward's is reused
A refusal left the Newton tangent installed	our own test drives a refusal, so a later FORWARD solve would silently run Newton; now try/finally
Stale rotation undetected	the handle guard caught a destroyed Q, not a live and wrong one; the boundary SET is compared
Docstring claims	μ is zero in the wall-normal COMPONENT not the rows; rotated nodes are not Dirichlet-zeroed; only the iterative path projects the null space

The gradient fixture turns out to have 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.

Known, recorded, not fixed here

• The adjoint builds a fresh KSP/PC every call (ctx=None), paying a full PCSetUp on the one path an inversion calls in a loop. An adjoint workspace cache keyed like the forward's would fix that and the prolongation leak together.
• Under mpirun -n 2 the new test file aborts about half the time in _jitextension with "JIT C-source hash differs across MPI ranks". The abort lands in the fixture's FIRST FORWARD solve, before any adjoint runs; test_0021 is stable at np=2 and the same expressions outside pytest never diverge, so this is 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. Filed separately.

Tests

tests/test_0022_rotated_adjoint.py (new, 7 tests). Adjoint suite test_0018–test_0022 and the forward rotated free-slip suite (test_1018, test_1061, test_1065): 75 passed, forward behaviour unchanged.

test_0019's refusal test used rotated free-slip as its stand-in for "refuses"; repointed at fault contact. 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