diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 668947b7..abe8c2e6 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/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). diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 38dc9275..f8f67640 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 @@ -10115,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 ------- @@ -10170,11 +10194,30 @@ 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()) - 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 81ccd457..8ba07077 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1596,6 +1596,191 @@ 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. + + **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 + 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 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: + 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.") + # 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)) + + # 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() + + 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) + + if use_lu: + 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: + # 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 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) + # 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) + + 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 0bbfda38..f2422186 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 31d5e267..d0e7f96c 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 00000000..c68fe86f --- /dev/null +++ b/tests/test_0022_rotated_adjoint.py @@ -0,0 +1,348 @@ +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, 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 + 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). + +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 — issue #752. +""" + +import math + +import numpy as np +import pytest +import sympy + +from petsc4py import PETSc + +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 +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) + # 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] + 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.""" + # 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 +@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"] + # 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"] + + +@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)) + + +@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()