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 (GJI125, 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-integerm — continuous deformation
• half-integerm — 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 andcos(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.

[gthyagi]

Additional BH92 near-tip comparison

I compared the finest completed BH92 UW3 rectangular-fault-zone model with the BH96 half-integer singular mode implemented in this PR.

BH92 numerical model

• plane strain, n = 1, P2 velocity / continuous P1 pressure;
• rectangular finite-width fault zone (RFZ), lc = 0.01, fault half-width = 0.0005, 8 MPI ranks;
• R0 = 10 km, boundary slip u0 = 10 mm/yr, and eta = 10^20 Pa s.

BH96 normalisation

• R0 = 10 km;
• U0 = 5 mm/yr, because the BH96 perimeter slip is 2 U0 = u0;
• only the half-integer singular component was used. The regular integer modes satisfy the BH96 circular boundary data and should not match the square BH92 model globally.

For the singular-exponent comparison, I fitted along the fault extension (theta = 180 degrees) over 0.2 <= r <= 1 km. This interval lies outside the finite-width/mesh-scale tip regularisation while remaining close enough to the tip that the nonsingular square-domain field is small.

field	BH92 UW3 near-tip exponent	BH96 exact exponent
abs(eps_xy)	-0.4983	-0.5000
abs(tau_xy)	-0.4983	-0.5000

At r = 1 km, the normalized angular abs(eps_xy) and abs(tau_xy) profiles are essentially coincident with the BH96 singular solution after excluding 1.397 degrees on each side of the finite-width fault. The radial amplitudes are also close.

For clarity, a single-power-law fit over the wider 1-5 km interval gives an effective slope of -0.4521. That interval increasingly includes the regular simple-shear and outer-boundary contribution of the square BH92 problem, so it should not be interpreted as the singular exponent.

This provides an independent consistency check: the BH92 n = 1 numerical model reproduces the analytical r^(-1/2) exponent and angular structure in the resolved near-tip region. It is deliberately not presented as a full-domain equality because BH92 and the BH96 Appendix use different domains and external boundary conditions.

[gthyagi]

Updated figure for the additional BH92 near-tip comparison. Panels (a) and (b) now fit the resolved singular interval 0.2 <= r <= 1 km (shaded), giving UW3 -0.498 against the exact BH96 -0.500. Blue: finest completed BH92 UW3 RFZ model at n = 1. Dashed black: BH96 half-integer singular mode.

[lmoresi]

@gthyagi — thank you for the verification; it is the most thorough external check anything in this suite has had. The A9b finding (printed sign andcos(3θ/2) where sin(3θ/2) belongs) and your BH92 near-tip comparison at −0.498 against the exact −0.500 are now recorded in the module docstring, the docs section and the errata, with attribution.

This PR went CONFLICTING through no fault of its own: #571 landed the uw.analytic package this morning and replaced function/analytic.pyx, which this branch edits. The rebase is done and green, on branch feature/pr550-integrate — pushed separately rather than force-pushed over this PR, so nothing shifts under your comparison work. Swap it in whenever you are finished.

What the integration required:

• analytic/barr_houseman.py, class FaultedMedium, registered in available()/describe(); tests → tests/test_1029_analytic_faulted_medium.py. (Naming follows kramer.py/CylindricalStokes — author-named module, problem-named class. Trivially revertible to BarrHouseman if preferred.)
• apply_boundary_conditions REFUSES rather than inheriting a wall mixin. Both mixins loop over every labelled boundary applying one condition; the fault faces are labelled surfaces, so inheriting FixedWalls would have prescribed the exact velocity on the fault — imposing the slip instead of letting it emerge, and making the comparison against 2 U0 √(r/R0) tautological. It now raises, names 2-D fault zones cannot reach a boundary — the outcrop machinery exists but is never wired into the ribbon path #549, and hands over boundary_velocity() / fault_normal_velocity() / slip().
• A pressure-convention erratum (§3.8): BH is extension-positive, UW3 is not, so fn_pressure = −(A9c), flipped once and asserted against the printed form.
• A numerical fix worth knowing about: in Cartesian coordinates r − x cancels near the fault and θ loses accuracy as 1/θ². Scaling both atan2 arguments by (r + x) via (r−x)(r+x) = y² takes the momentum residual from 3e-6 to 3e-15 at 1e-6 rad off the fault, so sample points can sit just off each face.
• The body-force negative control cannot fire here (there is no body force), so it was replaced with one that can: flipping the pressure sign moves the momentum residual to 1.06 while tracelessness stays at 5.7e-16 — that one is structural, set_fields cancels it.
• A compiled-path test: evaluate() through the JIT agrees with the NumPy evaluator to 1.7e-15, so the Piecewise/atan2 branch cut survives code generation.

Gates on the installed package, at defaults and at U0=1.3, R0=2.0, eta=0.7, tip=(0.3,−0.2): momentum 3.6e-16 / 3.7e-15, div u 1.0e-15 / 5.6e-16, tracelessness 5.7e-16 / 4.2e-16, strain-rate 4.7e-16 / 4.2e-16.

Tests: 18 in the file (every original check preserved, including the symbolic-parameter strengthening from 006d0ae5); analytic fast tier 193 passed; tests/analytic_full/ 206 passed; scripts/test.sh --p 2 green end to end. It is cheap enough for the per-PR tier (6.3s warm) rather than the full-family tier.

Not merging — waiting on your word that the comparison is complete.

Underworld development team with AI support from Claude Code

[gthyagi]

I have updated the green integration branch in commit e4f1d9bb to rename the module from uw.analytic.barr_houseman to uw.analytic.barr_houseman_96.

The explicit path is now:

uw.analytic.barr_houseman_96.FaultedMedium

The year suffix distinguishes the closed-form Barr & Houseman (1996) solution from the related BH92 study, which does not provide the corresponding full Cartesian analytical field. The registered convenience API remains uw.analytic.FaultedMedium.

Verification on the rebuilt amr-dev package:

• qualified module path resolves to the registered class;
• tests/test_1029_analytic_faulted_medium.py: 18 passed.

This is a naming-only change; the analytical expressions and the completed BH92/BH96 comparison are unchanged.

[lmoresi]

@gthyagi — heads up on two changes to feature/pr550-integrate, since you are working on it.

Your rename is in (e4f1d9bb, barr_houseman_96.py). I merged it rather than pushing over it; nothing of yours was rewritten.

development moved under us.#578 landed: the FreeSlipWalls / FixedWalls mixins are gone, replaced by three composed functions in analytic/_base.py that each take the boundaries they apply to:

free_slip(solver, boundaries, normal=None)
prescribed_velocity(solver, boundaries, velocity)
prescribed_scalar(solver, boundaries, field)

FaultedMedium was unaffected in substance — it declares AnalyticSolution directly and its apply_boundary_conditions refuses — but its refusal message and the subsystem docs both explained the refusal by naming the two mixins. Those now say what the solution actually needs (a component-wise Dirichlet condition on an internal boundary) instead of naming classes that no longer exist. That is the only change to your work.

Worth knowing for the comparison you are running: under the mixins, a solution that refuses was an exception to the pattern. After #578 it is the ordinary shape — every solution states its own conditions, and this one's answer is "not yet these".

Branch is at db2c3749, rebuilt and green on the analytic tier (120 passed). Still holding the merge until you confirm your numerical comparison is complete.

[gthyagi]

Numerical checks completed

I rebuilt the updated integration branch at db2c3749 and reran the analytical and numerical verification.

• qualified API check: uw.analytic.barr_houseman_96.FaultedMedium resolves correctly;
• analytical test set: 329 passed, 15 skipped (all skips require the optional assess package);
• BH96 Gmsh slit-disc Stokes model, P2 velocity / continuous P1 pressure, h = 0.1:
  • serial: SNES reason 3;
  • 8 MPI ranks: SNES reason 3 and parallel HDF5/XDMF output written successfully.

Serial / 8-rank relative errors:

metric	serial	8 ranks	absolute difference
velocity L2	0.005757886311	0.005757886317	5.77e-12
pressure L2	0.055002026425	0.055002026437	1.24e-11
fault-slip L2	0.009747029060	0.009747029070	9.71e-12
normal-velocity jump RMS	3.14e-17	3.14e-17	0

For the MPI comparison, field errors were evaluated on globally gathered and coordinate-deduplicated velocity and pressure DOFs, so partition-local interpolation does not enter the result.

The updated branch preserves the analytical result and the serial and 8-rank numerical solutions are equivalent to numerical precision. My numerical comparison is complete.