Barr & Houseman faulted-medium analytic solution in uw.analytic

[lmoresi]

The linear (n = 1) plane-strain analytic solution for a fault terminating inside a viscous medium, from the Appendix of Barr & Houseman 1996 (GJI 125, 473–490).

Implementation follows that of @gthyagi, who has been using this solution for fault benchmarking — hence the review request. A second independent implementation is the strongest check available on this, so please do compare against yours.

Why

Every fault result in the current campaign is measured against a gmsh-union ribbon control — that is, against another discretisation. There is no absolute standard. This is one, and it has a real fault in it: an internal boundary carrying zero shear traction, continuous normal velocity and continuous normal stress, with the tip inside the medium.

The structure is the interesting part

In polar coordinates about the tip, the stream function separates into a Fourier series in m = q/2:

• whole-integer m — continuous deformation
• half-integer m — the fault discontinuity itself

and boundedness of the velocity at r = 0 admits only one negative index, m = -1/2. That single mode carries the entire singularity, which is why slip goes as sqrt(r) and stress as 1/sqrt(r). The exponents are a property of the fault's own Fourier mode rather than an assumption, and one of the tests asserts exactly that — drop the half-integer term and the slip vanishes.

The tests verify it IS a solution

Not a comparison against a stored answer, and not against a UW3 solve — symbolic checks that the field satisfies the equations:

check
incompressibility	div u = 0
momentum, eta lap(u) + grad(p)	0 in both components
zero shear traction on the fault	tau_r_theta = 0 at theta = 0 and 2 pi
normal velocity continuous	u_theta equal on both faces
slip normalisation	2 U0 sqrt(r/R0), so 2 U0 at the perimeter — the paper's own anchor

If those hold simultaneously the transcription is the solution, whatever a solver later does with it.

A transcription trap, settled by mathematics

The half-integer sine terms of u_theta appear with one sign in the paper's boundary datum (A8b) and the opposite sign in its solution (A9b), and the same mismatch appears in the plane-stress pair (A12b vs A13b).

Incompressibility settles it without anyone adjudicating a scanned minus sign: for u_r = A sqrt(R) f(theta) and u_theta = sqrt(R) g(theta), div u = 0 forces g' = -(3/2) A f, which integrates to (A8b)'s sign. Verified by flipping it — the Stokes test then fails with a divergence of (0.75 cos(t/2) + 2.25 cos(3t/2))/sqrt(r), non-zero purely in the fault modes.

@gthyagi — this is the specific point worth checking against your implementation.

Scope and conventions

• Plane strain only. The paper's thin-viscous-sheet solution (A10–A13) has non-zero in-plane divergence — a different equation set from our incompressible Stokes, so not a benchmark for this solver.
• Their pressure is extension-positive, so their force balance is d_j tau_ij + d_i p = 0. The module follows the paper and documents it; negate to compare against a compression-positive solver.
• B = 2 eta_0 at n = 1 (their eq 1), so tau = 2 eta edot as usual.
• The branch cut is placed on the fault via arctan2(y, x) mod 2*pi. A bare atan2 would put it on the negative x axis and silently return the wrong face; there is a test for that.

Not yet wired to a solve

