Analytic solutions declare their own boundary conditions

[lmoresi]

Analytic solutions declare their own boundary conditions

The ruling, and why

uw.analytic shipped with two boundary-condition mixins:

classFreeSlipWalls:
defapply_boundary_conditions(self, solver):
forboundaryinself.boundaries:
solver.add_rotated_freeslip_bc(0.0, boundary)
solver.petsc_use_pressure_nullspace=True

Each loops over every boundary of the solution and applies one condition
to all of them. That encodes "every boundary is a wall of the same kind". It is
true of the classical box benchmarks the suite started with, and false of most of
what it has to serve next — a spherical shell or annulus with different
conditions on the two radii, a faulted disc, a channel driven at one end.
FreeSlipWalls would have had to grow a parameter for each of those geometries,
and inheritance advertises the choice as if it were part of what the solution
is rather than something it decides.

The ruling: delete the mixins, keep the shared bodies as functions.

What replaced them

Every solution now writes its own apply_boundary_conditions and composes three
module-level helpers in analytic/_base.py. Each takes the boundaries it
applies to — an explicit list, not self.boundaries:

free_slip(solver, boundaries, normal=None) # strong rotated u.n = 0prescribed_velocity(solver, boundaries, velocity) # Dirichlet velocityprescribed_scalar(solver, boundaries, field) # Dirichlet scalar

A Velic solution is then three lines, and the composition is a choice the reader
can see rather than something a base class does on its behalf:

defapply_boundary_conditions(self, solver):
"""Free slip on all four walls; the enclosed box has a pressure nullspace."""free_slip(solver, self.boundaries)
solver.petsc_use_pressure_nullspace=True

Per-boundary granularity is the whole point, and it now costs nothing:

defapply_boundary_conditions(self, solver):
free_slip(solver, ["Upper"])
prescribed_velocity(solver, ["Lower"], self.fn_velocity)

Three further things fall out of this:

The pressure nullspace is now a statement, not a side effect. It is a
property of the domain — enclosed, so the pressure is determined only up to a
constant — and not of any one boundary's condition. Two different mixins both
setting it hid that. Each solution now says
solver.petsc_use_pressure_nullspace = True for itself, one line, in view.

The third copy of the Dirichlet loop is gone._Transport and _Gardner
each carried their own scalar version because neither Stokes mixin fitted a
scalar solution; both now call prescribed_scalar.

free_slip takes normal=. A curved-geometry solution can now ask for an
analytic normal. CylindricalStokes — the one curved case in the suite — keeps
the geometric default, and that is now written down as deliberate rather than
left as an omission: per #561 the geometric normal is measure-weighted to match
the straight-facet integral the assembler evaluates, so the constant pressure
stays a null vector to machine precision, where X/|X| is exact for the true
circle but keeps a consistency error that grows with facet non-uniformity. See
"Which normal to use" in docs/developer/subsystems/rotated-freeslip.md.

The equivalence evidence

This is a pure refactor, and it was verified rather than asserted.

A recording mock solver stands in for the solver and captures the ordered
sequence of method calls and attribute sets that apply_boundary_conditions
produces — sympy arguments compared by srepr, floats by their exact hex form,
so the recording is a fingerprint and not a pretty-printed approximation. It was
run over every registered solution on the merge base and again on this branch:

• 21 cases: all 20 registered solutions (assess present, so
  CylindricalStokes is covered) plus both of its boundary="free" /
  boundary="zero" branches, which are the only case in the suite where the
  conditions depend on a constructor argument.
• diff of the two recordings: empty. Byte-identical, every solution.

What the equivalence check exposed

One latent defect, left in place, filed as #577 and marked with a TODO(BUG)
at the exact location (analytic/kramer.py,
CylindricalStokes.apply_boundary_conditions):

CylindricalStokes with boundary="zero" returns without removing the pressure
nullspace, where every other enclosed case in the suite removes it. An annulus
held at zero velocity on both arcs is enclosed, so its pressure is determined
only up to a constant, and a direct solve on the singular saddle can return a
quiet, wrong answer — which is precisely the failure this suite exists to catch.
It is preserved here because this PR is behaviour-preserving by contract; the fix
changes an answer and belongs under its own regression test.

The old code made this easy to miss: the nullspace was set inside the free-slip
branch, so the zero-slip branch's return skipped it. With the nullspace stated
by the solution rather than buried in a wall type, the omission is visible on the
page.

No deprecation is owed

__init__.py re-exported FreeSlipWalls and FixedWalls; both are removed
outright, with no shim. The package landed on 2026-08-15 as #571, so there are
no external users to break. free_slip, prescribed_velocity and
prescribed_scalar take their place in __all__ — a solution written outside
this package has to be able to reach the helpers it composes, so they are public
API rather than a private convenience (Style Charter §6, no deep-import-only
features).

