diff --git a/CLAUDE.md b/CLAUDE.md index 0f4a1272..b5f471dd 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -326,7 +326,7 @@ The PETSc-based solvers are carefully optimized and validated. **NO CHANGES with ## Boundary Conditions: Free-slip -**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is the only implemented datum) +**Prefer rotated strong free-slip** (`solver.add_rotated_freeslip_bc(conds, boundary, normal=None)`, value-first: `conds=0` is free-slip; a non-zero scalar/expression `conds` prescribes the wall-normal datum `u·n̂ = ũ_n` strongly) to impose `v·n̂ = 0`: - Enforces zero wall-normal flow to **machine precision** (Nitsche / penalty leak ~1e-3). diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 39caa892..2a8f28f2 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,55 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026) + +**Rotated strong free-slip now takes a prescribed wall-normal velocity datum +(`add_rotated_freeslip_bc(conds, boundary, ...)` with non-zero `conds`) through +the full nonlinear Newton machinery, and the separate linear and nonlinear +solve paths have been unified into one** (#438; #403 items 2 and 4). + +The rotated constraint `u·n̂ = ũ_n` is the primitive behind both the held free-slip +lid and the free-surface material-boundary condition — they differ only in the +constraint right-hand side. Previously the non-zero datum existed only on a +linear one-shot path, so a power-law or anisotropic rheology silently fell back +to a weak penalty (which leaks worst exactly where anisotropy makes it matter). +Now every rotated solve runs a single Newton/Picard loop in which accepted +iterates carry the datum exactly; a cold start imposes it through the first +increment's affine lift at the rest-state tangent (snapping a zero state onto a +datum creates a boundary strain state a shear-thinning tangent cannot recover +from — measured, not assumed); a linear model simply converges after that first +increment, so the up-front nonlinearity probe is gone from the dispatch and +every linear rotated solve saves two Jacobian assemblies. + +Three latent solver defects were found and fixed by making the loop report +honestly along the way: + +- Branching on rank-local datum bookkeeping desynchronised the ranks' + collective sequences (an np>1 deadlock class); the datum-activity decision is + now a collective PETSc reduction. +- A rigid-rotation mode pinned by an essential condition on *another* boundary + could still enter the solver null space (only the rotated rows were checked), + silently projecting an irreducible component out of every increment — Newton + converged superlinearly and then floored, far above tolerance. Candidate + modes are now verified as null vectors of the assembled operator, which + catches any form of pinning. +- A tiny Newton step was reported as convergence even when the residual was + still large (a stiff tangent also produces tiny steps); the step-norm exit is + now verified against the problem's rest-state residual scale, which is also + the convergence reference for warm starts (rtol relative to a good warm + start's own small initial residual demands ever-more absolute accuracy). + +The free-surface manager's `consistent_constraint="strong"` therefore works for +nonlinear rheologies: the penalty fallback is removed, and the consistent solve +warm-starts from the free solve's converged fields. Acceptance: power-law +annulus free-surface convection holds the material boundary at the strong-datum +level (5e-3, versus 4e-2 for the penalty) with net surface flux 1e-4; a +transversely isotropic fault-bearing smoke test converges through the same +path. Direct LU per Newton increment remains available as a serial, +preconditioner-free diagnostic (`solver._rotated_use_lu`). + +New subsystem documentation: `subsystems/rotated-freeslip.md`. + ### Local Interpolation That Reproduces Linear Fields (July 2026) **The local scattered-point interpolator now has a linear-reproduction diff --git a/docs/developer/index.md b/docs/developer/index.md index 28749336..cc663caf 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -36,6 +36,7 @@ same topic are reference or historical material subordinate to the governing doc | Coding style & API conventions | [UW3 Style Charter](UW3_STYLE_CHARTER.md) (detailed reference: [Style Guide](UW3_Style_and_Patterns_Guide.md)) | | Data access | [subsystems/data-access.md](subsystems/data-access.md) (internals reference: [NDArray System](UW3_Developers_NDArrays.md)) | | Local scattered-point interpolation | [subsystems/interpolation.md](subsystems/interpolation.md) | +| Rotated free-slip & wall-normal datum | [subsystems/rotated-freeslip.md](subsystems/rotated-freeslip.md) | | Units | [design/UNITS_SIMPLIFIED_DESIGN_2025-11.md](design/UNITS_SIMPLIFIED_DESIGN_2025-11.md) | | Testing tiers | [TESTING-RELIABILITY-SYSTEM.md](TESTING-RELIABILITY-SYSTEM.md) | | Branching & releases | [guides/branching-strategy.md](guides/branching-strategy.md) | @@ -186,6 +187,7 @@ subsystems/mesh-shape-relaxation subsystems/discretisation subsystems/solvers subsystems/boundary-stress-and-projection-postprocessing +subsystems/rotated-freeslip subsystems/petsc-jacobian-layout subsystems/constitutive-models subsystems/constitutive-models-theory diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md new file mode 100644 index 00000000..7e30897e --- /dev/null +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -0,0 +1,130 @@ +# Rotated strong free-slip and the prescribed-normal datum + +The rotated free-slip boundary condition imposes + +$$\mathbf{u}\cdot\hat{\mathbf{n}} = \tilde u_n$$ + +strongly — as an exact per-node constraint, not a weak (Nitsche/penalty) term — +by rotating each boundary node's velocity components into a (normal, +tangential) frame and constraining the rotated normal component. `conds = 0` +(the default datum) is free-slip; a non-zero scalar `conds` prescribes the +wall-normal velocity, which is what the free-surface manager uses to keep the +surface a material boundary. + +```python +stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) # free-slip +stokes.add_rotated_freeslip_bc(h_dot.sym[0], "Upper", normal=nhat) # u·n̂ = field +``` + +`normal=None` uses the geometric facet normal; a sympy `1×dim` matrix in +`mesh.X` supplies an analytic normal (exact `X/|X|` on curved boundaries — the +preferred choice there); a constant array is also accepted. The datum must be a +*scalar* (a number, an expression of `mesh.X`, or a scalar field read); on an +enclosed boundary it must be discretely flux-free for incompressibility. A +corner or 3D-edge node shared between rotated boundaries has no single normal +and stays at the free-slip pinning (the datum is ignored there). + +Why strong rather than Nitsche/penalty: the constraint holds to machine +precision (a penalty leaks ~1e-3, and the leak grows exactly where anisotropy +makes the boundary condition matter), it is correct on curved/tilted/deformed +boundaries, and the constraint **reaction is the boundary normal traction** +σ_nn (`solver.boundary_normal_traction`, `solver.dynamic_topography`) — the +quantity the free-surface machinery consumes. Reserve Nitsche for conditions +that must morph in time (Dirichlet→Neumann ramps). + +**Implementation**: `src/underworld3/utilities/rotated_bc.py`; registration and +dispatch in `petsc_generic_snes_solvers.pyx` (`add_rotated_freeslip_bc`). + +## One solve path + +There is a single driver, `solve_rotated_freeslip`: a manual outer +Newton/Picard loop that rotates the residual and tangent every iteration +(`F̂ = Q F`, `Ĵ = Q J Qᵀ`), imposes the constraint on the rotated normal rows, +and solves each increment with a self-contained fieldsplit-Schur KSP (geometric +FMG on the custom prolongation when a hierarchy is registered, else GAMG; the +native 1/μ pressure mass as the Schur preconditioner). There is no separate +linear path and no up-front nonlinearity probe: a linear model converges after +its first increment and the loop self-terminates. + +The manual loop exists because the rotated operator `Q A Qᵀ` carries no DM +field information, so PETSc's DM-coupled fieldsplit cannot precondition it; +the increment is solved by an IS-built fieldsplit instead, driven from the +loop. The loop honours `zero_init_guess`, `picard`, and the solver's +`consistent_jacobian` tangent policy (frozen / Newton / continuation). + +### How the datum passes through Newton + +The constraint is affine, and the design keeps every **accepted iterate +feasible** (`u·n̂ = ũ_n` exactly, via an affine snap in the rotated frame), so +each Newton increment satisfies the *homogeneous* constraint `n̂·δ = 0` — the +datum never touches the tangent, the increment right-hand side, the FMG +prolongation, or the null-space handling. + +The one deliberate exception is a **cold start with a non-zero datum**: the +first increment carries the datum jump through the affine lift +(`zeroRowsColumns(rows, diag, x̂, b̂)`) at the rest-state tangent, and is +accepted without a line search. Snapping the zero state onto the datum instead +manufactures an extreme boundary-strip strain state whose shear-thinning +tangent produces a first step orders of magnitude too large — no line-searched +step descends (measured). A warm start is assumed smooth and takes the exact +snap directly. + +### Convergence semantics + +The loop converges on `‖F̂‖ ≤ rtol·ref + atol` with +`ref = max(‖F̂(u₀)‖, ‖F̂(0)‖)`: the **rest-state residual is the intrinsic +forcing scale**, so a good warm start (whose own initial residual is small) is +not punished with an ever-stricter absolute target. A step-norm exit (tiny +Newton step) is *verified* against the same reference before it may report +convergence — a stiff tangent also produces tiny steps far from the solution, +and an unverified tiny step goes back through the line search (a slow crawl +continues; a genuine stall ends the loop with an explicit warning). An +unconverged exit always warns; the fields hold the last iterate. + +### Null-space handling + +Rigid-rotation candidates (one mode on a closed circle, three on a spherical +shell) are admitted to the increment null space and the post-solve gauge +removal only if they pass **two** tests: tangential to every rotated +constraint row (`Q·m ≈ 0` there), **and** a null vector of the assembled +operator (`‖J·m‖ ≈ 0` against the velocity diagonal scale). The second test is +what catches pinning the first cannot see — an essential condition on another +boundary leaves the mode tangential to the rotated wall while the eliminated +operator is *not* null on it, and admitting it projects an irreducible +component out of every increment (the residual then floors far above +tolerance instead of converging). The null space is therefore built after the +first Jacobian assembly. + +### Parallel notes + +- Every branch the loop takes is decided **collectively**. The datum-activity + flag rides a PETSc `Vec` norm of the lift vector (the datum bookkeeping + itself is rank-local — only owners of datum boundary nodes hold entries — + and branching on it per-rank desynchronises the collective sequences: the + np>1 deadlock class). +- All row surgery on vectors uses ownership-relative indices + (`_zero_rows_local` / `_set_rows_local`); indexing a local slice with global + rows is the np>1 crash class. +- Direct LU per increment (`solver._rotated_use_lu = True`) is a **serial** + preconditioner-free diagnostic (the pressure-gauge pin is a naive per-rank + scan, marked `TODO(BUG)`); use it to separate "the operator/constraint is + wrong" from "the preconditioner is struggling". + +## Result and reaction + +The solve fills the solver's fields and stores a result dict +(`solver._rotated_freeslip_info`) with the rotation, the constrained rows, the +per-increment KSP iteration counts, the convergence verdict, and the +**reaction** — the converged Cartesian residual `F(u)`, which for a linear +residual equals `A·u − b` exactly. `boundary_normal_traction` projects the +nodal reaction onto the boundary normal (corner-correct) and de-smears it with +the shared boundary-mass machinery in `utilities/boundary_flux.py`; +`dynamic_topography_field` writes `h = −(σ_nn − σ̄_nn)/(Δρ g)` onto a surface +field for the free-surface integrator. + +## 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). diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 79f960fb..4be2d14c 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -2802,6 +2802,15 @@ def flux(self): class TransverseIsotropicFlowModel(ViscousFlowModel): + # TODO(BUG): the CONSISTENT (Newton) tangent of this model is inconsistent + # with its residual when the anisotropy is active (eta_1 != eta_0): PETSc + # -snes_test_jacobian at bounded curvature reads ||J-Jfd||/||J|| ~ 2e-3-5e-3 + # in developed flow where the isotropic control is FD-limited (~4e-6) — the + # dC/d(eps_II) director-coupled terms are missing or wrong, so TI "Newton" + # runs at Picard pace. Native-path evidence (no rotated machinery); the + # isotropic limit through this class is clean at rest but shares the defect + # once eta_1 differs. See issue #457 (incl. the checker-regularisation + # methodology: at 1e-12 the FD reference itself is invalid for power-law). r""" Transversely isotropic (anisotropic) viscous flow model. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 96703b51..5ff3ce6b 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -1176,15 +1176,13 @@ class SolverBaseClass(uw_object): ``ksp.solve()`` on the rotated operator (``utilities/rotated_bc.py``) and never touches ``self.snes`` — so the generic reader would see stale SNES state. - Handles BOTH rotated result dicts: the linear one-shot solve (``ksp_its`` an - int, one outer solve, no outer verdict) and the manual nonlinear Newton loop - (``nonlinear_iterations``, ``ksp_its`` a per-increment LIST, and an outer - ``converged`` flag from the residual/step-norm tests). ``ksp_reason`` is a KSP - code on both paths (the nonlinear dict carries the LAST increment's), so it is - rendered with the KSP reason table — the SNES table shares the integers but - names different outcomes. A malformed dict raises: this reader and the dicts in - ``rotated_bc.py`` are a contract, and a silent fallback here previously masked - a broken one.""" + The rotated result dict is the manual Newton loop's (``nonlinear_iterations``, + ``ksp_its`` a per-increment LIST, an outer ``converged`` flag from the + residual/step-norm tests, and ``ksp_reason`` the LAST increment's KSP code) — + rendered with the KSP reason table, since the SNES table shares the integers + but names different outcomes. A malformed dict raises: this reader and the + dict in ``rotated_bc.py`` are a contract, and a silent fallback here + previously masked a broken one.""" from underworld3.systems.solve_report import SolveReport, ksp_reason_string info = info or {} @@ -5658,9 +5656,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._block_constraint_bcs = [] # Rotated strong free-slip BCs: [(boundary, normal), ...]. Registered via # add_rotated_freeslip_bc; when non-empty, solve() delegates to - # underworld3.utilities.rotated_bc (per-node DOF rotation + strong v_n=0 + - # reaction = sigma_nn). Empty by default → the solve path is unchanged. + # underworld3.utilities.rotated_bc (per-node DOF rotation + strong v_n = u_n + # + reaction = sigma_nn). Empty by default → the solve path is unchanged. + # _rotated_freeslip_datum maps boundary → prescribed wall-normal datum + # (the non-zero `conds` of add_rotated_freeslip_bc; absent ⇒ u.n = 0). self._rotated_freeslip_bcs = [] + self._rotated_freeslip_datum = {} self._rotated_freeslip_info = None # Give the Lagrange-multiplier (lambda) block its own viscosity-scaled # Schur preconditioner. The constraint Schur complement S_lambda = C A^-1 C^T @@ -5833,14 +5834,19 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Parameters ---------- - conds : float or None, optional - Prescribed wall-normal velocity datum, in the canonical value-first - BC order (Style Charter, API conventions). Only the homogeneous - free-slip constraint :math:`\mathbf{u}\cdot\hat{\mathbf n}=0` is - implemented, so this must be zero (or ``None``, meaning zero); a - non-zero datum raises ``NotImplementedError`` (use - :meth:`add_nitsche_bc` or ``add_constraint_bc`` for prescribed - normal in/outflow). + conds : float, scalar sympy expression, or None, optional + Prescribed wall-normal velocity datum :math:`\tilde u_n` (the SCALAR + component along the outward normal), in the canonical value-first BC + order (Style Charter, API conventions). Zero (or ``None``) is pure + free-slip :math:`\mathbf{u}\cdot\hat{\mathbf n}=0`. A non-zero value + — a number, a sympy expression of ``mesh.X``, or a scalar field read + such as ``h_dot.sym[0]`` — is imposed strongly (machine precision) as + :math:`\mathbf{u}\cdot\hat{\mathbf n}=\tilde u_n`, evaluated at the + boundary nodes at each ``solve()``. On an enclosed boundary the + datum must be discretely flux-free (:math:`\oint \tilde u_n = 0`) + for incompressibility. A corner/edge node shared between rotated + boundaries has no single normal and stays at the free-slip pinning + (datum ignored there). Vector/matrix values are rejected. boundary : str Boundary label to constrain. normal : None or sympy 1×dim Matrix or array, optional @@ -5884,14 +5890,17 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): raise TypeError( f"add_rotated_freeslip_bc() requires a boundary label string; " f"got {type(boundary).__name__}") - # Value comparison, not sympy's structural ==: Float(0.0) != Integer(0) - # structurally, but both are zero data. is_zero is True only when sympy - # can PROVE zero, so unprovable symbolic data is rejected too (#336). - if conds is not None and sympy.sympify(conds).is_zero is not True: - raise NotImplementedError( - "add_rotated_freeslip_bc imposes u.n = 0 only; a non-zero " - "wall-normal datum is not implemented (use add_nitsche_bc or " - "add_constraint_bc for prescribed normal in/outflow)") + if conds is not None: + if isinstance(conds, (sympy.MatrixBase, list, tuple, np.ndarray)): + raise TypeError( + "add_rotated_freeslip_bc(conds=...) takes the SCALAR " + "wall-normal datum u.n; got a vector/matrix value") + # Value comparison, not sympy's structural ==: Float(0.0) != Integer(0) + # structurally, but both are zero data (the free-slip member of the + # family). is_zero is True only when sympy can PROVE zero, so a + # symbolic datum (field read, expression) is correctly kept. + if sympy.sympify(conds).is_zero is not True: + self._rotated_freeslip_datum[boundary] = conds self._rotated_freeslip_bcs.append((boundary, normal)) self.is_setup = False return @@ -5911,11 +5920,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): viscosity ⇒ ``J`` independent of ``v`` ⇒ the two assemblies are bit-identical ⇒ linear. - Used to fail-fast on the rotated-free-slip path, which is a single linear - solve (assemble ``J(0)``, ``F(0)`` once) and would otherwise SILENTLY - return one Newton linearisation from ``u=0`` for a nonlinear model. The - caller must have run the pre-solve preamble (auxiliary vector + constants) - so the assembly sees the correct coefficients. + Used as a LAZY guard in the rotated-free-slip loop's picard + pure-Newton + corner (a Picard warmup is meaningless for a linear residual and impossible + for pure Newton) — the loop itself needs no up-front probe, it + self-terminates. The caller must have run the pre-solve preamble + (auxiliary vector + constants) so the assembly sees the correct + coefficients. """ snes = self.snes dm = self.dm @@ -8683,17 +8693,6 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.dm.setAuxiliaryVec(self.mesh.lvec, None) self._update_constants() - # Two paths, keyed on whether the model is nonlinear (detected by probing - # the assembled Jacobian for solution-dependence — a symbolic test cannot - # see the JIT-substituted strain-rate term): - # * LINEAR model → the validated one-shot solve (assemble J(0),F(0); - # rotate; single self-contained fieldsplit KSP). Fast, and the linear - # solution is exact, so warm-start / Picard warmup add nothing and are - # correctly ignored here. - # * NONLINEAR model → the manual outer Newton/Picard loop that rotates - # F(u), J(u) and the v_n=0 constraint every iteration. It honours - # zero_init_guess (warm start), the Picard warmup count, and the - # consistent_jacobian tangent (Picard / Newton / continuation). # guard() refuses rotated free-slip, but the BC can be added AFTER arming. # Re-check here: this path never reaches the instrumentation, so an armed # guard would be silently inert -- the exact state guard() exists to refuse. @@ -8704,15 +8703,18 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): "free-slip path, which runs its own Krylov loop outside self.snes " "and cannot be bounded by it. Call unguard() first." ) - from underworld3.utilities.rotated_bc import ( - solve_rotated_freeslip, solve_rotated_freeslip_nonlinear) - if not self._residual_is_nonlinear(): - self._rotated_freeslip_info = solve_rotated_freeslip( - self, self._rotated_freeslip_bcs, verbose=verbose) - else: - self._rotated_freeslip_info = solve_rotated_freeslip_nonlinear( - self, self._rotated_freeslip_bcs, verbose=verbose, - zero_init_guess=zero_init_guess, picard=picard) + # ONE path for linear and nonlinear models: the manual outer + # Newton/Picard loop that rotates F(u), J(u) and the strong + # v_n = u_n constraint every iteration. It honours zero_init_guess + # (warm start), the Picard warmup count, and the consistent_jacobian + # tangent (Picard / Newton / continuation); a linear model converges + # after its first increment, so no up-front nonlinearity probe is + # paid (the loop self-terminates; _residual_is_nonlinear survives as + # a lazy guard inside the picard+pure-Newton corner). + from underworld3.utilities.rotated_bc import solve_rotated_freeslip + self._rotated_freeslip_info = solve_rotated_freeslip( + self, self._rotated_freeslip_bcs, verbose=verbose, + zero_init_guess=zero_init_guess, picard=picard) # This path solves via ksp.solve on the rotated operator (not self.snes), # so give it a report from the rotated result rather than leaving a stale one. _rotated_report = self._capture_rotated_report(self._rotated_freeslip_info) diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index ebe32542..1fbbbb4d 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -493,10 +493,11 @@ def _build_consistent(self): A rotated per-node constraint — the same primitive as the held lid, differing only in the constraint right-hand side, so free-slip is the ``\tilde u_n = 0`` member of the same family. It enforces the datum to machine precision (no - penalty leak) and routes the solve through the rotated LINEAR path (a direct - KSP solve) rather than a Newton line search, so an isoviscous flow does not - thrash ``newtonls``. It requires the datum to be discretely flux-free, which - :meth:`_surface_mean` now guarantees by weighting the demean by arc length. + penalty leak) through the unified rotated Newton loop: an isoviscous flow + converges in one increment (the cold-start affine lift IS the linear solve), + and a nonlinear rheology iterates with the datum held exactly at every + accepted iterate. It requires the datum to be discretely flux-free, which + :meth:`_surface_mean` guarantees by weighting the demean by arc length. ``"penalty"`` A weak natural BC. It tolerates a datum that carries a small net flux, at the @@ -515,9 +516,10 @@ def _build_consistent(self): if self.consistent_constraint == "strong": # The datum is read along the rotated per-node normal — the same deform # direction the rate was measured along, so no spurious tangential-slope term - # on a bumpy or rotating surface. - self.consistent.add_rotated_freeslip_bc(0.0, self.surface, normal=self.normal) - self.consistent._rotated_freeslip_datum = {self.surface: self._un_target.sym[0]} + # on a bumpy or rotating surface. Value-first: the ũ_n field read IS the + # conds datum, re-evaluated at the boundary nodes at each solve. + self.consistent.add_rotated_freeslip_bc( + self._un_target.sym[0], self.surface, normal=self.normal) else: n_hat = (self.normal if self.normal is not None else self.mesh.boundary_normal(self.surface)) @@ -958,19 +960,15 @@ def _solve_consistent(self, increment, dt): u_tilde = self._demean(increment / dt) self._un_target.array[...] = 0.0 self._un_target.array[self._un_target_rows, 0, 0] = u_tilde - try: - self.consistent.solve(zero_init_guess=True) - except NotImplementedError as exc: - # The rotated datum is implemented on the LINEAR rotated path only; the solve - # dispatches by a measured residual probe, so this fires exactly when the - # rheology is genuinely nonlinear. Name the knob rather than leaving the - # caller with the primitive's message. - raise NotImplementedError( - "FreeSurface: the strong rotated constraint cannot prescribe u.n = ũ_n " - "for a nonlinear rheology (the rotated datum is implemented on the linear " - "path only). Construct the manager with consistent_constraint='penalty' " - "to impose the material-surface rate weakly instead." - ) from exc + # Warm-start from the free solve: the consistent solution IS the free + # solution with the (small) material-boundary datum imposed, and the free + # solve has already converged this step. Starting there keeps a power-law + # tangent at physical strain rates — a cold start puts it at the + # regularisation floor, where the Newton line search stalls at O(0.1) + # relative residual (measured, power-law annulus acceptance run). + self.consistent.u.array[...] = self.free.u.array + self.consistent.p.array[...] = self.free.p.array + self.consistent.solve(zero_init_guess=False) def _conserve_composition(self): r"""Hold :math:`\int` (conserve integrand) fixed by a uniform shift of the diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index c0eb262f..448a7297 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -1,10 +1,14 @@ """Rotated strong free-slip for the Stokes saddle (the implementation behind ``solver.add_rotated_freeslip_bc``): build a per-node rotation Q from boundary -normals, rotate the assembled saddle Â=Q A Qᵀ / b̂=Q b, impose v_n=0 on the -rotated normal rows, solve, rotate back u=Qᵀû, remove the rigid-rotation gauge, -and expose σ_nn as the constraint reaction. - -The rotated saddle is solved by a self-contained fieldsplit-Schur KSP by default: the +normals, then drive ONE manual Newton/Picard loop that rotates the residual and +tangent every iteration (Ĵ=Q J Qᵀ, F̂=Q F), imposes v_n = ũ_n strongly on the +rotated normal rows (ũ_n = 0 is pure free-slip; a datum enters cold starts via the +first increment's affine lift, feasible iterates thereafter), rotates back u=Qᵀû, +removes the rigid-rotation gauge, and exposes σ_nn as the constraint reaction. +A linear model is simply the loop converging after its first increment — there is +no separate linear path and no up-front nonlinearity probe. + +Each increment is solved by a self-contained fieldsplit-Schur KSP by default: the velocity block is geometric FMG on the custom prolongation (``set_custom_fmg``) when a hierarchy is registered (the PREFERRED route), else GAMG tuned to the native path's settings; the Schur complement is preconditioned by the native 1/mu pressure mass @@ -318,169 +322,29 @@ def _set_rows_local(vec, row_val_map): # --------------------------------------------------------------------------- # # The rotated solve # --------------------------------------------------------------------------- # -def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, verbose=False): - """Assemble + solve the rotated strong-free-slip Stokes saddle. Fills the - solver's velocity/pressure fields with the (rotated-back, gauge-removed) - solution. - - Called from ``SNES_Stokes_SaddlePt.solve`` after ``_build`` (so the SNES/DM - exist); when used standalone it builds the solver itself. - - Returns - ------- - dict - The solve result consumed by ``boundary_normal_traction`` / - ``dynamic_topography_field`` (as their ``solve_result`` argument), with keys: - - * ``"Q"``, ``"Qt"`` — the rotation and its transpose (PETSc Mats); - * ``"A"``, ``"b"`` — the unrotated operator and RHS, kept so the reaction - ``A·u − b`` can be reconstructed (linear path only); - * ``"U"``, ``"Uhat"`` — the Cartesian and rotated-frame solutions; - * ``"normal_rows"`` — global rows of the constrained normal components; - * ``"boundaries"`` — the boundary specs the rotation was built from; - * ``"rotation_gauge_removed"`` — whether a rigid-rotation gauge was projected out; - * ``"ksp_reason"``, ``"ksp_its"`` — outer KSP converged-reason and iteration count; - * ``"rnorm"`` — true rotated residual ‖·û − b̂‖ (feeds the solve report). - """ - if getattr(solver, "snes", None) is None: - solver._setup_pointwise_functions() - solver._setup_discretisation() - solver._setup_solver() - dm = solver.dm - snes = solver.snes - - # Assemble the operator FIRST so its parallel row layout is final before we build - # Q against it (A = exact Jacobian at 0 — linear; b = -F(0)). The Pmat is - # assembled alongside: its p-p block is the native 1/mu pressure mass (the - # DS JacobianPreconditioner term) that preconditions the Schur complement — - # petsc4py's computeJacobian(x, J) would silently pass J as its own Pmat and - # the mass block would never be assembled. - snes.setUp() - U0 = dm.getGlobalVec() - U0.set(0.0) - J, Jp = snes.getJacobian()[:2] - snes.computeJacobian(U0, J, Jp) - Aorig = J.copy() - F0 = dm.getGlobalVec() - snes.computeFunction(U0, F0) - b = F0.copy() - b.scale(-1.0) - # borrowed temporaries → return to pool - dm.restoreGlobalVec(U0) - dm.restoreGlobalVec(F0) - - datum_specs = getattr(solver, "_rotated_freeslip_datum", None) - Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) - - # rotate:  = Q A Qᵀ, b̂ = Q b - Ahat = Aorig.ptap(Qt) - bhat = b.duplicate() - Q.mult(b, bhat) - - # constrain rotated normal rows (v_n = datum, datum=0 for pure free-slip): zero the - # matrix rows/cols and set the RHS at those rows to the datum. The constraint - # diagonal is the mean |diag(A_vv)| rather than 1.0: unit diagonals amid O(eta/h^2) - # viscous entries put a spectrum outlier on every rotated boundary node, poisoning - # diagonal-based Schur approximations and MG smoothing in the boundary strip. Any - # positive diagonal is exact — the solution rows are set explicitly. - diag = _velocity_diag_scale(Ahat, solver) - if datum_map: - # v_n = ũ_n: pass the datum (rotated-frame normal component = v_n) as x so - # zeroRowsColumns(rows, diag, x, b) sets b[rows]=diag·x[rows] AND subtracts the - # eliminated columns' contribution A[:,cols]·x from the other rows. x=0 - # reproduces the pure free-slip RHS exactly, so datum=0 is bit-identical. - xhat = bhat.duplicate() - xhat.zeroEntries() - rs, re = xhat.getOwnershipRange() - xa = xhat.getArray() - for g, val in datum_map.items(): - if rs <= g < re: - xa[g - rs] = val - xhat.setArray(xa) - Ahat.zeroRowsColumns(normal_rows, diag=diag, x=xhat, b=bhat) - xhat.destroy() - else: - # zeroRowsColumns takes GLOBAL row indices; the RHS write goes through - # _zero_rows_local (ownership-relative indexing — the np>1 crash class). - Ahat.zeroRowsColumns(normal_rows, diag=diag) - _zero_rows_local(bhat, normal_rows) - - # ITERATIVE by default (LU is almost never right): a self-contained fieldsplit- - # Schur solve whose velocity block is geometric FMG on the custom prolongation - # when a hierarchy is registered (set_custom_fmg), else GAMG. Direct LU only when - # explicitly opted in via solver._rotated_use_lu. - if getattr(solver, "_rotated_use_lu", False): - # Pin one pressure DOF (datum) — row only, keeps the B^T coupling. Only the - # direct solve needs it; the iterative path fixes the pressure gauge with the - # constant-pressure null space instead. - # TODO(BUG): the datum search is a naive per-rank scan of the rank's OWN - # global-section chart, so at np>1 each rank pins a different pressure DOF - # (or none, pin=None → zeroRows([None]) fails). Parallel-unsafe; tolerated - # only because LU is opt-in via solver._rotated_use_lu. - gsec = dm.getGlobalSection() - pin = None - pS, pE = gsec.getChart() - for q in range(pS, pE): - if gsec.getFieldDof(q, _PRESSURE_FIELD) > 0 \ - and gsec.getFieldOffset(q, _PRESSURE_FIELD) >= 0: - pin = gsec.getFieldOffset(q, _PRESSURE_FIELD) - break - Ahat.zeroRows([pin], diag=1.0) - if pin is not None: - _zero_rows_local(bhat, [pin]) - ksp = PETSc.KSP().create() - ksp.setOperators(Ahat) - ksp.setType("preonly") - pc = ksp.getPC() - pc.setType("lu") - pc.setFactorSolverType("mumps") - Uhat = dm.createGlobalVec() # returned in the result dict → own it - ksp.solve(bhat, Uhat) - ksp_reason = ksp.getConvergedReason() - ksp_its = ksp.getIterationNumber() - _warn_if_ksp_diverged(ksp, kind="rotated direct-LU") - else: - Mp = _pressure_mass_schur_pmat(solver) - Uhat, ksp_reason, ctx = _solve_rotated_iterative( - solver, Ahat, bhat, Q, Qt, normal_rows, verbose=verbose, Mp=Mp) - ksp_its = ctx["ksp"].getIterationNumber() - _destroy_rotated_ksp_ctx(ctx) - - # a prescribed normal datum: the iterative path zeros the rotated normal rows of - # the solution (its exact v_n=0 cleanup); restore v_n = ũ_n before rotating back - # (the constrained matrix already drove the interior — only these rows were reset). - # Restored BEFORE the residual report below, so the reported norm describes the - # solution actually returned (the datum rows are consistent with b̂ by construction). - if datum_map: - _set_rows_local(Uhat, datum_map) - - # True rotated residual ‖·û − b̂‖ for the solve report (one matvec). Computed - # explicitly rather than read off the KSP: preonly/LU never computes a norm, and - # the iterative norm can be the preconditioned one. - _res = bhat.duplicate() - Ahat.mult(Uhat, _res) - _res.axpy(-1.0, bhat) - rnorm = float(_res.norm()) - _res.destroy() - - # rotate back u = Qᵀ û (U is returned in the result dict → create, don't - # borrow from the pool) - U = dm.createGlobalVec() - Qt.mult(Uhat, U) - - removed = _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) - - return {"Q": Q, "Qt": Qt, "A": Aorig, "b": b, "U": U, "Uhat": Uhat, - "normal_rows": normal_rows, "boundaries": list(boundaries), - "rotation_gauge_removed": removed, "ksp_reason": ksp_reason, - "ksp_its": ksp_its, "rnorm": rnorm} +def _naive_pressure_pin(dm): + """One owned pressure DOF (datum) for the direct-LU gauge pin — row only, so + the B^T coupling is kept. Only the direct solve needs it; the iterative path + fixes the pressure gauge with the constant-pressure null space instead. + + TODO(BUG): a naive per-rank scan of the rank's OWN global-section chart, so at + np>1 each rank pins a different pressure DOF (or none → the caller must skip). + Parallel-unsafe; tolerated only because LU is opt-in via + ``solver._rotated_use_lu`` (a serial PC-free diagnostic).""" + gsec = dm.getGlobalSection() + pS, pE = gsec.getChart() + for q in range(pS, pE): + if gsec.getFieldDof(q, _PRESSURE_FIELD) > 0 \ + and gsec.getFieldOffset(q, _PRESSURE_FIELD) >= 0: + return gsec.getFieldOffset(q, _PRESSURE_FIELD) + return None def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge): """Remove the rigid-rotation gauge (if it is a genuine null space of the constrained problem), scatter the composite global vector ``U`` into the velocity/pressure fields, and refresh the enhanced-variable caches. Shared by - the linear one-shot and the nonlinear driver. Returns whether the gauge was + every exit of the rotated solve. Returns whether the gauge was removed.""" dm = solver.dm # Remove the rigid-rotation gauge — every mode that is a genuine null space of @@ -536,7 +400,7 @@ def _finalize_rotated_solution(solver, U, Q, normal_rows, remove_rotation_gauge) def _gather_fields_to_global(solver): """Composite global vector built from the solver's current velocity/pressure - field values (the warm-start initial guess for the nonlinear driver).""" + field values (the warm-start initial guess for the rotated Newton loop).""" dm = solver.dm U = dm.createGlobalVec() U.set(0.0) @@ -547,25 +411,32 @@ def _gather_fields_to_global(solver): return U -def _project_out_normal_component(u, Q, Qt, normal_rows): - """Impose the strong ``v_n = 0`` constraint exactly on the composite vector - ``u``, in place: rotate to the boundary frame (``û = Q u``), zero the - constrained normal rows, rotate back (``u = Qᵀ û``). Q is orthogonal, so this - is the exact projection onto the constraint-satisfying subspace.""" +def _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map=None): + """Impose the strong affine constraint ``v_n = ũ_n`` exactly on the composite + vector ``u``, in place: rotate to the boundary frame (``û = Q u``), zero the + constrained normal rows, set the prescribed-datum rows to ``sgn·ũ_n`` + (``datum_map``; empty/None ⇒ pure free-slip ``v_n = 0``), rotate back + (``u = Qᵀ û``). Q is orthogonal, so this is the exact affine projection onto + the constraint-satisfying set. Keeping every Newton iterate feasible this way + is what lets the increment problem stay HOMOGENEOUS (``n̂·δ = 0``) — the + datum never touches the tangent or the increment RHS.""" uh = u.duplicate() # transient projection buffer Q.mult(u, uh) _zero_rows_local(uh, normal_rows) + if datum_map: + _set_rows_local(uh, datum_map) Qt.mult(uh, u) uh.destroy() def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, - max_halvings=8): + datum_map=None, max_halvings=8): """Backtracking line search on ‖F̂‖ for the rotated Newton/Picard update ``u + α d`` (full step α=1 first, halved on failure). Cheap insurance far from the solution; α=1 is accepted immediately near it. Every trial iterate is - projected back onto the strong v_n=0 constraint before its residual is - measured. + snapped back onto the strong ``v_n = ũ_n`` constraint (``datum_map``; free-slip + when empty) before its residual is measured, so ‖F̂‖ is always evaluated at a + feasible point. Owns all its temporaries' destroys: on acceptance the input ``u`` is destroyed and replaced by the accepted iterate. Returns ``(u, improved)``; ``improved`` @@ -576,7 +447,7 @@ def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, for _ls in range(max_halvings): utry = u.copy() utry.axpy(alpha, d) - _project_out_normal_component(utry, Q, Qt, normal_rows) + _impose_normal_constraint(utry, Q, Qt, normal_rows, datum_map) Ftry = rotated_residual(utry) fnorm = Ftry.norm() Ftry.destroy() @@ -588,15 +459,25 @@ def _backtracking_line_search(u, d, rnorm, rotated_residual, Q, Qt, normal_rows, return u, False -def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=True, - verbose=False, zero_init_guess=True, picard=0, - rtol=None, atol=1.0e-11, stol=1.0e-8, max_it=50): - """Nonlinear rotated strong-free-slip solve: a manual outer Newton/Picard loop - that rotates the residual F(u), the Jacobian J(u) and the v_n=0 constraint EVERY +def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, + verbose=False, zero_init_guess=True, picard=0, + rtol=None, atol=1.0e-11, stol=1.0e-8, max_it=50): + """THE rotated strong-free-slip solve (linear and nonlinear models alike): a + manual outer Newton/Picard loop that rotates the residual F(u), the Jacobian + J(u) and the strong ``v_n = ũ_n`` constraint (``ũ_n = 0`` for pure free-slip; + a prescribed wall-normal datum via ``solver._rotated_freeslip_datum``) EVERY iteration, reusing the validated self-contained rotated fieldsplit-Schur solve (``_solve_rotated_iterative``, incl. custom geometric FMG / GAMG velocity block and the rotated coupled null space) for each Newton increment. + A LINEAR model needs no separate path: the first increment from a cold start is + assembled at rest and already solves the problem (with a datum it is exactly the + affine-lift solve), so the loop converges at the next residual check — one + linear solve, plus two cheap residual evaluations. There is therefore no + up-front nonlinearity probe; the loop self-terminates. Direct LU per increment + (a serial PC-free diagnostic arbiter, e.g. for TI anisotropy experiments) is + opt-in via ``solver._rotated_use_lu``. + Why a manual loop rather than ``snes.solve()``: the rotated operator ``Q A Qᵀ`` (a ``ptap`` result) carries no DM field information, so PETSc's DM-coupled fieldsplit + geometric-MG cannot precondition it (SUBPC_ERROR); the increment @@ -605,14 +486,20 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T strong constraint exactly at every iterate. Each iteration (unknown carried in the CARTESIAN frame ``u``; the increment is - solved in the rotated frame ``û = Q u``): + solved in the rotated frame ``û = Q u``). Every ACCEPTED iterate is FEASIBLE + (``v_n = ũ_n`` imposed exactly), so increments satisfy the homogeneous + constraint ``n̂·δ = 0`` and the datum never enters the tangent. The one + exception is a cold start with a non-zero datum, where the first increment + carries the datum jump through the affine lift (see the initial-guess comment + in the body — snapping the zero state onto the datum creates a boundary-strip + strain state whose nonlinear tangent is unusable): * ``F = computeFunction(u)`` → Cartesian residual (native essential BCs on other boundaries already applied by the DM); * ``F̂ = Q F``, zero ``F̂`` at the constrained normal rows (the constraint - residual is 0 there); converge on ‖F̂‖; + residual is 0 there — the iterate is feasible); converge on ‖F̂‖; * ``Ĵ = Q J(u) Qᵀ`` with ``zeroRowsColumns(normal_rows)``; * solve ``Ĵ δ̂ = −F̂``, ``δ = Qᵀ δ̂``, with a ‖F̂‖ backtracking line search; - * ``u += α δ`` and re-impose ``v_n = 0`` exactly. + * ``u += α δ`` and re-impose ``v_n = ũ_n`` exactly. The tangent used by ``computeJacobian`` is the solver's own (``consistent_jacobian`` → Picard / Newton / continuation), so the rotated loop inherits the same tangent @@ -643,7 +530,8 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T * ``"Q"``, ``"Qt"`` — the rotation and its transpose (PETSc Mats); * ``"reaction"`` — the converged Cartesian residual F(u), stashed as the - constraint reaction for σ_nn recovery (this path has no ``"A"``/``"b"``); + constraint reaction for σ_nn recovery (for a linear residual this equals + the assembled ``A·u − b`` exactly); * ``"U"`` — the Cartesian solution; * ``"normal_rows"`` — global rows of the constrained normal components; * ``"boundaries"`` — the boundary specs the rotation was built from; @@ -672,37 +560,75 @@ def solve_rotated_freeslip_nonlinear(solver, boundaries, remove_rotation_gauge=T mode = getattr(solver, "consistent_jacobian", False) continuation = (mode == "continuation") if picard and not continuation and mode is True: - raise NotImplementedError( - f"rotated free-slip: a Picard->Newton warmup (picard={picard}) with the " - "consistent Newton tangent (consistent_jacobian=True) requires " - "consistent_jacobian='continuation' — with pure Newton the frozen (Picard) " - "tangent needed for the warmup is not compiled in. Use " - "consistent_jacobian='continuation' (staged Picard then Newton), or drop " - "picard to run pure Newton.") + # The ONLY place linearity must be known in advance (hence the lazy probe — + # two trial assemblies, paid only by this contradictory configuration): a + # Picard warmup is meaningless for a linear residual, so ignore it there; + # for a genuinely nonlinear model pure Newton has no frozen tangent + # compiled in to warm up with, so fail loudly. + if solver._residual_is_nonlinear(): + raise NotImplementedError( + f"rotated free-slip: a Picard->Newton warmup (picard={picard}) with the " + "consistent Newton tangent (consistent_jacobian=True) requires " + "consistent_jacobian='continuation' — with pure Newton the frozen (Picard) " + "tangent needed for the warmup is not compiled in. Use " + "consistent_jacobian='continuation' (staged Picard then Newton), or drop " + "picard to run pure Newton.") + picard = 0 switch_rtol = max(float(getattr(solver, "newton_switch_rtol", 1.0e-2)), rtol) if continuation: solver._set_newton_alpha(0.0) # start in the Picard phase - # Q, the custom-FMG prolongation and the coupled null space depend only on the - # geometry / normals (NOT the solution), so build them ONCE and reuse each step. - if getattr(solver, "_rotated_freeslip_datum", None): - raise NotImplementedError( - "a prescribed non-zero wall-normal datum (u.n = ũ_n) on rotated free-slip is " - "implemented only for the LINEAR solve path (solve_rotated_freeslip); the " - "nonlinear path still imposes u.n = 0. Use a linear model for a prescribed " - "normal velocity, or extend this path.") - Q, Qt, normal_rows, _ = build_rotation(solver, boundaries) - custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) - nsp = _rotated_nullspace(solver, Q, normal_rows) - - # initial guess (cartesian, composite): warm-start from the fields or zero, then - # impose v_n=0 exactly on it so the iteration starts feasible. + # Q, the custom-FMG prolongation, the coupled null space and the datum values + # depend only on the geometry / normals / prescribed surface field (NOT the + # solution), so build them ONCE and reuse each step. The datum ũ_n enters the + # iteration ONLY through the feasible-iterate projection below: every iterate + # carries v_n = ũ_n exactly, so the Newton increment satisfies the HOMOGENEOUS + # constraint n̂·δ = 0 and the per-iteration operator treatment (zeroRowsColumns + # with no lift, tangent transparency, custom_Pl, nullspace) is unchanged from + # pure free-slip. + datum_specs = getattr(solver, "_rotated_freeslip_datum", None) + Q, Qt, normal_rows, datum_map = build_rotation(solver, boundaries, datum_specs) + # Direct LU per increment: no FMG prolongation / null space to build — the + # gauge is fixed by the naive pressure pin instead (see _naive_pressure_pin). + use_lu = bool(getattr(solver, "_rotated_use_lu", False)) + custom_Pl = None if use_lu else _build_rotated_custom_Pl(solver, Q, normal_rows) + # The null space is built INSIDE the loop, after the first Jacobian assembly: + # _mode_satisfies_constraints verifies each candidate mode against the + # ASSEMBLED operator (‖J·m‖ ≈ 0), which an unassembled J cannot support. + nsp = None + + # initial guess (cartesian, composite): warm-start from the fields or zero. + # A prescribed datum on a COLD start is imposed through the FIRST increment's + # affine lift (zeroRowsColumns x/b at the rest-state tangent), + # NOT by snapping the zero state: v_n = ũ_n over a zero interior manufactures an + # extreme boundary-strip strain rate whose shear-thinning tangent linearises so + # badly that no line-searched step descends (measured: first step 4 orders too + # large, f(α) > f(0) down to α=2⁻¹¹). Assembling the first tangent at rest and + # letting the increment carry the datum jump gives the smooth regularised- + # viscosity response instead; every later iterate is feasible and the + # increments are homogeneous. A warm start is assumed smooth (the previous + # converged state) and takes the exact affine snap directly. + # x̂ (the prescribed rotated-frame datum as a composite vector) does double + # duty, as in the linear path: its collective PETSc norm is the rank-consistent + # datum-activity flag (datum_map is rank-local — per-rank branching on it + # desyncs the zeroRowsColumns variant and the line-search collectives, the + # np>1 deadlock class), and it is the ready-made lift vector for the + # cold-start first increment. + xhat = dm.createGlobalVec() + xhat.set(0.0) + _set_rows_local(xhat, datum_map) + datum_active = xhat.norm() > 0.0 + lift_datum = datum_active and zero_init_guess + if not lift_datum: + xhat.destroy() + xhat = None if zero_init_guess: u = dm.createGlobalVec() u.set(0.0) else: u = _gather_fields_to_global(solver) - _project_out_normal_component(u, Q, Qt, normal_rows) + if not lift_datum: + _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map) J, Jp = snes.getJacobian()[:2] pres_is = solver._subdict["pressure"][0] @@ -728,7 +654,24 @@ def rotated_residual(uvec, keep_cartesian=False): _zero_rows_local(Fh, normal_rows) return Fh + # Convergence reference: max(initial residual, REST-STATE residual ‖F̂(0)‖). + # A warm start's own initial residual is small, and rtol relative to it + # demands ever-more absolute accuracy the better the guess (the warm-start + # rtol trap: the FS time loop stalled at rel~2e-3 of its warm start while + # being far below tolerance on the physical scale). The rest-state residual + # is the problem's intrinsic forcing scale; a cold start's r0 IS that scale, + # so nothing changes there. + fnorm_rest = None + if not zero_init_guess: + z = dm.getGlobalVec() + z.set(0.0) + Fz = rotated_residual(z) + fnorm_rest = Fz.norm() + Fz.destroy() + dm.restoreGlobalVec(z) + r0 = None + ref = None last_reason = 0 iters = 0 converged = False @@ -738,12 +681,13 @@ def rotated_residual(uvec, keep_cartesian=False): rnorm = Fhat.norm() if r0 is None: r0 = rnorm + ref = max(r0, fnorm_rest) if fnorm_rest is not None else r0 if verbose: mpi.pprint(f"[rotated_bc] nonlinear iter {iters:2d} |F̂|={rnorm:.6e} " - f"rel={rnorm/(r0+1e-300):.3e} [{phase}]") - # residual convergence (relative to the initial residual, plus an absolute + f"rel={rnorm/(ref+1e-300):.3e} [{phase}]") + # residual convergence (relative to the reference scale, plus an absolute # floor so an already-converged warm start does not chase machine noise). - if rnorm <= rtol * r0 + atol: + if rnorm <= rtol * ref + atol: converged = True Fhat.destroy() break @@ -761,18 +705,52 @@ def rotated_residual(uvec, keep_cartesian=False): Ahat = J.ptap(Qt) else: J.ptap(Qt, result=Ahat) # same nonzero pattern → in-place refresh - if ctx is None: - Mp = _pressure_mass_schur_pmat(solver) - elif Mp is not None: - Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent + if not use_lu: + if ctx is None: + Mp = _pressure_mass_schur_pmat(solver) + nsp = _rotated_nullspace(solver, Q, normal_rows) # J assembled above + elif Mp is not None: + Jp.createSubMatrix(pres_is, pres_is, submat=Mp) # viscosity may be u-dependent if diag_scale is None: diag_scale = _velocity_diag_scale(Ahat, solver) - Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) bhat = Fhat.copy() bhat.scale(-1.0) - dhat, last_reason, ctx = _solve_rotated_iterative( - solver, Ahat, bhat, Q, Qt, normal_rows, - custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) + if lift_datum and iters == 0: + # cold-start lift: the increment FROM ZERO carries the datum jump — + # b̂[rows] = diag·x̂ and the eliminated columns' A[:,cols]·x̂ are + # subtracted from the other rows. For a LINEAR model this increment IS + # the whole solve; the loop converges at the next residual check. + Ahat.zeroRowsColumns(normal_rows, diag=diag_scale, x=xhat, b=bhat) + xhat.destroy() + xhat = None + else: + Ahat.zeroRowsColumns(normal_rows, diag=diag_scale) + if use_lu: + if ctx is None: + ksp_lu = PETSc.KSP().create(comm=dm.comm) + ksp_lu.setType("preonly") + pc_lu = ksp_lu.getPC() + pc_lu.setType("lu") + pc_lu.setFactorSolverType("mumps") + ctx = {"ksp": ksp_lu, "Mp": None, + "pin": _naive_pressure_pin(dm)} + pin = ctx["pin"] + if pin is not None: + Ahat.zeroRows([pin], diag=1.0) + _zero_rows_local(bhat, [pin]) + ctx["ksp"].setOperators(Ahat) # values changed in place → refactor + dhat = Ahat.createVecRight() + dhat.set(0.0) + ctx["ksp"].solve(bhat, dhat) + last_reason = ctx["ksp"].getConvergedReason() + _warn_if_ksp_diverged(ctx["ksp"], kind="rotated direct-LU") + # parity with the iterative path's exact v_n cleanup (LU is exact only + # to factorisation round-off at the decoupled constraint rows) + _zero_rows_local(dhat, normal_rows) + else: + dhat, last_reason, ctx = _solve_rotated_iterative( + solver, Ahat, bhat, Q, Qt, normal_rows, + custom_Pl=custom_Pl, nsp=nsp, Mp=Mp, verbose=False, ctx=ctx) lin_its.append(ctx["ksp"].getIterationNumber()) d = dm.createGlobalVec() Qt.mult(dhat, d) @@ -780,11 +758,34 @@ def rotated_residual(uvec, keep_cartesian=False): # are at the solution — the exit for a warm start that is already converged # (otherwise the relative test above, with a tiny r0, chatters near machine # level). ‖u‖=0 on a cold start ⇒ this never fires prematurely (d is large). + # A tiny step alone does NOT prove convergence — a stiff tangent (power-law + # at the regularisation floor) also produces tiny increments at a large + # residual — so the exit is VERIFIED against the reference scale (which + # includes the rest-state residual ‖F̂(0)‖ for warm starts): a warm start + # at the solution passes; a stagnating crawl does not. On verification + # failure the tiny step goes through the NORMAL line search instead: if it + # still improves, the crawl proceeds (bounded by max_it); if not, the + # stall exit ends the loop honestly. step_converged = d.norm() <= stol * (u.norm() + 1e-30) + if step_converged: + step_converged = rnorm <= rtol * ref + atol + if not step_converged: + mpi.pprint(f"[rotated_bc] step-norm at rel |F̂| = " + f"{rnorm / (ref + 1e-300):.2e} of the reference " + f"residual — stagnation-or-crawl, continuing to iterate.") improved = False - if not step_converged: + if lift_datum and iters == 0: + # accept the lift step unconditionally (it is the smooth rest-state + # response — the standard first-Picard warmup; a line search against + # ‖F̂(0)‖, which ignores the datum mismatch, would be meaningless) and + # snap the datum rows exactly (the KSP cleanup zeroed them). + u.axpy(1.0, d) + _impose_normal_constraint(u, Q, Qt, normal_rows, datum_map) + step_converged = False + improved = True + elif not step_converged: u, improved = _backtracking_line_search( - u, d, rnorm, rotated_residual, Q, Qt, normal_rows) + u, d, rnorm, rotated_residual, Q, Qt, normal_rows, datum_map) # every per-iteration temporary dies HERE, on all exit paths dhat.destroy() d.destroy() @@ -811,15 +812,17 @@ def rotated_residual(uvec, keep_cartesian=False): # meeting the residual / step-norm criteria. Warn — as the standard SNES path does # on divergence — so an unconverged iterate left in the fields is not silent. if not converged: - rel = (rnorm / (r0 + 1e-300)) if r0 is not None else float("nan") + rel = (rnorm / (ref + 1e-300)) if ref is not None else float("nan") mpi.pprint(f"[rotated_bc] WARNING: nonlinear rotated free-slip did NOT converge " - f"in {newton_its} iterations (rel |F̂| = {rel:.2e}); the fields hold " - f"the last (unconverged) iterate.") + f"in {newton_its} iterations (rel |F̂| = {rel:.2e} of the reference " + f"residual); the fields hold the last (unconverged) iterate.") Fc.destroy() # residual output buffer (reaction persists in the result dict) _destroy_rotated_ksp_ctx(ctx) # KSP/PC + the owned Schur pmat if Ahat is not None: Ahat.destroy() # the reused rotated operator + if xhat is not None: + xhat.destroy() # unused lift vector (loop exited before it) removed = _finalize_rotated_solution(solver, u, Q, normal_rows, remove_rotation_gauge) return {"Q": Q, "Qt": Qt, "reaction": reaction, "U": u, @@ -833,7 +836,7 @@ def rotated_residual(uvec, keep_cartesian=False): 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 no hierarchy is registered. Depends only on Q and - the mesh (NOT the solution), so the nonlinear driver builds it ONCE and reuses + the mesh (NOT the solution), so the rotated Newton loop builds it ONCE and reuses it across Newton iterations (the prolongation build is the expensive part).""" if getattr(solver, "_custom_mg", None) is None: return None @@ -915,8 +918,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal complement for enclosed domains (the hand-built IS fieldsplit does not inherit it from the operator). - ``custom_Pl`` / ``nsp`` may be PREBUILT (nonlinear driver: build once, reuse each - Newton step); when None they are built here (linear one-shot). + ``custom_Pl`` / ``nsp`` are built by the CALLER (once, reused each Newton step; + the null-space mode verification needs the assembled Jacobian, which only the + caller can guarantee); None simply means "no hierarchy" / "no null modes". Returns ``(Uhat, reason, ctx)``. Passing ``ctx`` back in reuses the KSP/PC across Newton iterations — the fieldsplit ISs, Schur USER pmat and FMG @@ -929,12 +933,9 @@ def _solve_rotated_iterative(solver, Ahat, bhat, Q, Qt, normal_rows, verbose=Fal pres_is = solver._subdict["pressure"][0] if ctx is None: - if custom_Pl is None: - custom_Pl = _build_rotated_custom_Pl(solver, Q, normal_rows) - - # rotated coupled null space (pressure-const ⊕ Q·rotation) on the operator - if nsp is None: - nsp = _rotated_nullspace(solver, Q, normal_rows) + # rotated coupled null space (pressure-const ⊕ Q·rotation) on the operator — + # built by the CALLER (the Newton loop, after the first Jacobian assembly: + # the mode verification needs the assembled operator). if nsp is not None: Ahat.setNullSpace(nsp) Ahat.setTransposeNullSpace(nsp) @@ -1116,10 +1117,23 @@ def _rotated_nullspace(solver, Q, normal_rows): def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): - """True iff the rigid-body mode ``tg`` satisfies all rotated v_n=0 - constraints — i.e. Q·tg is ~0 on every constrained normal row. (A closed - circular boundary admits its one rotation; a full spherical shell admits all - three; straight/partial walls pin them.) + """True iff the rigid-body mode ``tg`` is a genuine null mode of the + constrained problem: it satisfies all rotated v_n=0 constraints (Q·tg ~0 on + every constrained normal row) AND it is a null vector of the assembled + operator (‖J·tg‖ ~0 against the velocity diagonal scale). The operator test + is what catches pinning the rotated rows cannot see — an ESSENTIAL no-slip + on another boundary leaves the mode tangential to the rotated wall (first + test passes) while the eliminated operator is NOT null on it; admitting it + then projects an irreducible component out of every increment RHS, and the + Newton residual floors at that component's magnitude instead of converging + (measured: rel ~2e-5 plateau on the essential-inner + rotated-outer annulus). + (A closed circular free-slip boundary admits its one rotation; a full + spherical shell admits all three; straight/partial walls and any essential + BC pin them.) + + The caller must have ASSEMBLED the solver's Jacobian at some iterate; an + unassembled J (norm 0) skips the operator test rather than passing every + mode through a vacuous 0 ≈ 0. COLLECTIVE: every rank runs the same global-vector ops. Do NOT early-return on a per-rank ``not normal_rows`` — in parallel a rank may own no boundary node @@ -1141,14 +1155,28 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8): # transient duplicates tr.destroy() trc.destroy() - return viol < tol + if viol >= tol: + return False + # operator-nullity: rigid rotation has exactly zero strain in the discrete + # space (P2 contains linear fields, affine cells integrate the form exactly), + # so a genuine null mode gives assembly round-off; a pinned mode leaves O(1) + # boundary-strip rows. + J = solver.snes.getJacobian()[0] + Jm = tg.duplicate() + J.mult(tg, Jm) + jn = Jm.norm() + Jm.destroy() + if jn == 0.0 and J.norm() == 0.0: # J never assembled → cannot verify + return True + op_viol = jn / (_velocity_diag_scale(J, solver) * (tg.norm() + 1e-30)) + return op_viol < tol def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): """Boundary normal traction σ_nn on `boundary` from the constraint reaction of the last rotated-free-slip solve (``solve_result`` is the dict returned by - ``solve_rotated_freeslip`` / ``solve_rotated_freeslip_nonlinear`` — see their - Returns sections for the keys). Returned mean-removed (the ρg·h gauge), as + ``solve_rotated_freeslip`` — see its + Returns section for the keys). Returned mean-removed (the ρg·h gauge), as ``(xs, sigma)`` with one entry per boundary velocity node on this rank. σ_nn is recovered from the CARTESIAN nodal reaction r_c = A·u − b: the nodal load @@ -1173,20 +1201,10 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): """ dm = solver.dm dim = solver.mesh.dim - # Cartesian nodal reaction r_c = F(u) at the converged state. The nonlinear - # driver stashes it directly (the final ``computeFunction`` residual); the linear - # one-shot reconstructs it as A·u−b (with A=J(0), b=−F(0), F affine ⇒ A·u−b=F(u)). - if solve_result.get("reaction") is not None: - rc = solve_result["reaction"] # owned by the result dict — do NOT destroy - own_rc = False - else: - A = solve_result["A"] - b = solve_result["b"] - U = solve_result["U"] - rc = A.createVecLeft() - A.mult(U, rc) - rc.axpy(-1.0, b) - own_rc = True + # Cartesian nodal reaction r_c = F(u) at the converged state, stashed by the + # solve driver (the final ``computeFunction`` residual; for a linear residual + # this equals the assembled A·u−b exactly). + rc = solve_result["reaction"] # owned by the result dict — do NOT destroy rcl = dm.getLocalVec() dm.globalToLocal(rc, rcl) rca = np.asarray(rcl.getArray()) @@ -1205,8 +1223,6 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): xs.append(_point_coord(dm, dim, cvec, csec, v0, v1, q)) Rn.append(float(np.dot(nrm, rcv))) # R_i = n̂·r_c (corner-correct) dm.restoreLocalVec(rcl) - if own_rc: - rc.destroy() # reconstructed A·u−b (linear path) — free it xs = np.array(xs) Rn = np.array(Rn) # σ_nn = −(nodal reaction), de-smeared by the SHARED boundary-mass primitive and diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index a11a5815..f15e9ef8 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -41,8 +41,7 @@ def test_rotated_datum_prescribed_normal_partition_independent(): blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation) - s.add_rotated_freeslip_bc(0.0, "Upper", normal=nhat) - s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) # u.n = cos(theta), mean-zero s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" s.tolerance = 1.0e-9 @@ -57,3 +56,44 @@ def test_rotated_datum_prescribed_normal_partition_independent(): vv = float(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate()) assert abs(vv - _INT_VV_REF) / _INT_VV_REF < 1.0e-6, \ f"solve energy is partition-dependent: ∫v·v={vv:.8e} vs ref {_INT_VV_REF}" + + +def test_rotated_datum_nonlinear_parallel(): + """The NONLINEAR rotated datum path in parallel: a power-law annulus with + u.n = cos(theta) through the manual Newton loop (feasible iterates, cold-start + lift, collective datum-activity flag). Guards the np>1 collective-sequence + class: the lift branch, the line-search residual evaluations and the + nullspace-mode verification are all collective, and a rank-divergent branch + deadlocks rather than failing cleanly (hence the module timeout).""" + RI, RO = 0.5, 1.0 + mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x ** 2 + y ** 2) + nhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable("Vnl64", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable("Pnl64", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + gm = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + e = 0.5 * (gm + gm.T) + eII = sympy.sqrt(0.5 * (e[0, 0] ** 2 + e[1, 1] ** 2) + e[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) + s.consistent_jacobian = True + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-7 + s.solve() + + info = s._rotated_freeslip_info + assert info["converged"], "nonlinear parallel datum solve did not converge" + assert info["nonlinear_iterations"] > 1, "did not genuinely iterate" + # datum imposed, measured with parallel-safe boundary quadrature only + vn = (v.sym[0] * x + v.sym[1] * y) / r + num = float(uw.maths.BdIntegral(mesh=mesh, fn=(vn - x / r) ** 2, boundary="Upper").evaluate()) + den = float(uw.maths.BdIntegral(mesh=mesh, fn=(x / r) ** 2, boundary="Upper").evaluate()) + relL2 = np.sqrt(max(num, 0.0) / den) + assert relL2 < 1.0e-3, f"nonlinear u.n=cos(theta) not imposed (relL2={relL2:.3e})" diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index b74ee522..02d81ada 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -98,8 +98,8 @@ def test_rotated_freeslip_spherical_shell_3d(): info = s._rotated_freeslip_info assert info["ksp_reason"] > 0, f"rotated KSP diverged: {info['ksp_reason']}" - assert info["ksp_its"] <= 25, ( - f"Schur iteration blow-out: {info['ksp_its']} outer its " + assert max(info["ksp_its"]) <= 25, ( + f"Schur iteration blow-out: {info['ksp_its']} outer its per increment " "(1/mu-mass Schur preconditioning regressed?)") assert info["rotation_gauge_removed"], "3D rotation gauge not detected/removed" @@ -611,11 +611,11 @@ def test_rotated_freeslip_dynamic_topography_field(): # --- Prescribed non-zero wall-normal datum (u.n = ũ_n) --------------------------- -# The rotated constraint imposes u.n = datum strongly: datum=0 is pure free-slip (the -# held lid); a non-zero datum is the "consistent" material-surface velocity. Both share -# the same rotated matrix and differ only in the constraint RHS, so datum=0 must be -# BIT-IDENTICAL to plain free-slip. (Set via solver._rotated_freeslip_datum until the -# add_rotated_freeslip_bc datum argument lands — underworldcode/underworld3 tracking.) +# The rotated constraint imposes u.n = datum strongly, via the value-first conds +# argument: conds=0 is pure free-slip (the held lid); a non-zero conds is the +# "consistent" material-surface velocity. Both share the same rotated matrix and +# differ only in the constraint RHS, so an explicit zero datum must reproduce plain +# free-slip exactly. def _annulus_datum_solve(mode, tag): RI, RO = 0.5, 1.0 @@ -631,11 +631,8 @@ def _annulus_datum_solve(mode, tag): blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation gauge) - s.add_rotated_freeslip_bc(0, "Upper", normal=nhat) - if mode == "zero": - s._rotated_freeslip_datum = {"Upper": 0.0} - elif mode == "cos": - s._rotated_freeslip_datum = {"Upper": x / r} # u.n = cos(theta), mean-zero => ∮u.n=0 + datum = {"plain": 0, "zero": 0.0, "cos": x / r}[mode] # cos: mean-zero => ∮u.n=0 + s.add_rotated_freeslip_bc(datum, "Upper", normal=nhat) s.petsc_use_pressure_nullspace = True s.petsc_options["snes_type"] = "ksponly" s.tolerance = 1e-9 @@ -668,3 +665,46 @@ def test_rotated_freeslip_prescribed_normal_datum(): err = np.abs(vn - target).max() assert err < 1e-8, f"prescribed u.n=cos(theta) not imposed: max nodal error {err:.2e}" assert vn.max() > 0.9 and vn.min() < -0.9, "prescribed normal velocity magnitude wrong" + + +def test_rotated_freeslip_nonlinear_prescribed_normal_datum(): + """The prescribed wall-normal datum u.n = cos(theta) through the NONLINEAR rotated + path (power-law annulus): the loop genuinely iterates and every iterate is kept + feasible, so the converged solution carries the datum to the same (machine) + precision as the linear path — the constraint is exact independent of where the + nonlinear iteration stops.""" + RI, RO = 0.5, 1.0 + mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.2, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x**2 + y**2) + nhat = sympy.Matrix([[x / r, y / r]]) + v = uw.discretisation.MeshVariable("Vdn", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable("Pdn", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + g = sympy.Matrix([[v.sym[0].diff(x), v.sym[0].diff(y)], + [v.sym[1].diff(x), v.sym[1].diff(y)]]) + e = 0.5 * (g + g.T) + eII = sympy.sqrt(0.5 * (e[0, 0] ** 2 + e[1, 1] ** 2) + e[0, 1] ** 2 + 1.0e-12) + s.constitutive_model.Parameters.shear_viscosity_0 = eII ** (1.0 / 3.0 - 1.0) + blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2) / 0.05)) + s.bodyforce = sympy.Matrix([[50.0 * blob * x / r, 50.0 * blob * y / r]]) + s.add_essential_bc((0.0, 0.0), "Lower") # no-slip inner (pins rotation gauge) + s.add_rotated_freeslip_bc(x / r, "Upper", normal=nhat) # u.n = cos(theta), mean-zero flux + s.consistent_jacobian = True # Newton tangent (few iterations) + s.petsc_use_pressure_nullspace = True + s.tolerance = 1e-7 + s.solve() + + info = s._rotated_freeslip_info + assert info["nonlinear_iterations"] > 1, "nonlinear datum solve did not genuinely iterate" + assert info["converged"], "nonlinear datum solve did not converge" + vc = v.coords + rr = np.hypot(vc[:, 0], vc[:, 1]) + outer = np.abs(rr - RO) < 4.0e-2 # outer velocity nodes (incl. edge mids) + rhat = vc[outer] / rr[outer, None] + vn = np.einsum("ij,ij->i", v.data[outer], rhat) + target = vc[outer, 0] / rr[outer] # cos(theta) at the same node coords + err = np.abs(vn - target).max() + assert err < 1e-8, f"nonlinear u.n=cos(theta) not imposed: max nodal error {err:.2e}" + assert vn.max() > 0.9 and vn.min() < -0.9, "prescribed normal velocity magnitude wrong"