The analytic domain is a disc with the fault running from the centre to the perimeter — the fault reaches the boundary, which 2-D UW3 cannot do yet (#549). The solution and its verification stand on their own until that lands, and the benchmark can be connected immediately afterwards.

Tests: test_0210_barr_houseman_analytic (7). test_1015_analytic_solcx and test_1062_constrained_solcx pass unchanged (12 total); style gate clean.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

Adds a new analytic benchmark solution (Barr & Houseman 1996, plane strain, Newtonian n=1) to uw.function.analytic, with tests that symbolically verify the expressions satisfy incompressible Stokes and the fault interface conditions. This provides an absolute (closed-form) reference for fault benchmarking, independent of discretisation.

Changes:

• Introduces BarrHouseman analytic solution implementation (SymPy form + NumPy evaluator) and re-exports it under underworld3.function.analytic.
• Adds a focused test suite that verifies the solution by symbolic PDE/traction checks plus branch-cut/slip behaviour checks.
• Adds an import hook in the Cython analytic module to keep a single public analytic namespace.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.

File	Description
tests/test_0210_barr_houseman_analytic.py	New symbolic + numeric verification tests for the Barr–Houseman fault-tip solution.
src/underworld3/function/analytic.pyx	Re-exports BarrHouseman into the existing uw.function.analytic namespace.
src/underworld3/function/_barr_houseman.py	Implements the Barr–Houseman plane-strain (n=1) analytic solution with SymPy expressions and a NumPy evaluator.

Suppressed comments (2)

src/underworld3/function/_barr_houseman.py:211

• slip() will currently return nan for negative r without an explicit error, which can hide caller mistakes (unlike evaluate(), which validates its inputs). Adding a simple non-negativity check makes the API fail-fast and keeps downstream results from silently becoming nan.

        r = np.asarray(r, dtype=float)
        return 2.0 * self.U0 * np.sqrt(r / self.R0)

src/underworld3/function/_barr_houseman.py:188

• The closed-form expressions are duplicated in both _polar() (SymPy) and evaluate() (NumPy). Even with the cross-check test, this is a drift hazard (a future edit might update one side only). Consider extracting the shared formula into a single helper that takes {sin, cos, sqrt} callables (or lambdifying the SymPy expressions once and caching the resulting NumPy callables) so there is only one source of truth.

        u_r = (self.U0 / 4) * (
            R**2 * (np.sin(t) - np.sin(3 * t))
            - R**3 * (2 * np.sin(2 * t) - 2 * np.sin(4 * t))
            + np.sqrt(R) * (np.cos(t / 2) + 3 * np.cos(3 * t / 2))
        )

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[lmoresi]

Adversarial review

We tried to break the transcription independently rather than re-running its own tests. Two coverage gaps found, both of which we closed by independent verification — the solution itself survives.

1. The symbolic Stokes check pins U0 = R0 = eta = 1, so a parameter-scaling error in the pressure would pass. A transcription with, say, a missing R0 power in p satisfies the unit-parameter identity and fails the physics. We re-derived div u, both momentum components, with all three parameters fully symbolic: all identically zero. Recommend making the test's parameters symbolic (sympy.symbols("U0 R0 eta", positive=True)) — same runtime, strictly stronger statement.

2. The PR text claims three fault conditions; the tests assert two. Zero shear traction and normal-velocity continuity are tested; normal-stress continuity is claimed and never asserted. We checked it: with the extension-positive convention, sigma_tt = p + 2*eta*e_tt has identically zero jump across theta = 0 vs 2*pi. It holds — add the assertion so the claim in the description is load-bearing rather than decorative.

3. tier_a on a new test file. The testing doctrine reserves tier A for explicitly reviewed tests. These are deterministic symbolic identities, and with this review plus the two additions above we would endorse the marking; without them it is self-assigned.

What we could not fault: the branch-cut test is a genuine negative control (a bare atan2 fails it in exactly the advertised way); the A8b-vs-A9b sign resolution by incompressibility is the right way to settle a paper's internal inconsistency and is itself tested by its failure mode; the slip test's whole-integer-modes control discharges the vacuous-pass risk. The plane-stress exclusion (A10-A13) is correctly reasoned — different equation set, not a benchmark for this solver.

Cross-PR note: the disc-with-daylighting-fault configuration this ultimately drives needs #552 (ribbon outcrop, now green) and #553 (transversal slit) — the interior verification here stands alone, as the description says.

Underworld development team with AI support from Claude Code

[lmoresi]

Both findings accepted and fixed in 006d0ae. Thank you — the first one was a real blind spot, not a theoretical one.

1. Parameter scaling. We did not take this on faith; we built the corruption you described and measured it. Giving the singular pressure term the wrong power of R0:

	momentum residual
unit parameters (U0 = R0 = eta = 1)	0 — passes silently
symbolic parameters	U0*eta*(1 - R0)*cos(3*theta/2) / (2*sqrt(R0)*r^(3/2)) — caught

So the test was blind in exactly the way you said. BarrHouseman now accepts SymPy parameters on the symbolic paths and refuses them, naming the reason, on evaluate and slip. The Stokes and fault-condition tests run fully symbolic. Same runtime.

2. Normal-stress continuity. Correct — claimed in the description, asserted nowhere. Now asserted, and it holds as you found: with the extension-positive convention sigma_tt = 2*eta*e_tt + p has identically zero jump across theta = 0 vs 2*pi. The test module's docstring now lists all three fault conditions against the checks that cover them, so the description and the tests cannot drift apart again.

3. tier_a. Taking the endorsement as conditional on the two above, which are now in. If you would rather it sat at tier_b until @gthyagi has compared it against his implementation, say so and we will move it — the marking is not load-bearing for anything else.

On the cross-PR note: #552 and #553 are better news than we knew when this was opened. The description's "not yet wired to a solve" paragraph is now the only thing standing between this and the full manufactured benchmark, and it stops being true as soon as #552 lands and #553 is decided. We have left the paragraph as-is rather than pre-announcing a capability that has not merged.

Underworld development team with AI support from Claude Code

[gthyagi]

I independently verified the analytical solution on current head edfd64b against both Barr & Houseman papers and our BH92 benchmark analysis.

The important transcription point is slightly stronger than the PR description: BH96 equation A8b has -3 sin(theta/2) - 3 sin(3 theta/2), whereas printed A9b has both the opposite sign and cos(3 theta/2) in place of sin(3 theta/2). The A9b form used in this PR is correct:

• it reproduces A8b exactly at r = R0;
• div(u) = 0;
• both Stokes momentum residuals vanish symbolically with U0, R0, and eta left symbolic;
• zero fault shear traction, continuous normal velocity, continuous normal stress, and the published slip 2 U0 sqrt(r/R0) all hold;
• the printed A9b fails both boundary matching and incompressibility, so it is a paper typo rather than an alternative convention.

BH96 is an extension of BH92, but its Appendix disc is a separate exact n = 1 test problem rather than the same whole-domain BH92 box. Their near-tip physics is consistent: the A9 half-integer mode gives velocity proportional to r^(1/2) and strain rate/stress proportional to r^(-1/2), exactly matching the BH92 n = 1 exponents and our benchmark script.

The pressure sign in the PR also matches the papers' extension-positive convention; it must be negated only when comparing with a compression-positive numerical pressure.

I built the current branch as an isolated macOS arm64 wheel and ran an independent SymPy verification script plus test_0210_barr_houseman_analytic.py, test_1015_analytic_solcx.py, and test_1062_constrained_solcx.py.

Result: 12 passed. I have verified the analytical implementation.

[gthyagi]

I have verified only the analytical solution for BH96 so far. I’m currently building the numerical model and will compare its results with the analytical solution. Please keep this PR open until I complete the comparison.

[gthyagi]

Follow-up: I also verified this against a numerical UW3 Stokes model with an actual mesh and MeshVariables.

Because split_fault() still refuses a fault reaching the domain boundary (#549), I generated the BH96 unit-disc slit topology directly in Gmsh. The two fault faces have coincident coordinates but independent node/DOF identities, share the tip node, and meet separate copies of the outer-boundary endpoint. The negative-x seam remains welded, so this is one connected slit disc rather than two disconnected half-discs.

Numerical setup:

• P2 Velocity and continuous P1 Pressure MeshVariables, matching the BH96 Appendix discretisation;
• U0 = R0 = eta = 1;
• published A8 velocity on the outer circle;
• published common fault-normal velocity imposed on both fault faces, with tangential traction left natural (zero shear traction);
• UW3 pressure compared with -p_BH96, after removing only the constant pressure gauge;
• errors sampled on the fixed annulus 0.2 <= r/R0 <= 0.95, away from the singular tip.

h	triangles	velocity rel. L2	pressure rel. L2	slip rel. L2	RMS normal-velocity jump
0.20	202	1.920%	13.104%	1.700%	1.94e-17
0.10	780	0.744%	6.909%	0.975%	3.14e-17
0.05	3004	0.396%	3.281%	0.526%	1.30e-17

All three solves converged (SNES reason = 3), all fields were finite, and the independent coincident fault P2 DOFs reproduced the published 2 U0 sqrt(r/R0) slip with decreasing error.

The pressure sign is also confirmed numerically: using the required opposite sign gives the convergent errors above, while comparing UW3 pressure directly with the paper's extension-positive pressure gives approximately 199% error at every resolution.

Conclusion: the Mesh/MeshVariable Stokes solution converges to the PR's analytic field. This validates the analytical implementation numerically, but it does not claim that the current split_fault() boundary-intersection path works; that remains #549.