Composes with the Barr & Houseman work

feature/pr550-integrate (#550, not yet merged) has an apply_boundary_conditions
that refuses, with a message naming #549. Under the mixins that was an
exception to the pattern; after this refactor it is the ordinary shape — a
solution stating its own boundary conditions, which in that case is "not yet
these". It composes cleanly and needs only a trivial rebase (its class declaration
does not name either mixin).

Tests

• tests/test_1016_analytic_contract.py — the two mixin tests are rewritten
  against the helpers, and three tests are added: that two boundaries can carry
  different conditions (the reason for the change), that normal= reaches
  the solver when given and is absent when not, and that the helpers are exported.
• Analytic fast tier (test_1015–test_1028): 309 passed, 4m13s.
• tests/analytic_full/ — the whole family, every gate: 189 passed, 9m11s.
• scripts/test.sh --p 2 end to end, serial batches plus np=2: 2163 passed,
  0 failures, 0 errors, exit 0.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review

We went looking for behaviour this refactor changes without saying so, and for places where the new shape is weaker than the one it replaces.

The equivalence claim holds, and it is the right claim

The recording mock captures the ordered call sequence with sympy compared by srepr and floats by float.hex, which is a fingerprint rather than a pretty-print. Both recordings are 75038 bytes and diff is empty over 21 cases. That covers the only thing the mixins supplied.

Worth stating explicitly, because deleting a base class is not usually free: both mixins carried exactly one method each and no attributes, so removing them from the MRO cannot change anything except apply_boundary_conditions. We checked the class bodies rather than inferring it from the test result.

Findings

1. prescribed_velocity's docstring is narrower than what it accepts. It documents velocity : sympy Matrix, shape (1, dim), and CylindricalStokes calls it with a plain tuple (0.0, 0.0). That call is correct — add_dirichlet_bc takes a sequence — but the one in-tree caller that is not a Matrix is already outside the documented type. Either widen the docstring or the next author writing a solution with a constant wall velocity will build a Matrix they did not need.

2. test_free_slip_passes_an_analytic_normal_through asserts against the solver's default, not the helper's behaviour.assert registered["Left"] is None passes because add_rotated_freeslip_bc's own default happens to be None. The claim the test is making — "we did not pass normal when none was chosen" — is real and worth pinning, but as written it fails if the solver's default ever changes, in a test file about the analytic contract. Asserting on the recorded call instead of the stored value would test the helper.

3. The asymmetry between the two Dirichlet helpers is unexplained in the code.prescribed_scalar wraps its argument (add_dirichlet_bc([field], ...)) and prescribed_velocity does not. Both are right — the scalar solvers want a one-component sequence — but the two sit six lines apart and read like an inconsistency. One line saying why would stop someone "fixing" it.

None of these change a result.

The defect the check exposed

CylindricalStokes(boundary="zero") skipping the pressure nullspace is filed as #577 and left in place, correctly: this PR is byte-equivalent by contract and the fix changes an answer.

We note what made it findable. Under the mixins the nullspace was set inside a wall type, so the zero-slip branch's early return skipped it silently and no reading of apply_boundary_conditions would have shown the omission. It is visible on the page now because the nullspace is stated by the solution. That is a small piece of evidence for the design change beyond the argument from geometry.

Coverage

scripts/test.sh --p 2 end to end (2163 passed, 0 failed, exit 0), the analytic fast tier, and tests/analytic_full/. Given how many defects this session traced to tests that were never wired in, running the full script rather than the touched files is the right level.

pixi run docs-build was not run. The change is markdown-only inside an existing page with no new toctree entry; we confirmed the one new cross-reference resolves (docs/developer/subsystems/rotated-freeslip.md, section "Which normal to use") and that no MyST directives were added.

[lmoresi]

All three review findings are addressed in df1b7b2.

Finding 2 was the one worth acting on. test_free_slip_passes_an_analytic_normal_through now records the call rather than reading the value back off the solver:

assertsolver.calls== [
(0.0, "Left", {}),
(0.0, "Right", {"normal": radial}),
]

That is the claim the helper actually makes — not passing normal is different from passing normal=None, and only the recorded call distinguishes them. A separate one-line test confirms the recorded call is one a real Stokes accepts, so the mock cannot drift away from the interface it stands in for.

Findings 1 and 3 are docstrings: prescribed_velocity documents a dim-long sequence rather than a Matrix, and prescribed_scalar says why it brackets its argument.

Equivalence re-verified after the edits — the recording is byte-identical to the merge-base baseline, 119 lines, diff empty. test_1016: 23 passed.