Repository navigation
The discrete adjoint of a rotated free-slip solve (re-opens #751, which never reached development) - #846
Open
lmoresi wants to merge 6 commits into
Open
The discrete adjoint of a rotated free-slip solve (re-opens #751, which never reached development)#846lmoresi wants to merge 6 commits into
lmoresi wants to merge 6 commits into
Conversation
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
…d 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
Underworld development team with AI support from Claude Code
replacement feature/discrete-adjoint has moved a long way since this branch forked. It now assembles the adjoint operator from TRANSPOSED KERNELS where it can (_adjoint_by_kernels / _install_adjoint_kernels), falling back to transposing the assembled matrix only where a boundary Jacobian is in play. The single conflict is where that choice is made, and resolving it by keeping either side would have been wrong. The rotated dispatch is now a third case ahead of the other two, and it takes the MATRIX route by construction: solve_rotated_adjoint builds Q K Qᵀ from the assembled FORWARD Jacobian and transposes the result, so it needs K. Under the kernel route J already IS K^T, and rotating that would rotate the wrong matrix and produce a plausible, wrong gradient. _adjoint_by_kernels() therefore returns False whenever a rotated constraint is present, which is what makes the J reaching solve_rotated_adjoint the forward operator. Said at both ends, because the two facts are only correct together. Teardown is now in `finally` across all three routes, and picks the right one: _uninstall_adjoint_kernels for the kernel route, _restore_tangent otherwise. The rotated path refuses by design — a released rotation, a changed boundary set — and a refusal that left the solver on the Newton tangent would make the next FORWARD solve run Newton on a solver configured for Picard. The merge had also left a duplicated _restore_tangent call; removed. test_0022_rotated_adjoint 7 passed; adjoint suite test_0018-0026 45 passed; deprecated-pattern scan clean. Underworld development team with AI support from Claude Code
The discrete adjoint of a rotated free-slip solve
9 tasks done
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Re-opens #751 against
development. #751 is marked MERGED but its work never reacheddevelopment, so the rotated free-slip adjoint is not in the tree.What happened
feature/discrete-adjoint→developmentf414da63feature/adjoint-rotated-bc→feature/discrete-adjoint69ccc320Thirty-nine seconds apart, in that order. #751's base was the branch #744 had just merged away, so its merge landed on a branch nothing downstream reads:
Five commits were stranded there. GitHub reports #751 as MERGED, which is why this was invisible — the PR looks done and the code is absent. Nothing was lost, and nothing had to be cherry-picked: with #744 in,
feature/discrete-adjointminusdevelopmentis now exactly #751's own delta.Contents
Unchanged from #751 — the same five commits, same diff:
The adjoint transposes the ROTATED operator rather than K, which is what unblocks the annulus. #751 carries the review that got it there, including the pin transpose, the prolongation leak and the vacuous guard.
Verification
Rebuilt from clean (
rm -rf build/, a.pyxchange) at this head and run againstdevelopmentf414da63:test_0022_rotated_adjoint,test_0018_adjoint_support_record,test_0019_adjoint_solve-k adjoint)test_1018_rotated_freeslip,test_1060_nitsche_freeslip,test_1063_constrained_traction,test_0023_field_realness_and_derivatives,test_1066_stokes_jacobian_layoutThe last group is there because this branch rewrites
rotated_bc.pyand the shared solver.pyx.Worth noting for next time
A stacked PR whose base merges first does not follow it. GitHub retargets an open PR when its base branch is deleted, but #751 was merged, not retargeted, and the merge silently went to the orphaned base. When a stack lands, check the child's merge commit is an ancestor of
developmentrather than trusting the MERGED badge.Underworld development team with AI support from Claude Code
🤖 Generated with Claude Code
https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv