From 4c8caba14638cd78a2867956ee9ce5fb37224557 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 08:45:14 +1000 Subject: [PATCH 1/5] Barr & Houseman faulted-medium analytic solution in uw.analytic MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The linear (n=1) plane-strain 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. Why it is worth having: every fault result we have is currently measured against another discretisation — a gmsh-union ribbon control — so there is no absolute standard. This is one, and it has a 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 is continuous deformation, half-integer m IS the fault discontinuity, and boundedness at r=0 admits only one negative index, m = -1/2. That single mode carries the whole 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 the test asserting the slip is carried entirely by that mode says so. The tests verify the field IS a Stokes solution rather than comparing it to a stored answer: div u = 0, the momentum balance vanishes identically in both components, zero shear traction on both faces of the fault, normal velocity continuous across it, and slip = 2 U0 sqrt(r/R0) — the paper's own normalisation. If those hold together the transcription is the solution, whatever a solver later does with it. That mattered here. 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). Incompressibility settles it: for u_r = A sqrt(R) f and u_theta = sqrt(R) g, div u = 0 forces g' = -(3/2) A f, which is (A8b)'s sign. Verified by flipping it — the Stokes test fails with a divergence of (0.75 cos(t/2) + 2.25 cos(3t/2))/sqrt(r). Plane strain only. The paper's thin-viscous-sheet solution (A10-A13) has non-zero in-plane divergence and is a different equation set from our incompressible Stokes, so it is not a benchmark for this solver. Not yet wired to a UW3 solve: the analytic domain is a disc with the fault running from the centre to the perimeter, and 2-D faults cannot reach a boundary yet (#549). The solution and its verification stand on their own until that lands. Underworld development team with AI support from Claude Code --- src/underworld3/function/_barr_houseman.py | 211 +++++++++++++++++++++ src/underworld3/function/analytic.pyx | 5 + tests/test_0210_barr_houseman_analytic.py | 170 +++++++++++++++++ 3 files changed, 386 insertions(+) create mode 100644 src/underworld3/function/_barr_houseman.py create mode 100644 tests/test_0210_barr_houseman_analytic.py diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py new file mode 100644 index 00000000..e97fad4f --- /dev/null +++ b/src/underworld3/function/_barr_houseman.py @@ -0,0 +1,211 @@ +r"""Barr & Houseman (1996) analytic solution for a fault embedded in a viscous +medium — the linear (:math:`n = 1`) plane-strain case. + +Reference +--------- +T. D. Barr & G. A. Houseman, *Deformation fields around a fault embedded in a +non-linear ductile medium*, Geophys. J. Int. **125**, 473-490 (1996); +Appendix, equations (A1)-(A9). The companion letter, Barr & Houseman, +Geophys. Res. Lett. **19**, 1145-1148 (1992), gives the near-tip asymptotics. + +Implementation follows that of @gthyagi, who has been using this solution for +fault benchmarking. + +Why this solution is unusual, and useful +---------------------------------------- +The deformation around a terminating fault is written in polar coordinates +with the **origin at the fault tip** and the fault along :math:`\theta = 0`. +The stream function separates into a Fourier series in :math:`m = q/2`, and the +two halves of that series mean different things: + +* **whole-integer** :math:`m` — continuous deformation, no fault; +* **half-integer** :math:`m` — the fault-type discontinuity. + +Boundedness of the velocity at :math:`r = 0` admits only one negative index, +:math:`m = -1/2`. That single mode carries the whole fault singularity, and it +is why the slip goes as :math:`\sqrt{r}` and the stress as +:math:`1/\sqrt{r}` — the exponents are a property of the fault's own Fourier +mode, not an assumption. + +The solution below is their plane-strain test problem: prescribe the velocity +(:meth:`boundary_velocity`) on the perimeter of a disc of radius :math:`R_0`, +impose the fault conditions on :math:`\theta = 0`, and the interior field is +exact. + +Only plane strain is implemented. The paper also gives a thin-viscous-sheet +(plane-stress) solution, equations (A10)-(A13), in which the in-plane +divergence is NOT zero — a different equation set from Underworld's +incompressible Stokes, so it is not a benchmark for this solver. + +Conventions +----------- +Their constitutive relation is :math:`\tau_{ij} = B \dot E^{(1/n - 1)} +\dot\varepsilon_{ij}` with :math:`B = 2\eta_0` at :math:`n = 1`, so +:math:`\tau = 2 \eta \dot\varepsilon` as usual. + +Their pressure takes **extension as positive**, so their force balance is +:math:`\partial_j \tau_{ij} + \partial_i p = 0` — a sign opposite to the more +common convention. :attr:`pressure` follows the paper. Negate it to compare +against a solver whose pressure is compression-positive. +""" +import numpy as np +import sympy + + +class BarrHouseman: + r"""The linear plane-strain fault-tip solution on a disc of radius ``R0``. + + The fault occupies :math:`\theta = 0` from the tip at the origin to the + perimeter. Slip is the jump in the fault-parallel velocity across it, and + is :math:`2 U_0 \sqrt{r/R_0}` — so the relative slip velocity at the + perimeter is :math:`2 U_0`, the paper's normalisation. + + Parameters + ---------- + U0 : float + Velocity scale. The slip at the perimeter is ``2 * U0``. + R0 : float + Radius of the disc, and the fault's length. + eta : float + Newtonian viscosity of the medium. + + Examples + -------- + The solution and the boundary datum that produces it: + + >>> sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + >>> float(sol.slip(1.0)) + 2.0 + + Notes + ----- + The field is multivalued around the tip — that is what a fault is — so it + is a function of :math:`(r, \theta)` with :math:`\theta \in [0, 2\pi)`, + NOT of the Cartesian coordinates alone. The branch cut lies **on the + fault**. :meth:`evaluate` places it there by taking + ``arctan2(y, x) mod 2*pi``; a Cartesian expression using a bare + ``atan2`` would put the cut on the negative :math:`x` axis instead and + silently return the wrong side of the fault. + """ + + def __init__(self, U0=1.0, R0=1.0, eta=1.0): + if not (float(R0) > 0.0 and float(eta) > 0.0): + raise ValueError("R0 and eta must be positive.") + self.U0 = float(U0) + self.R0 = float(R0) + self.eta = float(eta) + + # ------------------------------------------------------------------ sympy + @property + def symbols(self): + """The polar symbols ``(r, theta)`` the expressions are written in.""" + return sympy.symbols("r theta", positive=True) + + def _polar(self): + r, t = self.symbols + R = r / self.R0 + U0, eta, R0 = self.U0, self.eta, self.R0 + + # Whole-integer (continuous) modes, then the half-integer (fault) mode. + # The half-integer group is the entire singular content: sqrt(R) in the + # velocity, 1/sqrt(R) in the pressure. + u_r = (U0 / 4) * ( + R**2 * (sympy.sin(t) - sympy.sin(3 * t)) + - R**3 * (2 * sympy.sin(2 * t) - 2 * sympy.sin(4 * t)) + + sympy.sqrt(R) * (sympy.cos(t / 2) + 3 * sympy.cos(3 * t / 2)) + ) + u_t = (U0 / 4) * ( + R**2 * (3 * sympy.cos(t) - sympy.cos(3 * t)) + - R**3 * (4 * sympy.cos(2 * t) - 2 * sympy.cos(4 * t)) + - sympy.sqrt(R) * (3 * sympy.sin(t / 2) + 3 * sympy.sin(3 * t / 2)) + ) + p = (eta * U0 / R0) * ( + -2 * R * sympy.sin(t) + + 3 * R**2 * sympy.sin(2 * t) + + sympy.cos(t / 2) / sympy.sqrt(R) + ) + return u_r, u_t, p + + @property + def velocity_polar(self): + """``(u_r, u_theta)`` as SymPy expressions in ``r`` and ``theta``.""" + u_r, u_t, _p = self._polar() + return u_r, u_t + + @property + def pressure_polar(self): + """Pressure as a SymPy expression, EXTENSION POSITIVE (see module doc).""" + return self._polar()[2] + + def boundary_velocity(self): + r"""The velocity datum on :math:`r = R_0` that produces the solution. + + Their equations (A8a, A8b). Returned as ``(U_r, U_theta)`` SymPy + expressions in ``theta``; this is what a solver's Dirichlet condition + on the disc perimeter must impose. + """ + _r, t = self.symbols + u_r, u_t = self.velocity_polar + return (sympy.simplify(u_r.subs(_r, self.R0)), + sympy.simplify(u_t.subs(_r, self.R0))) + + # ------------------------------------------------------------------ numpy + def evaluate(self, coords): + r"""Velocity and pressure at Cartesian ``coords`` measured FROM THE TIP. + + Parameters + ---------- + coords : array_like, shape (N, 2) + Points relative to the fault tip, with the fault along ``+x``. + + Returns + ------- + velocity : ndarray, shape (N, 2) + Cartesian components. + pressure : ndarray, shape (N,) + + Notes + ----- + ``theta`` is taken as ``arctan2(y, x) mod 2*pi`` so the branch cut sits + ON the fault, which is where the field is genuinely discontinuous. A + point exactly on the fault returns the ``theta = 0`` side; approach + from ``y < 0`` to obtain the other. + """ + X = np.asarray(coords, dtype=float) + if X.ndim != 2 or X.shape[1] != 2: + raise ValueError("coords must have shape (N, 2)") + r = np.hypot(X[:, 0], X[:, 1]) + t = np.mod(np.arctan2(X[:, 1], X[:, 0]), 2.0 * np.pi) + if np.any(r == 0.0): + raise ValueError( + "the solution is singular at the fault tip; exclude r = 0") + + R = r / self.R0 + 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)) + ) + u_t = (self.U0 / 4) * ( + R**2 * (3 * np.cos(t) - np.cos(3 * t)) + - R**3 * (4 * np.cos(2 * t) - 2 * np.cos(4 * t)) + - np.sqrt(R) * (3 * np.sin(t / 2) + 3 * np.sin(3 * t / 2)) + ) + p = (self.eta * self.U0 / self.R0) * ( + -2 * R * np.sin(t) + 3 * R**2 * np.sin(2 * t) + + np.cos(t / 2) / np.sqrt(R) + ) + velocity = np.column_stack([u_r * np.cos(t) - u_t * np.sin(t), + u_r * np.sin(t) + u_t * np.cos(t)]) + return velocity, p + + def slip(self, r): + r"""Fault slip :math:`2 U_0 \sqrt{r/R_0}` at radius ``r`` from the tip. + + The jump in fault-parallel velocity between the two faces of the fault. + The :math:`\sqrt{r}` dependence is the :math:`m = -1/2` mode and is the + quantity a discrete model can be asked to reproduce — unlike the + stress, which is singular at the tip. + """ + r = np.asarray(r, dtype=float) + return 2.0 * self.U0 * np.sqrt(r / self.R0) diff --git a/src/underworld3/function/analytic.pyx b/src/underworld3/function/analytic.pyx index 93ddb1c3..876613f4 100644 --- a/src/underworld3/function/analytic.pyx +++ b/src/underworld3/function/analytic.pyx @@ -2,6 +2,11 @@ import os import sympy import underworld3 +# The Barr & Houseman faulted-medium solution is elementary (no compiled +# kernel), so it lives in a plain module and is re-exported here to keep +# one analytic namespace. +from underworld3.function._barr_houseman import BarrHouseman + # Add info for linking against the Cython compiled module which contains symbols defined below. libdir = os.path.dirname(__file__) libfile = os.path.basename(__file__) diff --git a/tests/test_0210_barr_houseman_analytic.py b/tests/test_0210_barr_houseman_analytic.py new file mode 100644 index 00000000..7ccef8f4 --- /dev/null +++ b/tests/test_0210_barr_houseman_analytic.py @@ -0,0 +1,170 @@ +"""Barr & Houseman (1996) faulted-medium solution — verified as a solution. + +These tests do not compare the expressions against a stored answer or against +a UW3 solve. They check that the field IS a Stokes solution satisfying the +fault conditions, symbolically: + + div u = 0 + eta * lap(u) + grad(p) = 0 (their eq 3; extension-positive pressure) + tau_r_theta = 0 on both faces of the fault + u_theta continuous across the fault + slip = 2 U0 sqrt(r/R0) + +If those hold simultaneously, the transcription is the solution, whatever any +solver later does with it. That is a stronger statement than a regression test +and it is what makes this usable as a benchmark. + +The transcription needed it: 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). Incompressibility settles it — 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 is +(A8b)'s sign. The test below would fail on the other choice. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function.analytic import BarrHouseman + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _polar_operators(u_r, u_t, p, r, t, eta): + """(div u, momentum_x, momentum_y) for a field given in polar form.""" + div = sympy.diff(r * u_r, r) / r + sympy.diff(u_t, t) / r + + u_x = u_r * sympy.cos(t) - u_t * sympy.sin(t) + u_y = u_r * sympy.sin(t) + u_t * sympy.cos(t) + + def lap(f): + return sympy.diff(r * sympy.diff(f, r), r) / r + sympy.diff(f, t, 2) / r**2 + + p_r, p_t = sympy.diff(p, r), sympy.diff(p, t) + grad_p_x = p_r * sympy.cos(t) - p_t * sympy.sin(t) / r + grad_p_y = p_r * sympy.sin(t) + p_t * sympy.cos(t) / r + return div, eta * lap(u_x) + grad_p_x, eta * lap(u_y) + grad_p_y + + +def test_the_field_is_an_incompressible_stokes_solution(): + """div u = 0 and the momentum balance vanishes identically.""" + sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + r, t = sol.symbols + u_r, u_t = sol.velocity_polar + p = sol.pressure_polar + + div, mom_x, mom_y = _polar_operators(u_r, u_t, p, r, t, sol.eta) + assert sympy.simplify(div) == 0 + assert sympy.simplify(mom_x) == 0 + assert sympy.simplify(mom_y) == 0 + + +def test_the_fault_conditions_hold_on_both_faces(): + """Zero shear traction, and normal velocity continuous across the fault.""" + sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + r, t = sol.symbols + u_r, u_t = sol.velocity_polar + + tau_rt = sol.eta * (r * sympy.diff(u_t / r, r) + sympy.diff(u_r, t) / r) + assert sympy.simplify(tau_rt.subs(t, 0)) == 0 + assert sympy.simplify(tau_rt.subs(t, 2 * sympy.pi)) == 0 + + # u_theta is the fault-NORMAL component and must not jump; u_r is the + # fault-parallel one and must (that jump is the slip). + assert sympy.simplify(u_t.subs(t, 0) - u_t.subs(t, 2 * sympy.pi)) == 0 + + +def test_the_slip_is_the_published_normalisation(): + """slip = 2 U0 sqrt(r/R0), so 2 U0 at the perimeter — the paper's anchor. + + Negative control: the whole-integer modes alone are continuous, so a + solution without the half-integer mode would give zero slip and pass a + weaker test vacuously. + """ + sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + r, t = sol.symbols + u_r, _u_t = sol.velocity_polar + + jump = sympy.simplify(u_r.subs(t, 0) - u_r.subs(t, 2 * sympy.pi)) + assert sympy.simplify(jump - 2 * sol.U0 * sympy.sqrt(r / sol.R0)) == 0 + assert float(jump.subs(r, sol.R0)) == pytest.approx(2.0 * sol.U0) + assert float(sol.slip(sol.R0)) == pytest.approx(2.0 * sol.U0) + + # The slip is carried entirely by the half-integer mode: drop the sqrt + # term and the fault disappears. + continuous_only = u_r - (sol.U0 / 4) * sympy.sqrt(r / sol.R0) * ( + sympy.cos(t / 2) + 3 * sympy.cos(3 * t / 2)) + assert sympy.simplify(continuous_only.subs(t, 0) + - continuous_only.subs(t, 2 * sympy.pi)) == 0 + + +def test_the_numpy_evaluator_agrees_with_the_symbolic_form(): + """The Cartesian evaluator must reproduce the polar expressions. + + Includes points either side of the fault, which is where the branch cut + lives and where a bare ``atan2`` would silently return the wrong face. + """ + sol = BarrHouseman(U0=1.3, R0=2.0, eta=0.7) + r_sym, t_sym = sol.symbols + u_r_sym, u_t_sym = sol.velocity_polar + p_sym = sol.pressure_polar + + rng = np.random.default_rng(7) + radii = rng.uniform(0.2, 1.9, 12) + angles = np.r_[rng.uniform(0.05, 2 * np.pi - 0.05, 10), 0.02, + 2 * np.pi - 0.02] + pts = np.column_stack([radii * np.cos(angles), radii * np.sin(angles)]) + + velocity, pressure = sol.evaluate(pts) + for k, (rr, tt) in enumerate(zip(radii, angles)): + subs = {r_sym: float(rr), t_sym: float(tt)} + ur = float(u_r_sym.subs(subs)) + ut = float(u_t_sym.subs(subs)) + expect = np.array([ur * np.cos(tt) - ut * np.sin(tt), + ur * np.sin(tt) + ut * np.cos(tt)]) + assert velocity[k] == pytest.approx(expect, rel=1e-10, abs=1e-12) + assert pressure[k] == pytest.approx(float(p_sym.subs(subs)), + rel=1e-10, abs=1e-12) + + +def test_the_branch_cut_lies_on_the_fault(): + """Straddling the fault must show the slip; straddling +x elsewhere must not. + + This is the test that a bare ``atan2`` fails: it would put the cut on the + negative x axis, reporting a jump where the medium is continuous and none + where the fault is. + """ + sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + eps = 1e-7 + + above, _ = sol.evaluate(np.array([[0.5, +eps]])) + below, _ = sol.evaluate(np.array([[0.5, -eps]])) + assert (above[0, 0] - below[0, 0]) == pytest.approx( + float(sol.slip(0.5)), rel=1e-4), "no slip across the fault" + + left_up, _ = sol.evaluate(np.array([[-0.5, +eps]])) + left_dn, _ = sol.evaluate(np.array([[-0.5, -eps]])) + assert np.allclose(left_up, left_dn, atol=1e-5), ( + "the medium is continuous on the fault's projection; a jump here " + "means the branch cut is in the wrong place") + + +def test_the_boundary_datum_reproduces_the_solution_at_the_perimeter(): + """The Dirichlet datum a solver would impose is the solution at r = R0.""" + sol = BarrHouseman(U0=1.0, R0=1.5, eta=1.0) + r, t = sol.symbols + u_r, u_t = sol.velocity_polar + U_r, U_t = sol.boundary_velocity() + + assert sympy.simplify(U_r - u_r.subs(r, sol.R0)) == 0 + assert sympy.simplify(U_t - u_t.subs(r, sol.R0)) == 0 + assert not U_r.free_symbols - {t}, "the datum depends on theta only" + + +def test_a_degenerate_geometry_is_refused(): + with pytest.raises(ValueError, match="positive"): + BarrHouseman(R0=0.0) + with pytest.raises(ValueError, match="positive"): + BarrHouseman(eta=-1.0) + with pytest.raises(ValueError, match="singular at the fault tip"): + BarrHouseman().evaluate(np.array([[0.0, 0.0]])) From 006d0ae50ad67fbb1a8b517da90e88e6c9af4a16 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 12:31:56 +1000 Subject: [PATCH 2/5] Close two verification gaps found in review of the analytic solution MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both from the adversarial review on #550. The symbolic Stokes and fault-condition checks pinned U0 = R0 = eta = 1, so a parameter-scaling error could pass. Demonstrated rather than assumed: give the singular pressure term the wrong power of R0 and the unit-parameter momentum residual is exactly zero, while the same transcription with the parameters free gives U0 eta (1 - R0) cos(3 theta / 2) / (2 sqrt(R0) r^(3/2)) The checks now run with the parameters symbolic, which costs nothing and is a strictly stronger statement. BarrHouseman accepts SymPy parameters for the symbolic paths and refuses them, naming the reason, on the numeric ones. The paper states three fault conditions and the tests asserted two: normal-stress continuity was claimed in the description and never checked. It holds — with the extension-positive convention, sigma_tt = 2 eta e_tt + p has identically zero jump across the fault — and is now asserted, so the claim is load-bearing. Underworld development team with AI support from Claude Code --- src/underworld3/function/_barr_houseman.py | 28 ++++++++++--- tests/test_0210_barr_houseman_analytic.py | 47 +++++++++++++++++++--- 2 files changed, 64 insertions(+), 11 deletions(-) diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py index e97fad4f..c2252b30 100644 --- a/src/underworld3/function/_barr_houseman.py +++ b/src/underworld3/function/_barr_houseman.py @@ -89,11 +89,21 @@ class BarrHouseman: """ def __init__(self, U0=1.0, R0=1.0, eta=1.0): - if not (float(R0) > 0.0 and float(eta) > 0.0): - raise ValueError("R0 and eta must be positive.") - self.U0 = float(U0) - self.R0 = float(R0) - self.eta = float(eta) + # The parameters may be SymPy symbols. That is not a convenience: a + # symbolic check with U0 = R0 = eta = 1 would be satisfied by a + # transcription carrying the wrong power of R0 in the pressure, so the + # verification is strictly stronger with them left free. + for name, value in (("R0", R0), ("eta", eta)): + if not isinstance(value, sympy.Basic) and not float(value) > 0.0: + raise ValueError(f"{name} must be positive.") + self.U0 = U0 if isinstance(U0, sympy.Basic) else float(U0) + self.R0 = R0 if isinstance(R0, sympy.Basic) else float(R0) + self.eta = eta if isinstance(eta, sympy.Basic) else float(eta) + + @property + def _is_symbolic(self): + return any(isinstance(v, sympy.Basic) + for v in (self.U0, self.R0, self.eta)) # ------------------------------------------------------------------ sympy @property @@ -171,6 +181,10 @@ def evaluate(self, coords): point exactly on the fault returns the ``theta = 0`` side; approach from ``y < 0`` to obtain the other. """ + if self._is_symbolic: + raise ValueError( + "this solution was built with symbolic parameters; give U0, " + "R0 and eta numeric values to evaluate it") X = np.asarray(coords, dtype=float) if X.ndim != 2 or X.shape[1] != 2: raise ValueError("coords must have shape (N, 2)") @@ -207,5 +221,9 @@ def slip(self, r): quantity a discrete model can be asked to reproduce — unlike the stress, which is singular at the tip. """ + if self._is_symbolic: + raise ValueError( + "this solution was built with symbolic parameters; give U0, " + "R0 and eta numeric values to evaluate it") r = np.asarray(r, dtype=float) return 2.0 * self.U0 * np.sqrt(r / self.R0) diff --git a/tests/test_0210_barr_houseman_analytic.py b/tests/test_0210_barr_houseman_analytic.py index 7ccef8f4..47cdb794 100644 --- a/tests/test_0210_barr_houseman_analytic.py +++ b/tests/test_0210_barr_houseman_analytic.py @@ -6,10 +6,17 @@ div u = 0 eta * lap(u) + grad(p) = 0 (their eq 3; extension-positive pressure) - tau_r_theta = 0 on both faces of the fault - u_theta continuous across the fault + tau_r_theta = 0 on both faces of the fault fault condition 3 + u_theta continuous across the fault fault condition 2 + sigma_theta_theta continuous across the fault fault condition 1 slip = 2 U0 sqrt(r/R0) +The Stokes and fault-condition checks run with the parameters left SYMBOLIC. +That is not decoration: with U0 = R0 = eta = 1, a transcription carrying the +wrong power of R0 in the singular pressure term gives a momentum residual of +exactly zero, and the same transcription with the parameters free gives +U0 eta (1 - R0) cos(3 theta / 2) / (2 sqrt(R0) r^(3/2)). Measured. + If those hold simultaneously, the transcription is the solution, whatever any solver later does with it. That is a stronger statement than a regression test and it is what makes this usable as a benchmark. @@ -47,8 +54,15 @@ def lap(f): def test_the_field_is_an_incompressible_stokes_solution(): - """div u = 0 and the momentum balance vanishes identically.""" - sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + """div u = 0 and the momentum balance vanishes identically. + + The parameters are left SYMBOLIC on purpose. With U0 = R0 = eta = 1 a + transcription carrying the wrong power of R0 in the pressure satisfies the + identity and still fails the physics, so the unit-parameter version of this + test is strictly weaker for the same runtime. + """ + U0, R0, eta = sympy.symbols("U_0 R_0 eta", positive=True) + sol = BarrHouseman(U0=U0, R0=R0, eta=eta) r, t = sol.symbols u_r, u_t = sol.velocity_polar p = sol.pressure_polar @@ -60,8 +74,13 @@ def test_the_field_is_an_incompressible_stokes_solution(): def test_the_fault_conditions_hold_on_both_faces(): - """Zero shear traction, and normal velocity continuous across the fault.""" - sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + """All THREE of the paper's fault conditions, with symbolic parameters. + + Zero shear traction, continuous normal velocity, continuous normal stress. + The third was claimed in the original description and not asserted. + """ + U0, R0, eta = sympy.symbols("U_0 R_0 eta", positive=True) + sol = BarrHouseman(U0=U0, R0=R0, eta=eta) r, t = sol.symbols u_r, u_t = sol.velocity_polar @@ -73,6 +92,16 @@ def test_the_fault_conditions_hold_on_both_faces(): # fault-parallel one and must (that jump is the slip). assert sympy.simplify(u_t.subs(t, 0) - u_t.subs(t, 2 * sympy.pi)) == 0 + # Normal STRESS continuity — the third of the paper's three fault + # conditions. Extension-positive, so sigma = tau + p I, and the + # fault-normal component is sigma_tt = 2 eta e_tt + p with + # e_tt = (1/r) du_theta/dtheta + u_r/r. + p = sol.pressure_polar + e_tt = sympy.diff(u_t, t) / r + u_r / r + sigma_tt = 2 * sol.eta * e_tt + p + assert sympy.simplify(sigma_tt.subs(t, 0) + - sigma_tt.subs(t, 2 * sympy.pi)) == 0 + def test_the_slip_is_the_published_normalisation(): """slip = 2 U0 sqrt(r/R0), so 2 U0 at the perimeter — the paper's anchor. @@ -168,3 +197,9 @@ def test_a_degenerate_geometry_is_refused(): BarrHouseman(eta=-1.0) with pytest.raises(ValueError, match="singular at the fault tip"): BarrHouseman().evaluate(np.array([[0.0, 0.0]])) + + symbolic = BarrHouseman(U0=sympy.Symbol("U_0", positive=True)) + with pytest.raises(ValueError, match="symbolic parameters"): + symbolic.evaluate(np.array([[0.5, 0.1]])) + with pytest.raises(ValueError, match="symbolic parameters"): + symbolic.slip(0.5) From edfd64bef8ee29afcdf2c65dcdf4ed924527fcc1 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 13:49:51 +1000 Subject: [PATCH 3/5] Cache the polar symbols, and let theta be real rather than positive MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Review finding on #550. `theta` carried `positive=True` while the fault conditions are checked at `theta = 0`, which that assumption excludes — SymPy would be within its rights to simplify a substitution the assumption says cannot occur. It is now `real=True`; `r` stays positive, which is honest since the solution is singular at the origin. The symbols were also rebuilt on every property access, so identity across calls depended on SymPy's global symbol cache rather than on us. Cached on the instance instead. Underworld development team with AI support from Claude Code --- src/underworld3/function/_barr_houseman.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py index c2252b30..0826154d 100644 --- a/src/underworld3/function/_barr_houseman.py +++ b/src/underworld3/function/_barr_houseman.py @@ -99,6 +99,7 @@ def __init__(self, U0=1.0, R0=1.0, eta=1.0): self.U0 = U0 if isinstance(U0, sympy.Basic) else float(U0) self.R0 = R0 if isinstance(R0, sympy.Basic) else float(R0) self.eta = eta if isinstance(eta, sympy.Basic) else float(eta) + self._symbols = None @property def _is_symbolic(self): @@ -108,8 +109,22 @@ def _is_symbolic(self): # ------------------------------------------------------------------ sympy @property def symbols(self): - """The polar symbols ``(r, theta)`` the expressions are written in.""" - return sympy.symbols("r theta", positive=True) + """The polar symbols ``(r, theta)`` the expressions are written in. + + ``r`` is positive — the solution is singular at ``r = 0`` and never + evaluated there — but ``theta`` is only REAL. It runs over + :math:`[0, 2\pi)` and the fault conditions are checked at + :math:`\theta = 0`, which a ``positive=True`` assumption excludes; + SymPy would then be entitled to simplify a substitution that the + assumption says cannot happen. + + Cached on the instance, so repeated access returns the same objects + rather than relying on SymPy's global symbol cache for identity. + """ + if self._symbols is None: + self._symbols = (sympy.Symbol("r", positive=True), + sympy.Symbol("theta", real=True)) + return self._symbols def _polar(self): r, t = self.symbols From 18491a28ea5062dec5ae26404c5dd6f7496dc2f4 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 16:16:24 +1000 Subject: [PATCH 4/5] Split velocity and pressure evaluation: the tip is singular in one, not both Every term of the velocity carries a positive power of r, so the limit at the fault tip exists and is zero. Refusing r = 0 for the velocity was over-strict, and it bit immediately: a mesh node lands exactly on the tip in any benchmark that puts the tip inside the domain, which is the whole point of the solution. evaluate_velocity now accepts the tip and returns zero there; evaluate_pressure refuses it and says why, since the pressure carries the r^(-1/2) term of the m = -1/2 mode and genuinely diverges. evaluate() returns both and inherits the pressure's refusal. Underworld development team with AI support from Claude Code --- src/underworld3/function/_barr_houseman.py | 56 +++++++++++++++------- tests/test_0210_barr_houseman_analytic.py | 6 ++- 2 files changed, 44 insertions(+), 18 deletions(-) diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py index 0826154d..e65c65a7 100644 --- a/src/underworld3/function/_barr_houseman.py +++ b/src/underworld3/function/_barr_houseman.py @@ -196,20 +196,18 @@ def evaluate(self, coords): point exactly on the fault returns the ``theta = 0`` side; approach from ``y < 0`` to obtain the other. """ - if self._is_symbolic: - raise ValueError( - "this solution was built with symbolic parameters; give U0, " - "R0 and eta numeric values to evaluate it") - X = np.asarray(coords, dtype=float) - if X.ndim != 2 or X.shape[1] != 2: - raise ValueError("coords must have shape (N, 2)") - r = np.hypot(X[:, 0], X[:, 1]) - t = np.mod(np.arctan2(X[:, 1], X[:, 0]), 2.0 * np.pi) - if np.any(r == 0.0): - raise ValueError( - "the solution is singular at the fault tip; exclude r = 0") + velocity = self.evaluate_velocity(coords) + return velocity, self.evaluate_pressure(coords) - R = r / self.R0 + def evaluate_velocity(self, coords): + r"""Velocity at Cartesian ``coords`` measured from the tip. + + Defined AT the tip: every term of the velocity carries a positive + power of :math:`r`, so the limit is zero and is returned. It is the + pressure that is singular there, not the velocity — see + :meth:`evaluate_pressure`. + """ + r, t, R = self._polar_of(coords) 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)) @@ -220,13 +218,37 @@ def evaluate(self, coords): - R**3 * (4 * np.cos(2 * t) - 2 * np.cos(4 * t)) - np.sqrt(R) * (3 * np.sin(t / 2) + 3 * np.sin(3 * t / 2)) ) - p = (self.eta * self.U0 / self.R0) * ( + return np.column_stack([u_r * np.cos(t) - u_t * np.sin(t), + u_r * np.sin(t) + u_t * np.cos(t)]) + + def evaluate_pressure(self, coords): + r"""Pressure at Cartesian ``coords``; refuses the tip. + + The pressure carries the :math:`r^{-1/2}` term of the + :math:`m = -1/2` mode and genuinely diverges at :math:`r = 0`. + """ + r, t, R = self._polar_of(coords) + if np.any(r == 0.0): + raise ValueError( + "the pressure is singular at the fault tip; exclude r = 0 " + "(the velocity is defined there — use evaluate_velocity)") + return (self.eta * self.U0 / self.R0) * ( -2 * R * np.sin(t) + 3 * R**2 * np.sin(2 * t) + np.cos(t / 2) / np.sqrt(R) ) - velocity = np.column_stack([u_r * np.cos(t) - u_t * np.sin(t), - u_r * np.sin(t) + u_t * np.cos(t)]) - return velocity, p + + def _polar_of(self, coords): + """(r, theta, r/R0) from Cartesian coordinates, cut ON the fault.""" + if self._is_symbolic: + raise ValueError( + "this solution was built with symbolic parameters; give U0, " + "R0 and eta numeric values to evaluate it") + X = np.asarray(coords, dtype=float) + if X.ndim != 2 or X.shape[1] != 2: + raise ValueError("coords must have shape (N, 2)") + r = np.hypot(X[:, 0], X[:, 1]) + t = np.mod(np.arctan2(X[:, 1], X[:, 0]), 2.0 * np.pi) + return r, t, r / self.R0 def slip(self, r): r"""Fault slip :math:`2 U_0 \sqrt{r/R_0}` at radius ``r`` from the tip. diff --git a/tests/test_0210_barr_houseman_analytic.py b/tests/test_0210_barr_houseman_analytic.py index 47cdb794..1d657272 100644 --- a/tests/test_0210_barr_houseman_analytic.py +++ b/tests/test_0210_barr_houseman_analytic.py @@ -195,8 +195,12 @@ def test_a_degenerate_geometry_is_refused(): BarrHouseman(R0=0.0) with pytest.raises(ValueError, match="positive"): BarrHouseman(eta=-1.0) + # The PRESSURE is singular at the tip and refuses; the VELOCITY is + # defined there — every term carries a positive power of r — and is zero. with pytest.raises(ValueError, match="singular at the fault tip"): - BarrHouseman().evaluate(np.array([[0.0, 0.0]])) + BarrHouseman().evaluate_pressure(np.array([[0.0, 0.0]])) + assert np.allclose( + BarrHouseman().evaluate_velocity(np.array([[0.0, 0.0]])), 0.0) symbolic = BarrHouseman(U0=sympy.Symbol("U_0", positive=True)) with pytest.raises(ValueError, match="symbolic parameters"): From 18f36a205aa81efd511849b75e29755c8195129e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 08:38:59 +1000 Subject: [PATCH 5/5] Add the analytic traction, so a boundary can leave a component free MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Prompted by Louis: a benchmark boundary does not have to be pinned everywhere, and leaving the normal component free on one wall is worth doing. With velocity Dirichlet on EVERY wall the pressure is determined only up to a constant, and the datum must additionally carry exactly zero net flux; a traction condition on one component removes both requirements rather than patching them. It is also what the paper does — Barr & Houseman's left-hand boundary carries a constant normal stress, not a prescribed normal velocity. evaluate_traction returns sigma . n with sigma = tau + p I, extension positive as the paper has it, refusing the tip where the stress diverges. The stress is built from SymPy-derived strain rates, lambdified once and rotated into Cartesian. The test is a cross-check rather than a restatement: it arrives at the zero-shear fault condition by this different route, with a negative control off the fault. Underworld development team with AI support from Claude Code --- src/underworld3/function/_barr_houseman.py | 48 ++++++++++++++++++++++ tests/test_0210_barr_houseman_analytic.py | 29 +++++++++++++ 2 files changed, 77 insertions(+) diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py index e65c65a7..3482c551 100644 --- a/src/underworld3/function/_barr_houseman.py +++ b/src/underworld3/function/_barr_houseman.py @@ -100,6 +100,7 @@ def __init__(self, U0=1.0, R0=1.0, eta=1.0): self.R0 = R0 if isinstance(R0, sympy.Basic) else float(R0) self.eta = eta if isinstance(eta, sympy.Basic) else float(eta) self._symbols = None + self._stress_fn = None @property def _is_symbolic(self): @@ -250,6 +251,53 @@ def _polar_of(self, coords): t = np.mod(np.arctan2(X[:, 1], X[:, 0]), 2.0 * np.pi) return r, t, r / self.R0 + def evaluate_traction(self, coords, normal): + r"""Traction :math:`\sigma \cdot \hat n` at Cartesian ``coords``. + + Extension-positive, matching the paper: :math:`\sigma = \tau + p I` + with :math:`\tau = 2\eta\dot\varepsilon`. Refuses the tip, where the + stress diverges. + + This is what a boundary needs if its NORMAL velocity component is left + free rather than prescribed. Leaving one component free is worth doing: + with velocity Dirichlet on every wall the pressure is determined only + up to a constant AND the datum must carry exactly zero net flux, and a + traction condition removes both requirements at once. It is also what + Barr & Houseman do — their left-hand boundary carries a constant normal + stress, not a prescribed normal velocity. + """ + r, t, _R = self._polar_of(coords) + if np.any(r == 0.0): + raise ValueError( + "the stress is singular at the fault tip; exclude r = 0") + srr, srt, stt = self._stress_polar_numeric(r, t) + + n = np.asarray(normal, dtype=float) + if n.ndim == 1: + n = np.broadcast_to(n, (len(r), 2)) + c, s_ = np.cos(t), np.sin(t) + # rotate the polar stress into Cartesian, then contract with n + sxx = srr * c**2 - 2 * srt * c * s_ + stt * s_**2 + sxy = (srr - stt) * c * s_ + srt * (c**2 - s_**2) + syy = srr * s_**2 + 2 * srt * c * s_ + stt * c**2 + return np.column_stack([sxx * n[:, 0] + sxy * n[:, 1], + sxy * n[:, 0] + syy * n[:, 1]]) + + def _stress_polar_numeric(self, r, t): + """(sigma_rr, sigma_r_theta, sigma_theta_theta), built once via SymPy.""" + if self._stress_fn is None: + rs, ts = self.symbols + u_r, u_t = self.velocity_polar + p = self.pressure_polar + e_rr = sympy.diff(u_r, rs) + e_tt = sympy.diff(u_t, ts) / rs + u_r / rs + e_rt = (rs * sympy.diff(u_t / rs, rs) + sympy.diff(u_r, ts) / rs) / 2 + two_eta = 2 * self.eta + self._stress_fn = sympy.lambdify( + (rs, ts), [two_eta * e_rr + p, two_eta * e_rt, + two_eta * e_tt + p], "numpy") + return self._stress_fn(r, t) + def slip(self, r): r"""Fault slip :math:`2 U_0 \sqrt{r/R_0}` at radius ``r`` from the tip. diff --git a/tests/test_0210_barr_houseman_analytic.py b/tests/test_0210_barr_houseman_analytic.py index 1d657272..46be5082 100644 --- a/tests/test_0210_barr_houseman_analytic.py +++ b/tests/test_0210_barr_houseman_analytic.py @@ -207,3 +207,32 @@ def test_a_degenerate_geometry_is_refused(): symbolic.evaluate(np.array([[0.5, 0.1]])) with pytest.raises(ValueError, match="symbolic parameters"): symbolic.slip(0.5) + + +def test_the_traction_reproduces_the_zero_shear_fault_condition(): + """The traction machinery must independently give tau_r_theta = 0 on the fault. + + `evaluate_traction` builds the stress by a different route from the + symbolic fault-condition test — SymPy-derived strain rates lambdified and + rotated into Cartesian — so agreeing with it is a genuine cross-check + rather than a restatement. + + On the fault the outward normal of the upper face is -theta_hat, i.e. + (0, -1) in Cartesian along theta = 0. The SHEAR part of that traction is + its x-component, and it is the quantity the paper sets to zero. + """ + sol = BarrHouseman(U0=1.0, R0=1.0, eta=1.0) + x = np.array([0.1, 0.25, 0.4, 0.6]) + on_fault = np.column_stack([x, np.zeros_like(x)]) + + traction = sol.evaluate_traction(on_fault, [0.0, -1.0]) + assert np.allclose(traction[:, 0], 0.0, atol=1e-10), ( + f"shear traction on the fault is {traction[:, 0]}, not zero") + + # Negative control: off the fault it is emphatically NOT zero, so the + # assertion above is not passing for a trivial reason. + off_fault = np.column_stack([x, np.full_like(x, 0.15)]) + assert np.abs(sol.evaluate_traction(off_fault, [0.0, -1.0])[:, 0]).max() > 0.1 + + with pytest.raises(ValueError, match="singular at the fault tip"): + sol.evaluate_traction(np.array([[0.0, 0.0]]), [1.0, 0.0])