From 4c8caba14638cd78a2867956ee9ce5fb37224557 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 14 Aug 2026 08:45:14 +1000 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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]) From e035377132cb2b81071c8c992b4bc7c4bb7bc525 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 23:26:38 +1000 Subject: [PATCH 06/11] Barr & Houseman's faulted medium joins uw.analytic MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Rebase of PR #550 (commits 4c8caba1..18f36a20) onto the package that landed in #571. The mathematics is @gthyagi's implementation, unchanged and not re-derived; what moved is where it lives and which conventions it is stated in. function/_barr_houseman.py -> analytic/barr_houseman.py class BarrHouseman -> class FaultedMedium, registered The naming follows kramer.py, the package's other single-paper module: module named for the authors, class named for the problem. The old import path needs no shim — BarrHouseman was never in the released compiled extension, so #571's shim correctly does not forward it. It is a full member of the contract: fields built through set_fields rather than around it, symbolic on mesh.X, swept by test_1024 and test_1028. Three ways it does not fit, stated rather than bent: Not a function of position. The field is multivalued about the tip, so the class carries two representations of one solution — the polar fault-frame expressions, where the fault conditions can be stated exactly at theta = 0 and 2 pi, and the contract's fn_* on mesh.X. A test pins the second to the first. The branch cut has to lie ON the fault; from tan(theta/2) = (r-x)/y that is 2*atan2(r-x, y), but near the fault r - x cancels and theta loses accuracy as 1/theta^2. Since (r-x)(r+x) = y^2, scaling both atan2 arguments by the positive (r+x) leaves the angle alone and removes it: measured at 1e-6 radians off the fault, the momentum residual is 3e-6 formed directly and 3e-15 formed this way, so sample_points can sit just off each face. No walls. apply_boundary_conditions refuses, naming the reason and handing over boundary_velocity(), fault_normal_velocity() and slip(). The fault is an internal boundary with two coincident faces (#549) and its conditions are per-component; applying the perimeter datum and leaving the fault unconstrained would solve a different problem and report a plausible error. No body force. Boundary-driven, like EllipticalInclusion, so the family's negative control cannot fire and it is excluded by name in both sweeps with the claim asserted. For this solution the momentum residual certifies the pressure sign instead: 3.6e-16 with UW3's compression-positive pressure, 1.06 with the paper's extension-positive one, while tr(sigma) + d p sits at 5.7e-16 either way. Their pressure is extension-positive and ours is not, so fn_pressure is minus their (A9c) — flipped once, and asserted against the printed form rather than left to a comment. Recorded as an erratum with the two defects @gthyagi found in printed A9b. Per-PR tier: 23s for the test file, about 6s across the two sweeps. Underworld development team with AI support from Claude Code --- src/underworld3/analytic/__init__.py | 3 + src/underworld3/analytic/barr_houseman.py | 609 ++++++++++++++++++ .../test_analytic_full_family.py | 22 +- tests/test_1028_analytic_parameter_sweep.py | 28 +- tests/test_1029_analytic_faulted_medium.py | 502 +++++++++++++++ 5 files changed, 1151 insertions(+), 13 deletions(-) create mode 100644 src/underworld3/analytic/barr_houseman.py create mode 100644 tests/test_1029_analytic_faulted_medium.py diff --git a/src/underworld3/analytic/__init__.py b/src/underworld3/analytic/__init__.py index f7a45858..0b9c7655 100644 --- a/src/underworld3/analytic/__init__.py +++ b/src/underworld3/analytic/__init__.py @@ -28,6 +28,7 @@ from ._base import AnalyticSolution, FreeSlipWalls, FixedWalls +from .barr_houseman import FaultedMedium from .inclusion import EllipticalInclusion from .kramer import CylindricalStokes from .richards import GardnerSteady, GardnerTransient @@ -41,6 +42,7 @@ "AdvectedFront", "EllipticalInclusion", "ErfcDiffusion", + "FaultedMedium", "CylindricalStokes", "GardnerSteady", "GardnerTransient", @@ -72,6 +74,7 @@ "CylindricalStokes": CylindricalStokes, "EllipticalInclusion": EllipticalInclusion, "ErfcDiffusion": ErfcDiffusion, + "FaultedMedium": FaultedMedium, "GardnerSteady": GardnerSteady, "GardnerTransient": GardnerTransient, "Poisson1D": Poisson1D, diff --git a/src/underworld3/analytic/barr_houseman.py b/src/underworld3/analytic/barr_houseman.py new file mode 100644 index 00000000..1ea650db --- /dev/null +++ b/src/underworld3/analytic/barr_houseman.py @@ -0,0 +1,609 @@ +r"""Barr & Houseman (1996): flow around a fault that terminates inside a medium. + +The linear (:math:`n = 1`) plane-strain case of their Appendix. + +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 uses this solution for fault +benchmarking, and who verified it independently against both papers and against +two numerical models — see the module notes below and PR #550. + +Why this solution is unusual, and useful +---------------------------------------- +Every other solution in this suite is a smooth field on a box. This one has a +**fault** in it: an internal boundary carrying zero shear traction, with +continuous normal velocity and continuous normal stress, and with its tip inside +the medium. That makes it the only absolute standard here for a fault +calculation — everything else a fault model can be measured against is another +discretisation. + +The deformation 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:`FaultedMedium.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. + +Two transcription points, both settled by mathematics +----------------------------------------------------- +**The sign of the half-integer sine terms.** They appear with one sign in the +paper's boundary datum (A8b) and the opposite sign in its solution (A9b). +Incompressibility settles it without anyone adjudicating a scanned minus sign: +for :math:`u_r = A\sqrt{R}f(\theta)` and :math:`u_\theta = \sqrt{R}g(\theta)`, +:math:`\nabla\cdot\mathbf u = 0` forces :math:`g' = -\tfrac32 A f`, which +integrates to (A8b)'s sign. Flip it and the divergence becomes +:math:`(0.75\cos(\theta/2) + 2.25\cos(3\theta/2))/\sqrt r` — non-zero purely in +the fault modes. + +**A9b also has the wrong function.** @gthyagi's independent reading of both +papers found that printed A9b carries :math:`\cos(3\theta/2)` where +:math:`\sin(3\theta/2)` belongs, as well as the opposite sign; the form used +here reproduces A8b exactly at :math:`r = R_0` and satisfies the equations, +while printed A9b fails both. It is a typographical error in the paper rather +than an alternative convention. + +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 the paper's force balance is +:math:`\partial_j\tau_{ij} + \partial_i p = 0`. Underworld3 — and therefore this +suite — takes pressure positive in compression, so :attr:`fn_pressure` and +:attr:`FaultedMedium.pressure_polar` are the **negative** of the paper's (A9c). +That is a convention difference, not a defect in the paper, and it is recorded +with the family's other errata; it is also measured rather than asserted, in two +independent ways: + +* symbolically, by this suite's momentum residual, which is formed as + :math:`\nabla\cdot\sigma + \mathbf f` with :math:`\sigma = 2\eta\dot\varepsilon + - p\mathbf I`; +* numerically, by @gthyagi against a UW3 Stokes solve on a Gmsh slit disc — + compared against :math:`-p_{\rm BH96}` the pressure error converges (13.1%, + 6.9%, 3.3% at :math:`h` = 0.20, 0.10, 0.05), and compared against the paper's + own sign it is about 199% at every resolution. + +See Also +-------- +underworld3.analytic._base : the contract this solution satisfies. +""" + +import numpy as np +import sympy + +from ._base import AnalyticSolution + + +class FaultedMedium(AnalyticSolution): + r"""Barr & Houseman's fault-tip solution on a disc of radius :math:`R_0`. + + The fault occupies :math:`\theta = 0` from the tip to the perimeter. Slip is + the jump in fault-parallel velocity across it, :math:`2U_0\sqrt{r/R_0}`, so + the relative slip velocity at the perimeter is :math:`2U_0` — the paper's + normalisation. + + Parameters + ---------- + mesh : Mesh + Supplies the coordinate symbols the exact fields are written in. + U0 : float or sympy.Symbol + Velocity scale. The slip at the perimeter is ``2 * U0``. + R0 : float or sympy.Symbol + Radius of the disc, and the length of the fault. + eta : float or sympy.Symbol + Newtonian viscosity of the medium. + tip : sequence of 2 floats + Where the fault tip sits in the mesh's coordinates. The fault runs from + there along :math:`+x`. + + Notes + ----- + **The domain is the solution's, not the mesh's.** The disc of radius + :math:`R_0` about the tip is where this field means something; a mesh only + supplies the coordinate symbols, and :meth:`sample_points` returns points on + the disc rather than on the mesh. Validating the solution therefore does not + need the slit-disc mesh that solving it does — which is just as well, since + Underworld cannot yet build one (see :meth:`apply_boundary_conditions`). + + **The field is multivalued about 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**, and + both representations here put it there deliberately: the NumPy evaluators + take ``arctan2(y, x) mod 2*pi``, and the symbolic Cartesian fields use a + half-angle form of :math:`\theta` whose own cut is the positive :math:`x` + axis. A bare ``atan2`` would put the cut on the *negative* :math:`x` axis and + silently return the wrong face of the fault. + + The parameters may be SymPy symbols. That is not a convenience: a symbolic + check with :math:`U_0 = R_0 = \eta = 1` is satisfied by a transcription + carrying the wrong power of :math:`R_0` in the pressure, so leaving them free + makes the verification strictly stronger for the same runtime. The NumPy + evaluators refuse a symbolic instance, naming the reason. + """ + + dim = 2 + solves = "stokes" + nonlinear = False + + # Half a dozen terms per field: every gate on it runs in well under a + # second, so it belongs in the per-PR tier. + expensive_to_validate = False + + reference = ( + "T. D. Barr & G. A. Houseman, Geophys. J. Int. 125, 473-490 (1996), " + "Appendix (A1)-(A9), plane strain, n = 1. Implementation follows " + "@gthyagi; verified by him against both BH papers and against numerical " + "UW3 models (PR #550)." + ) + + eqn_velocity = ( + r"u_r = \tfrac{U_0}{4}\left[R^2(\sin\theta - \sin 3\theta) " + r"- 2R^3(\sin 2\theta - \sin 4\theta) " + r"+ \sqrt{R}\left(\cos\tfrac{\theta}{2} + 3\cos\tfrac{3\theta}{2}\right)\right]" + ) + eqn_pressure = ( + r"p = -\frac{\eta U_0}{R_0}\left[-2R\sin\theta + 3R^2\sin 2\theta " + r"+ R^{-1/2}\cos\tfrac{\theta}{2}\right] \quad (R = r/R_0)" + ) + eqn_viscosity = r"\eta = \text{const}" + eqn_bodyforce = r"\mathbf f = 0 \quad\text{(driven by the perimeter datum)}" + + def __init__(self, mesh, U0=1.0, R0=1.0, eta=1.0, tip=(0.0, 0.0)): + super().__init__(mesh) + + 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) + self.tip = np.asarray(tip, dtype=float) + + self._symbols = None + self._stress_fn = None + + velocity, pressure, strainrate = self._cartesian() + + self.set_fields( + velocity=velocity, + pressure=pressure, + viscosity=self.eta, + bodyforce=(0, 0), + strainrate=strainrate, + ) + + @property + def _is_symbolic(self): + return any( + isinstance(value, sympy.Basic) for value in (self.U0, self.R0, self.eta) + ) + + # --------------------------------------------------------------- polar form + @property + def symbols(self): + r"""The polar symbols ``(r, theta)`` the fault-frame expressions use. + + ``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 away a substitution 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): + """``(u_r, u_theta, p)`` in the fault frame, in UW3's pressure sign.""" + + 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)) + ) + + # Their (A9c) is extension-positive; UW3 is compression-positive, so the + # sign is flipped HERE, once, and nowhere else. See the module notes. + 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 + + def _strainrate_polar(self): + r"""``(e_rr, e_r_theta, e_theta_theta)`` from the polar velocity.""" + + r, t = self.symbols + u_r, u_t, _p = self._polar() + + e_rr = sympy.diff(u_r, r) + e_tt = sympy.diff(u_t, t) / r + u_r / r + e_rt = (r * sympy.diff(u_t / r, r) + sympy.diff(u_r, t) / r) / 2 + + return e_rr, e_rt, e_tt + + @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): + r"""Pressure in the fault frame, **positive in compression**. + + This is :math:`-1` times the paper's (A9c), which is extension-positive. + The suite follows Underworld3's convention rather than any paper's; see + the module notes for the two measurements that pin it. + """ + + 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)), + ) + + def fault_normal_velocity(self): + r"""The fault-normal velocity :math:`u_\theta` on :math:`\theta = 0`. + + A SymPy expression in ``r``. It is the *same* on both faces — that + continuity is one of the paper's three fault conditions — so it is the + datum a split-node model imposes on each face while leaving the + tangential traction natural. The tangential velocity is what jumps; see + :meth:`slip`. + """ + + _r, t = self.symbols + _u_r, u_t = self.velocity_polar + + return sympy.simplify(u_t.subs(t, 0)) + + def slip(self, r): + r"""Fault slip :math:`2U_0\sqrt{r/R_0}` at radius ``r`` from the tip. + + The jump in fault-parallel velocity between the two faces. The + :math:`\sqrt r` dependence is the :math:`m = -1/2` mode, and it is the + quantity a discrete model can be asked to reproduce — unlike the stress, + which is singular at the tip. + """ + + self._refuse_if_symbolic() + + return 2.0 * self.U0 * np.sqrt(np.asarray(r, dtype=float) / self.R0) + + # ----------------------------------------------------- Cartesian, on mesh.X + def _cartesian(self): + """The contract's fields, in the mesh coordinates. + + The polar expressions with :math:`(r, \\theta)` replaced by their + Cartesian forms about the tip, and the polar components rotated. The + rotation uses :math:`\\cos\\theta = x/r` and :math:`\\sin\\theta = y/r` + directly rather than trigonometry of an arctangent — same value, far + smaller expression. + """ + + r_sym, t_sym = self.symbols + + x = self.mesh.X[0] - self.tip[0] + y = self.mesh.X[1] - self.tip[1] + r = sympy.sqrt(x**2 + y**2) + + # theta in [0, 2*pi) with the branch cut ON the fault, from the + # half-angle identity tan(theta/2) = (r - x)/y. + # + # The two branches are the SAME function — (r - x)(r + x) = y^2, so the + # first is the second with both arguments of atan2 scaled by the + # positive quantity (r + x), which leaves the angle alone. The split is + # purely for conditioning: near the fault r - x is a difference of two + # nearly equal numbers, and forming it directly loses the relative + # accuracy of theta as 1/theta^2. Measured, at 1e-6 radians off the + # fault: momentum residual 3e-6 formed directly, 3e-15 formed this way. + # Both branches are also division-free, so the unselected one cannot + # raise on the axis where its denominator would have vanished. + theta = 2 * sympy.Piecewise( + (sympy.atan2(y**2, y * (r + x)), x > 0), + (sympy.atan2(r - x, y), True), + ) + + into_cartesian = [(r_sym, r), (t_sym, theta)] + + u_r, u_t, pressure = ( + expression.subs(into_cartesian, simultaneous=True) + for expression in self._polar() + ) + e_rr, e_rt, e_tt = ( + expression.subs(into_cartesian, simultaneous=True) + for expression in self._strainrate_polar() + ) + + cos_t, sin_t = x / r, y / r + + velocity = ( + u_r * cos_t - u_t * sin_t, + u_r * sin_t + u_t * cos_t, + ) + strainrate = self._rotate(e_rr, e_rt, e_tt, cos_t, sin_t) + + return velocity, pressure, strainrate + + @staticmethod + def _rotate(a_rr, a_rt, a_tt, cos_t, sin_t): + """A symmetric polar tensor as its ``[[xx, xy], [xy, yy]]`` components.""" + + return ( + ( + a_rr * cos_t**2 - 2 * a_rt * cos_t * sin_t + a_tt * sin_t**2, + (a_rr - a_tt) * cos_t * sin_t + a_rt * (cos_t**2 - sin_t**2), + ), + ( + (a_rr - a_tt) * cos_t * sin_t + a_rt * (cos_t**2 - sin_t**2), + a_rr * sin_t**2 + 2 * a_rt * cos_t * sin_t + a_tt * cos_t**2, + ), + ) + + def sample_points(self, count=12): + r"""Points on the disc, away from the tip and off both faces of the fault. + + The generic sampler is a box, which this solution is not posed on, and it + would land on both of the places the field is not a function of position: + the tip, where the pressure and stress diverge, and the fault itself, + where the velocity is genuinely discontinuous. + + The fixed points carry the sampling that matters: a decade of radii, the + perimeter, and a point :math:`10^{-6}` radians off *each* face — close + enough that a wrongly placed branch cut, or a badly conditioned + :math:`\theta`, shows up in the residual gates rather than only in a + dedicated test. + """ + + self._refuse_if_symbolic() + + rng = np.random.default_rng(20260815) + + radii = self.R0 * rng.uniform(0.15, 0.95, count) + angles = rng.uniform(1.0e-3, 2.0 * np.pi - 1.0e-3, count) + + near = 1.0e-6 + radii = np.r_[radii, self.R0 * np.array([0.05, 0.2, 0.5, 0.5, 1.0, 1.0])] + angles = np.r_[ + angles, + np.array([np.pi / 3, 3 * np.pi / 2, near, 2 * np.pi - near, np.pi / 2, np.pi]), + ] + + return np.column_stack( + [ + self.tip[0] + radii * np.cos(angles), + self.tip[1] + radii * np.sin(angles), + ] + ) + + # --------------------------------------------------- boundaries and solving + @property + def boundaries(self): + """The perimeter of the disc and the two faces of the fault.""" + + return ["Perimeter", "FaultUpper", "FaultLower"] + + def apply_boundary_conditions(self, solver): + """Refused: a fault is not a set of wall conditions. + + This solution is the one member of the family whose conditions the + mixins cannot express, and bending the contract to pretend otherwise + would be worse than saying so. Two reasons, both real: + + * The fault is an **internal** boundary whose two faces must be separate + degrees of freedom at the same coordinates. That is a property of the + mesh, not of the solver, and Underworld cannot yet build it for a fault + that reaches the domain boundary (#549). + * Its conditions are per-component — the fault-normal velocity is + prescribed on both faces while the tangential traction is left natural + — rather than the whole-velocity or whole-normal conditions + :class:`FixedWalls` and :class:`FreeSlipWalls` apply. + + The pieces a model needs are all here: + :meth:`boundary_velocity` (the perimeter datum, their A8), + :meth:`fault_normal_velocity` (the common datum on both faces) and + :meth:`slip` (the answer to check against). @gthyagi's Gmsh slit-disc + model in PR #550 is a worked example of assembling them. + """ + + raise NotImplementedError( + "FaultedMedium is posed on a slit disc: the fault is an internal " + "boundary with two coincident faces, and its conditions are " + "per-component (fault-normal velocity prescribed, tangential " + "traction natural). Neither mixin expresses that, and UW3 cannot yet " + "mesh a fault that reaches the domain boundary (#549). Build the " + "conditions from boundary_velocity(), fault_normal_velocity() and " + "slip() — see the docstring." + ) + + # ------------------------------------------------------- NumPy, fault frame + def evaluate_velocity(self, coords): + r"""Velocity at Cartesian ``coords``, in the mesh's coordinates. + + 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. + """ + + _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)) + + 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)) + ) + + 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. + + Positive in compression, as everywhere else in this suite. 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) + self._refuse_the_tip(r, "pressure") + + 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) + ) + + def evaluate_traction(self, coords, normal): + r"""Traction :math:`\sigma\cdot\hat n` at Cartesian ``coords``. + + With :math:`\sigma = 2\eta\dot\varepsilon - p\mathbf I` and the pressure + positive in compression — the suite's convention, and the same total + stress the paper writes as :math:`\tau + p_{\rm BH}\mathbf I`. 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) + self._refuse_the_tip(r, "stress") + + s_rr, s_rt, s_tt = self._stress_polar_numeric(r, t) + (sxx, sxy), (_syx, syy) = self._rotate(s_rr, s_rt, s_tt, np.cos(t), np.sin(t)) + + n = np.asarray(normal, dtype=float) + if n.ndim == 1: + n = np.broadcast_to(n, (len(r), 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: + r_sym, t_sym = self.symbols + e_rr, e_rt, e_tt = self._strainrate_polar() + pressure = self.pressure_polar + two_eta = 2 * self.eta + + self._stress_fn = sympy.lambdify( + (r_sym, t_sym), + [ + two_eta * e_rr - pressure, + two_eta * e_rt, + two_eta * e_tt - pressure, + ], + "numpy", + ) + + return self._stress_fn(r, t) + + def _polar_of(self, coords): + """``(r, theta, r/R0)`` about the tip, with the cut ON the fault. + + ``arctan2`` is accurate close to either face here — it is only the + *symbolic* Cartesian form that needs care there — so the NumPy path takes + the direct route and shifts the branch cut by a modulus. + """ + + self._refuse_if_symbolic() + + X = np.asarray(coords, dtype=float) + if X.ndim != 2 or X.shape[1] != 2: + raise ValueError("coords must have shape (N, 2)") + + x = X[:, 0] - self.tip[0] + y = X[:, 1] - self.tip[1] + + r = np.hypot(x, y) + t = np.mod(np.arctan2(y, x), 2.0 * np.pi) + + return r, t, r / self.R0 + + def _refuse_if_symbolic(self): + if self._is_symbolic: + raise ValueError( + "this solution was built with symbolic parameters; give U0, R0 " + "and eta numeric values to evaluate it" + ) + + @staticmethod + def _refuse_the_tip(r, quantity): + if np.any(r == 0.0): + raise ValueError( + f"the {quantity} is singular at the fault tip; exclude r = 0 " + f"(the velocity is defined there — use evaluate_velocity)" + ) diff --git a/tests/analytic_full/test_analytic_full_family.py b/tests/analytic_full/test_analytic_full_family.py index e83c2dae..7e6d7d2d 100644 --- a/tests/analytic_full/test_analytic_full_family.py +++ b/tests/analytic_full/test_analytic_full_family.py @@ -97,6 +97,11 @@ {"matrix_viscosity": 3.0}, {"matrix_viscosity": 0.5, "viscosity_ratio": 100.0}, ], + "FaultedMedium": [ + {"U0": 2.5}, + {"R0": 3.0, "eta": 0.25}, + {"U0": 1.3, "R0": 2.0, "eta": 0.7, "tip": (0.3, -0.2)}, + ], } # Defaults first, then every off-default case. @@ -206,7 +211,11 @@ def test_strain_rate_matches_the_velocity(cache, meshes, name, kw): assert _validation.strainrate_consistency(sol, sol.sample_points(count=8)) < 1.0e-8 -FORCED = [n for n in STOKES if n != "EllipticalInclusion"] +# Driven entirely by their boundaries: no body force to flip. Named rather than +# detected, and asserted below. +FORCE_FREE = {"EllipticalInclusion", "FaultedMedium"} + +FORCED = [n for n in STOKES if n not in FORCE_FREE] @pytest.mark.parametrize("name", FORCED) @@ -216,8 +225,8 @@ def test_flipping_the_body_force_breaks_the_momentum_balance(cache, meshes, name The negative control, run here over the whole family. Without it the momentum gate is only an assertion that a small number is small. - EllipticalInclusion is excluded because it is boundary-driven and has no body - force to flip — asserted below rather than skipped silently. + The FORCE_FREE solutions are excluded because they are boundary-driven and + have no body force to flip — asserted below rather than skipped silently. """ sol = _built(cache, meshes, name, {}) @@ -238,10 +247,11 @@ def test_flipping_the_body_force_breaks_the_momentum_balance(cache, meshes, name ) -def test_the_inclusion_really_has_no_body_force(cache, meshes): - """The stated grounds for excluding it from the negative control above.""" +@pytest.mark.parametrize("name", sorted(FORCE_FREE)) +def test_the_force_free_solutions_really_have_no_body_force(cache, meshes, name): + """The stated grounds for excluding them from the negative control above.""" - sol = _built(cache, meshes, "EllipticalInclusion", {}) + sol = _built(cache, meshes, name, {}) assert all(component == 0 for component in sol.fn_bodyforce) diff --git a/tests/test_1028_analytic_parameter_sweep.py b/tests/test_1028_analytic_parameter_sweep.py index b2db0402..a1ea6adb 100644 --- a/tests/test_1028_analytic_parameter_sweep.py +++ b/tests/test_1028_analytic_parameter_sweep.py @@ -46,6 +46,13 @@ {"matrix_viscosity": 3.0}, {"matrix_viscosity": 0.5, "viscosity_ratio": 100.0}, ], + # The tip position is this solution's equivalent of a viscosity off unity: + # at the default it sits on the origin, where a dropped offset cancels. + "FaultedMedium": [ + {"U0": 2.5}, + {"R0": 3.0, "eta": 0.25}, + {"U0": 1.3, "R0": 2.0, "eta": 0.7, "tip": (0.3, -0.2)}, + ], } ALL_STOKES = sorted( @@ -167,7 +174,12 @@ def test_strain_rate_matches_the_velocity_off_default(cache, meshes, name, kw): # -------------------------------------------------------------------------- -FORCED = [n for n in STOKES if n != "EllipticalInclusion"] +# Two solutions are driven entirely by their boundaries and have no body force +# to flip. That is a property of the problems, so they are named here and the +# claim is asserted below rather than skipped silently. +FORCE_FREE = {"EllipticalInclusion", "FaultedMedium"} + +FORCED = [n for n in STOKES if n not in FORCE_FREE] @pytest.mark.parametrize("name", FORCED) @@ -179,9 +191,10 @@ def test_flipping_the_body_force_breaks_the_momentum_balance(cache, meshes, name :math:`\nabla\cdot\sigma + \mathbf f = 0`; if the suite tolerated the other sign it would not be validating that. - EllipticalInclusion is excluded because it has no body force to flip — it is - driven entirely by its boundary. That is a property of the problem, so it is - excluded by name here and asserted below rather than skipped silently. + The solutions in FORCE_FREE are excluded because they have no body force to + flip — they are driven entirely by their boundaries. That is a property of + those problems, so they are excluded by name here and asserted below rather + than skipped silently. """ sol = _built(cache, meshes, name, {}) @@ -202,10 +215,11 @@ def test_flipping_the_body_force_breaks_the_momentum_balance(cache, meshes, name ) -def test_the_inclusion_really_has_no_body_force(cache, meshes): - """The stated grounds for excluding it from the negative control above.""" +@pytest.mark.parametrize("name", sorted(FORCE_FREE)) +def test_the_force_free_solutions_really_have_no_body_force(cache, meshes, name): + """The stated grounds for excluding them from the negative control above.""" - sol = _built(cache, meshes, "EllipticalInclusion", {}) + sol = _built(cache, meshes, name, {}) assert all(component == 0 for component in sol.fn_bodyforce) diff --git a/tests/test_1029_analytic_faulted_medium.py b/tests/test_1029_analytic_faulted_medium.py new file mode 100644 index 00000000..5291bcb7 --- /dev/null +++ b/tests/test_1029_analytic_faulted_medium.py @@ -0,0 +1,502 @@ +"""Barr & Houseman (1996) faulted medium — 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 + div sigma = eta lap(u) - grad(p) = 0 (their eq 3, in UW3's sign) + 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. + +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. + +Two things are new here relative to the fault-frame verification, and both come +from joining the family: the pressure is carried in UW3's sign rather than the +paper's (asserted below against the printed A9c), and the same field is also +exposed in the mesh coordinates, where the family's residual gates reach it. + +Run: pixi run python -m pytest tests/test_1029_analytic_faulted_medium.py -v +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.analytic import _validation + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +@pytest.fixture(scope="module") +def mesh(): + """A carrier for the coordinate symbols. + + The solution is posed on a disc about the fault tip, which is not this box — + and does not need to be. `sample_points` returns points on the disc, and the + gates evaluate symbolic expressions rather than anything the mesh holds. + """ + + return uw.meshing.StructuredQuadBox( + elementRes=(4, 4), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3 + ) + + +def _symbolic(mesh): + """The solution with U0, R0 and eta left free.""" + + U0, R0, eta = sympy.symbols("U_0 R_0 eta", positive=True) + + return uw.analytic.FaultedMedium(mesh, U0=U0, R0=R0, eta=eta) + + +def _polar_operators(u_r, u_t, p, r, t, eta): + """(div u, momentum_x, momentum_y) for a field given in polar form. + + The momentum operator is `eta lap(u) - grad(p)`: the divergence of + `sigma = 2 eta edot - p I` at constant viscosity, i.e. UW3's compression- + positive pressure. The paper writes `+ grad(p)` because its pressure is + extension-positive; the two statements are the same one. + """ + + 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(mesh): + """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. + """ + + sol = _symbolic(mesh) + 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(mesh): + """All THREE of the paper's fault conditions, with symbolic parameters. + + Zero shear traction, continuous normal velocity, continuous normal stress. + """ + + sol = _symbolic(mesh) + 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 + + # Normal STRESS continuity — the third of the paper's three fault + # conditions. sigma = 2 eta edot - p I with p positive in compression, so + # 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_pressure_is_the_negative_of_the_published_one(mesh): + """The one convention this solution does not share with its source. + + Barr & Houseman take extension as positive, so their force balance is + `d_j tau_ij + d_i p = 0`. UW3 — and so this suite — takes pressure positive + in compression. The flip is applied once, and asserted here against the + paper's printed (A9c) rather than left to a comment, because a sign that is + only documented is a sign that can drift. + + Independent measurement: @gthyagi's UW3 Stokes solve on a Gmsh slit disc + converges against -p_BH96 (13.1%, 6.9%, 3.3% at h = 0.20, 0.10, 0.05) and + sits at about 199% at every resolution against the paper's own sign. + """ + + sol = _symbolic(mesh) + r, t = sol.symbols + R = r / sol.R0 + + published_A9c = (sol.eta * sol.U0 / sol.R0) * ( + -2 * R * sympy.sin(t) + + 3 * R**2 * sympy.sin(2 * t) + + sympy.cos(t / 2) / sympy.sqrt(R) + ) + + assert sympy.simplify(sol.pressure_polar + published_A9c) == 0 + + +def test_the_slip_is_the_published_normalisation(mesh): + """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 = uw.analytic.FaultedMedium(mesh, 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(mesh): + """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 = uw.analytic.FaultedMedium(mesh, 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 = sol.evaluate_velocity(pts) + pressure = sol.evaluate_pressure(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(mesh): + """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 = uw.analytic.FaultedMedium(mesh) + eps = 1e-7 + + above = sol.evaluate_velocity(np.array([[0.5, +eps]])) + below = sol.evaluate_velocity(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_velocity(np.array([[-0.5, +eps]])) + left_dn = sol.evaluate_velocity(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(mesh): + """The Dirichlet datum a solver would impose is the solution at r = R0.""" + + sol = uw.analytic.FaultedMedium(mesh, 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_the_fault_normal_datum_is_common_to_both_faces(mesh): + """The datum a split-node model imposes on each face. + + It is `u_theta` on theta = 0, and the SAME expression on theta = 2 pi — that + equality is fault condition 2, and it is why one datum serves both faces. + """ + + sol = uw.analytic.FaultedMedium(mesh, U0=1.0, R0=1.5, eta=1.0) + r, t = sol.symbols + _u_r, u_t = sol.velocity_polar + + datum = sol.fault_normal_velocity() + + assert sympy.simplify(datum - u_t.subs(t, 0)) == 0 + assert sympy.simplify(datum - u_t.subs(t, 2 * sympy.pi)) == 0 + assert not datum.free_symbols - {r}, "the datum depends on r only" + + +def test_a_degenerate_geometry_is_refused(mesh): + with pytest.raises(ValueError, match="positive"): + uw.analytic.FaultedMedium(mesh, R0=0.0) + with pytest.raises(ValueError, match="positive"): + uw.analytic.FaultedMedium(mesh, 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. + sol = uw.analytic.FaultedMedium(mesh) + with pytest.raises(ValueError, match="singular at the fault tip"): + sol.evaluate_pressure(np.array([[0.0, 0.0]])) + assert np.allclose(sol.evaluate_velocity(np.array([[0.0, 0.0]])), 0.0) + + symbolic = uw.analytic.FaultedMedium(mesh, U0=sympy.Symbol("U_0", positive=True)) + with pytest.raises(ValueError, match="symbolic parameters"): + symbolic.evaluate_velocity(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(mesh): + """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 = uw.analytic.FaultedMedium(mesh) + 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]) + + +# --------------------------------------------------------------------------- +# Joining the family: the same field in the mesh coordinates +# --------------------------------------------------------------------------- + + +def test_the_cartesian_fields_are_the_fault_frame_ones(mesh): + """`fn_velocity` and `fn_pressure` must be the polar solution, rotated. + + The contract's fields are a second representation of the same solution, so + they need pinning to the first — this is the join the integration adds, and + the one place a rotation or a branch-cut slip could enter unseen. + """ + + sol = uw.analytic.FaultedMedium(mesh, 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 + + points = sol.sample_points(count=8) + radii = np.hypot(points[:, 0], points[:, 1]) + angles = np.mod(np.arctan2(points[:, 1], points[:, 0]), 2.0 * np.pi) + + for name, expression, reference in ( + ("velocity_x", sol.fn_velocity[0, 0], None), + ("velocity_y", sol.fn_velocity[0, 1], None), + ("pressure", sol.fn_pressure, p_sym), + ): + mine = _validation.sample(sol, expression, points) + + for k, (rr, tt) in enumerate(zip(radii, angles)): + subs = {r_sym: float(rr), t_sym: float(tt)} + if reference is None: + ur = float(u_r_sym.subs(subs)) + ut = float(u_t_sym.subs(subs)) + expect = ( + ur * np.cos(tt) - ut * np.sin(tt) + if name == "velocity_x" + else ur * np.sin(tt) + ut * np.cos(tt) + ) + else: + expect = float(reference.subs(subs)) + + assert mine[k] == pytest.approx(expect, rel=1e-9, abs=1e-11), name + + +def test_the_residual_gates_hold_with_the_tip_off_the_origin(mesh): + """The oracle-free gates, at a tip the conformance sweep never uses. + + The sweep builds every solution from its defaults, and the default tip is + the origin — where a dropped offset is invisible. Moving it is this + solution's equivalent of moving a viscosity off unity. + """ + + sol = uw.analytic.FaultedMedium(mesh, U0=1.3, R0=2.0, eta=0.7, tip=(0.3, -0.2)) + points = sol.sample_points(count=8) + + assert _validation.incompressibility_residual(sol, points) < 1.0e-8 + assert _validation.momentum_residual(sol, points) < 1.0e-8 + assert _validation.strainrate_consistency(sol, points) < 1.0e-8 + + # The offset is real: the fields are not the same function of position. + at_origin = uw.analytic.FaultedMedium(mesh, U0=1.3, R0=2.0, eta=0.7) + moved = _validation.sample(sol, sol.fn_velocity[0, 0], points) + unmoved = _validation.sample(at_origin, at_origin.fn_velocity[0, 0], points) + assert np.abs(moved - unmoved).max() > 0.1 * np.abs(moved).max() + + +def test_the_sample_points_avoid_the_tip_and_land_on_both_faces(mesh): + """The override is what keeps the gates off the two places this field is not + a function of position, and on the one place it is hardest. + """ + + sol = uw.analytic.FaultedMedium(mesh, R0=2.0, tip=(0.3, -0.2)) + points = sol.sample_points(count=8) + + offset = points - sol.tip + radii = np.hypot(offset[:, 0], offset[:, 1]) + angles = np.mod(np.arctan2(offset[:, 1], offset[:, 0]), 2.0 * np.pi) + + assert radii.min() > 0.0 + assert radii.max() <= sol.R0 * (1.0 + 1.0e-12) + + # One point just above the fault and one just below it. + assert angles.min() < 1.0e-5 + assert angles.max() > 2.0 * np.pi - 1.0e-5 + + +def test_it_is_registered_in_the_family(mesh): + """Listed, described, and declaring the same things the others declare.""" + + assert "FaultedMedium" in uw.analytic.available() + assert uw.analytic.is_available("FaultedMedium") + assert "Barr" in uw.analytic.describe("FaultedMedium") + + assert uw.analytic.FaultedMedium.solves == "stokes" + assert uw.analytic.FaultedMedium.symbolic is True + assert uw.analytic.FaultedMedium.expensive_to_validate is False + assert uw.analytic.FaultedMedium.dim == 2 + assert uw.analytic.FaultedMedium.reference + + +def test_it_refuses_to_pretend_the_fault_is_a_wall(mesh): + """The one member of the family that cannot use the mixins says so. + + Applying the perimeter datum and quietly leaving the fault unconstrained + would solve a different problem and report a plausible error, which is worse + than refusing. + """ + + sol = uw.analytic.FaultedMedium(mesh) + + assert sol.boundaries == ["Perimeter", "FaultUpper", "FaultLower"] + + velocity = uw.discretisation.MeshVariable("Ubh", mesh, 2, degree=2) + pressure = uw.discretisation.MeshVariable("Pbh", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=velocity, pressureField=pressure) + + with pytest.raises(NotImplementedError, match="slit disc"): + sol.apply_boundary_conditions(stokes) + + +def test_it_has_no_body_force(mesh): + """Boundary-driven, like the elliptical inclusion — and for the same reason + it is excluded from the body-force negative control in the sweeps. + """ + + sol = uw.analytic.FaultedMedium(mesh) + + assert all(component == 0 for component in sol.fn_bodyforce) + + +def test_flipping_the_pressure_sign_breaks_the_momentum_balance(mesh): + """The negative control this solution gets INSTEAD of the body-force one. + + The family pins its body-force convention by flipping the force and watching + the momentum residual go to order unity. There is no force here to flip, so + the gate would otherwise be an assertion that a small number is small. + + The pressure sign is what it certifies instead, and it certifies it just as + sharply: with the paper's extension-positive pressure the residual is 1.06, + with UW3's compression-positive one it is 3.6e-16. Note which gate fires — + `tr(sigma) + d p` is 5.7e-16 either way, because `set_fields` builds the + stress from the pressure and tracelessness cannot see the sign of a term it + cancels by construction. + """ + + sol = uw.analytic.FaultedMedium(mesh) + points = sol.sample_points(count=8) + + assert _validation.momentum_residual(sol, points) < 1.0e-8 + + flipped = uw.analytic.FaultedMedium(mesh) + flipped.set_fields( + velocity=list(sol.fn_velocity), + pressure=-sol.fn_pressure, + viscosity=sol.fn_viscosity, + bodyforce=(0, 0), + strainrate=sol.fn_strainrate.tolist(), + ) + + assert _validation.momentum_residual(flipped, points) > 1.0e-2, ( + "the momentum gate cannot see the pressure sign it is supposed to " + "certify for this solution" + ) From 421b4ab30c1327578be36ffc1f6469fcce8f2c5c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 23:26:49 +1000 Subject: [PATCH 07/11] docs: a solution with a fault in it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The analytic subsystem document gains the faulted medium: why it is worth having (the only absolute standard in the suite for a fault calculation), the three ways it does not fit the contract and what was done about each, the conditioning measurement behind the Cartesian theta, and the erratum. Also corrects the per-PR tier's file glob in the two-tiers table — it read test_102[0-8], and the new file is test_1029. The CI batch glob is tests/test_102*py, so the file was always covered; the table was not. Underworld development team with AI support from Claude Code --- .../subsystems/analytic-solutions.md | 88 ++++++++++++++++++- 1 file changed, 87 insertions(+), 1 deletion(-) diff --git a/docs/developer/subsystems/analytic-solutions.md b/docs/developer/subsystems/analytic-solutions.md index 2fe4ae6b..6f1f78ff 100644 --- a/docs/developer/subsystems/analytic-solutions.md +++ b/docs/developer/subsystems/analytic-solutions.md @@ -390,7 +390,7 @@ before this suite existed, so they cannot ride on every PR. | tier | what it covers | cost | where | |---|---|---|---| -| per-PR | every gate, on every solution that is cheap to validate; one canonical parameter case for SolKx/SolKz | **4m07s** | `tests/test_101[5-9]_analytic_*`, `tests/test_102[0-8]_analytic_*` — matched by the CI batch globs | +| per-PR | every gate, on every solution that is cheap to validate; one canonical parameter case for SolKx/SolKz | **4m07s** | `tests/test_101[5-9]_analytic_*`, `tests/test_102[0-9]_analytic_*` — matched by the CI batch globs | | full family | every gate, on every solution, over the whole parameter table | **8m34s** | `tests/analytic_full/` — matched by nothing in CI | ```bash @@ -509,6 +509,91 @@ the cost was imagined. | C ternary and `&&`/`\|\|` | SolH guards its zero modes with `(n!=0 \|\| m!=0) ? … : …` | | `resolve_branches` | guards on loop indices have an answer at transcription time. Left unresolved, `evaluate_block` reads every branch in order and each guarded variable keeps the *last* one — in SolH that silently zeroes two velocity components, which looks plausible rather than broken | +## A solution with a fault in it + +`FaultedMedium` (Barr & Houseman 1996, `analytic/barr_houseman.py`) is the linear +plane-strain case of their Appendix: a fault that terminates *inside* a viscous +medium, on a disc of radius $R_0$ with the tip at the centre and the fault +running out to the perimeter. It carries zero shear traction on the fault, +continuous normal velocity and continuous normal stress across it, and slip +$2U_0\sqrt{r/R_0}$. + +It is here because it is the only absolute standard in the suite for a fault +calculation. Everything else a fault model can be measured against is another +discretisation. + +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 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. Drop the half-integer term +and the slip vanishes; there is a test that does exactly that. + +### Three ways it does not fit the contract, and what was done about each + +**It is not a function of position.** The field is multivalued about the tip — +that is what a fault is — so it is a function of $(r,\theta)$ with +$\theta\in[0,2\pi)$, and the branch cut has to lie *on the fault*, which is where +the field is genuinely discontinuous. A bare `atan2` puts the cut on the negative +$x$ axis and silently returns the wrong face. + +The class therefore carries **two representations of one solution**: the polar +fault-frame expressions, in which the fault conditions can be stated exactly at +$\theta = 0$ and $2\pi$, and the contract's `fn_*` fields on `mesh.X`, which is +what the residual gates and a solver consume. A test pins the second to the +first. + +Getting $\theta$ into the Cartesian form is where the care went. From +$\tan(\theta/2) = (r-x)/y$, the cut of $2\,\mathrm{atan2}(r-x,\,y)$ falls on the +positive $x$ axis — the fault — which is right. But near the fault $r - x$ is a +difference of two nearly equal numbers, and the relative accuracy of $\theta$ +degrades as $1/\theta^2$. Since $(r-x)(r+x) = y^2$, scaling both arguments by the +positive quantity $(r+x)$ leaves the angle alone and removes the cancellation: +$2\,\mathrm{atan2}(y^2,\,y(r+x))$, used where $x > 0$. Measured, at $10^{-6}$ +radians off the fault, the momentum residual is **3e-6** formed directly and +**3e-15** formed this way — so `sample_points` can put a point just off each face +and the gates still mean something there. + +**It has no walls.** `apply_boundary_conditions` **refuses**, and that refusal is +deliberate. The fault is an internal boundary whose two faces must be separate +degrees of freedom at coincident coordinates — a property of the mesh, not of the +solver, and one Underworld cannot yet build for a fault that reaches the domain +boundary (#549). Its conditions are also per-component (fault-normal velocity +prescribed on both faces, tangential traction natural), which is not what +`FixedWalls` or `FreeSlipWalls` apply. Applying the perimeter datum and quietly +leaving the fault unconstrained would solve a different problem and report a +plausible error, which is worse than refusing. The pieces a model needs — +`boundary_velocity()`, `fault_normal_velocity()`, `slip()` — are exposed instead. + +**It has no body force.** Like `EllipticalInclusion` it is driven entirely by its +boundary, so the family's body-force negative control cannot fire and it is +excluded from that control by name in both sweeps. It is not left without one: +for this solution the momentum residual certifies the **pressure sign** instead, +and just as sharply — 3.6e-16 with UW3's compression-positive pressure, **1.06** +with the paper's extension-positive one, while $\mathrm{tr}\,\sigma + d\,p$ sits +at 5.7e-16 either way because `set_fields` cancels that term by construction. + +### The erratum + +Barr & Houseman's pressure is extension-positive; ours is not, so `fn_pressure` +is $-1$ times their (A9c). Their printed (A9b) also disagrees with their own +(A8b) in the sign of the half-integer sine terms, and carries $\cos(3\theta/2)$ +where $\sin(3\theta/2)$ belongs. All three are recorded, with the measurements +that settle them, in the conventions-and-errata note (§3.8) rather than applied +silently. + +### Credit + +The implementation follows @gthyagi's, and he verified it independently three +ways in PR #550: against both BH papers (which is how the A9b function error was +found), against a UW3 Stokes model on a Gmsh slit disc (velocity error 1.92% → +0.74% → 0.40% under refinement, normal-velocity jump at machine zero, and the +pressure sign confirmed numerically at ~199% error with the paper's own sign), +and against his BH92 rectangular-fault-zone model, whose near-tip exponent fits +**-0.498** against the exact **-0.500** on the resolved singular interval. + ## The scalar transport family The Stokes solutions solve for a velocity-and-pressure pair. The scalar solutions @@ -692,6 +777,7 @@ Each vendored reference kernel keeps its original copyright header. | Underworld2 `Velic_sol*` kernels | LGPL-3 | Yes — same licence | | PETSc `ex69.c`, `ex13.c`, `ex24.c`, `ex45.c` | BSD-2-Clause | Yes — permissive; retain the notice | | Schmid & Podladchikov MATLAB (`dwschmid/muskhelishvili`) | BSD-3-Clause | Yes — permissive; retain the notice | +| Barr & Houseman 1996 (GJI 125, 473-490), Appendix | Nothing vendored — transcribed from the published equations | n/a | | `assess` (Kramer et al. 2021) | External, optional dependency | Not vendored; wrapped lazily | ## Optional dependencies From 3c53e040b0ac5476bfec18906d7eb0e135e28c7e Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 23:31:02 +1000 Subject: [PATCH 08/11] Two gates the faulted medium needed and did not have MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The pressure-sign negative control. The family pins its body-force convention by flipping the force and watching the momentum residual go to order unity; this solution has no force to flip, so without a substitute its momentum gate was only an assertion that a small number is small. Flipping the PRESSURE instead moves it from 3.6e-16 to 1.06, so the gate is demonstrably able to fail. Note that tr(sigma) + d p does not move: set_fields builds the stress from the pressure, so tracelessness cannot see the sign of a term it cancels by construction. The compiled path. The family's argument for one SymPy form is that it compiles — the same expression is a Dirichlet datum, a JIT kernel and an error norm. Here that rests on a Piecewise around two atan2 calls, which is exactly what a C printer could get wrong, so evaluate() is checked against the independent NumPy evaluator rather than assumed: they agree to 1.7e-15. This also shows how to place the disc inside a mesh, which is what the eventual benchmark will need. Underworld development team with AI support from Claude Code --- src/underworld3/analytic/barr_houseman.py | 5 +++-- tests/test_1029_analytic_faulted_medium.py | 23 ++++++++++++++++++++++ 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/src/underworld3/analytic/barr_houseman.py b/src/underworld3/analytic/barr_houseman.py index 1ea650db..122e87a5 100644 --- a/src/underworld3/analytic/barr_houseman.py +++ b/src/underworld3/analytic/barr_houseman.py @@ -604,6 +604,7 @@ def _refuse_if_symbolic(self): def _refuse_the_tip(r, quantity): if np.any(r == 0.0): raise ValueError( - f"the {quantity} is singular at the fault tip; exclude r = 0 " - f"(the velocity is defined there — use evaluate_velocity)" + f"the {quantity} is singular at the fault tip; exclude r = 0. " + f"The velocity is defined there — evaluate_velocity returns the " + f"limit, which is zero." ) diff --git a/tests/test_1029_analytic_faulted_medium.py b/tests/test_1029_analytic_faulted_medium.py index 5291bcb7..cbdc2d27 100644 --- a/tests/test_1029_analytic_faulted_medium.py +++ b/tests/test_1029_analytic_faulted_medium.py @@ -381,6 +381,29 @@ def test_the_cartesian_fields_are_the_fault_frame_ones(mesh): assert mine[k] == pytest.approx(expect, rel=1e-9, abs=1e-11), name +def test_the_contract_evaluation_path_compiles(mesh): + """`evaluate` — and so `error` — must work, which means the branch cut has + to survive code generation. + + The whole family's argument for one SymPy form is that it compiles: the same + expression is a Dirichlet datum, a JIT kernel and an error norm. For this + solution that claim rests on a `Piecewise` around two `atan2` calls, which is + the part a C printer could plausibly get wrong, so it is exercised rather + than assumed. Agreement with the independent NumPy path is the check. + + The disc is placed INSIDE the mesh here, unlike everywhere else in this file: + `uw.function.evaluate` needs points the mesh actually contains, while the + residual gates only need the coordinate symbols. + """ + + sol = uw.analytic.FaultedMedium(mesh, R0=0.4, tip=(0.5, 0.5)) + points = sol.sample_points(count=4) + + through_the_jit = np.asarray(sol.evaluate("velocity", points)).reshape(-1, 2) + + assert np.abs(through_the_jit - sol.evaluate_velocity(points)).max() < 1.0e-12 + + def test_the_residual_gates_hold_with_the_tip_off_the_origin(mesh): """The oracle-free gates, at a tip the conformance sweep never uses. From e4f1d9bb891e44d6f67429b44d5856af99562459 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 15 Aug 2026 20:56:54 +0530 Subject: [PATCH 09/11] Rename Barr-Houseman analytic module for 1996 reference Rename analytic/barr_houseman.py to barr_houseman_96.py so the module path identifies the paper that supplies the closed-form Cartesian solution and cannot be confused with the related BH92 study. Update the package import and developer documentation, and add a focused assertion that uw.analytic.barr_houseman_96.FaultedMedium resolves to the registered public class. --- docs/developer/subsystems/analytic-solutions.md | 6 +++++- src/underworld3/analytic/__init__.py | 2 +- .../analytic/{barr_houseman.py => barr_houseman_96.py} | 0 tests/test_1029_analytic_faulted_medium.py | 1 + 4 files changed, 7 insertions(+), 2 deletions(-) rename src/underworld3/analytic/{barr_houseman.py => barr_houseman_96.py} (100%) diff --git a/docs/developer/subsystems/analytic-solutions.md b/docs/developer/subsystems/analytic-solutions.md index 6f1f78ff..5b02d2c7 100644 --- a/docs/developer/subsystems/analytic-solutions.md +++ b/docs/developer/subsystems/analytic-solutions.md @@ -511,13 +511,17 @@ the cost was imagined. ## A solution with a fault in it -`FaultedMedium` (Barr & Houseman 1996, `analytic/barr_houseman.py`) is the linear +`FaultedMedium` (Barr & Houseman 1996, `analytic/barr_houseman_96.py`) is the linear plane-strain case of their Appendix: a fault that terminates *inside* a viscous medium, on a disc of radius $R_0$ with the tip at the centre and the fault running out to the perimeter. It carries zero shear traction on the fault, continuous normal velocity and continuous normal stress across it, and slip $2U_0\sqrt{r/R_0}$. +The year remains in the module name to distinguish this closed-form 1996 +solution from the related 1992 study, which does not provide the corresponding +full Cartesian analytical field. + It is here because it is the only absolute standard in the suite for a fault calculation. Everything else a fault model can be measured against is another discretisation. diff --git a/src/underworld3/analytic/__init__.py b/src/underworld3/analytic/__init__.py index 0b9c7655..df4d4ff6 100644 --- a/src/underworld3/analytic/__init__.py +++ b/src/underworld3/analytic/__init__.py @@ -28,7 +28,7 @@ from ._base import AnalyticSolution, FreeSlipWalls, FixedWalls -from .barr_houseman import FaultedMedium +from .barr_houseman_96 import FaultedMedium from .inclusion import EllipticalInclusion from .kramer import CylindricalStokes from .richards import GardnerSteady, GardnerTransient diff --git a/src/underworld3/analytic/barr_houseman.py b/src/underworld3/analytic/barr_houseman_96.py similarity index 100% rename from src/underworld3/analytic/barr_houseman.py rename to src/underworld3/analytic/barr_houseman_96.py diff --git a/tests/test_1029_analytic_faulted_medium.py b/tests/test_1029_analytic_faulted_medium.py index cbdc2d27..a40e1a68 100644 --- a/tests/test_1029_analytic_faulted_medium.py +++ b/tests/test_1029_analytic_faulted_medium.py @@ -453,6 +453,7 @@ def test_it_is_registered_in_the_family(mesh): assert uw.analytic.is_available("FaultedMedium") assert "Barr" in uw.analytic.describe("FaultedMedium") + assert uw.analytic.barr_houseman_96.FaultedMedium is uw.analytic.FaultedMedium assert uw.analytic.FaultedMedium.solves == "stokes" assert uw.analytic.FaultedMedium.symbolic is True assert uw.analytic.FaultedMedium.expensive_to_validate is False From 7c8ab2168db1a5f5ad49a602137dabb100221e3c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sun, 16 Aug 2026 10:18:55 +1000 Subject: [PATCH 10/11] Drop the mixin references from the faulted medium's refusal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit #578 replaced the boundary-condition mixins with composed functions, so the refusal's explanation now names what the solution actually needs — a component-wise Dirichlet condition on an internal boundary — rather than the two classes that no longer exist. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/analytic-solutions.md | 4 ++-- src/underworld3/analytic/barr_houseman.py | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/developer/subsystems/analytic-solutions.md b/docs/developer/subsystems/analytic-solutions.md index 906dc4b7..8ff23aa1 100644 --- a/docs/developer/subsystems/analytic-solutions.md +++ b/docs/developer/subsystems/analytic-solutions.md @@ -617,8 +617,8 @@ deliberate. The fault is an internal boundary whose two faces must be separate degrees of freedom at coincident coordinates — a property of the mesh, not of the solver, and one Underworld cannot yet build for a fault that reaches the domain boundary (#549). Its conditions are also per-component (fault-normal velocity -prescribed on both faces, tangential traction natural), which is not what -`FixedWalls` or `FreeSlipWalls` apply. Applying the perimeter datum and quietly +prescribed on both faces, tangential traction natural), so none of the +whole-velocity helpers fits. Applying the perimeter datum and quietly leaving the fault unconstrained would solve a different problem and report a plausible error, which is worse than refusing. The pieces a model needs — `boundary_velocity()`, `fault_normal_velocity()`, `slip()` — are exposed instead. diff --git a/src/underworld3/analytic/barr_houseman.py b/src/underworld3/analytic/barr_houseman.py index 122e87a5..23e618de 100644 --- a/src/underworld3/analytic/barr_houseman.py +++ b/src/underworld3/analytic/barr_houseman.py @@ -443,11 +443,11 @@ def boundaries(self): return ["Perimeter", "FaultUpper", "FaultLower"] def apply_boundary_conditions(self, solver): - """Refused: a fault is not a set of wall conditions. + """Refused: this solution needs a mesh Underworld cannot yet build. - This solution is the one member of the family whose conditions the - mixins cannot express, and bending the contract to pretend otherwise - would be worse than saying so. Two reasons, both real: + Every other solution in the family states its conditions here by + composing :func:`~underworld3.analytic.free_slip` and friends. This one + refuses, and the reason is the mesh rather than the boundary conditions: * The fault is an **internal** boundary whose two faces must be separate degrees of freedom at the same coordinates. That is a property of the @@ -455,8 +455,8 @@ def apply_boundary_conditions(self, solver): that reaches the domain boundary (#549). * Its conditions are per-component — the fault-normal velocity is prescribed on both faces while the tangential traction is left natural - — rather than the whole-velocity or whole-normal conditions - :class:`FixedWalls` and :class:`FreeSlipWalls` apply. + — so they need a component-wise Dirichlet condition on an internal + boundary rather than any of the whole-velocity helpers. The pieces a model needs are all here: :meth:`boundary_velocity` (the perimeter datum, their A8), @@ -469,8 +469,8 @@ def apply_boundary_conditions(self, solver): "FaultedMedium is posed on a slit disc: the fault is an internal " "boundary with two coincident faces, and its conditions are " "per-component (fault-normal velocity prescribed, tangential " - "traction natural). Neither mixin expresses that, and UW3 cannot yet " - "mesh a fault that reaches the domain boundary (#549). Build the " + "traction natural). UW3 cannot yet mesh a fault that reaches the " + "domain boundary (#549). Build the " "conditions from boundary_velocity(), fault_normal_velocity() and " "slip() — see the docstring." ) From db55e8e29519d57aa5d2de174c00937a8c897a67 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Mon, 17 Aug 2026 19:35:02 +1000 Subject: [PATCH 11/11] Resolve the merge: the faulted medium lives in uw.analytic MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PR #550 was written against the pre-#571 layout, where the analytic solutions were a compiled extension at function/analytic.pyx and a new elementary solution was re-exported from there. #571 moved the suite to underworld3.analytic and left a deprecation shim, so the PR's edits landed on a file that no longer holds that role — git followed the rename and carried the re-export onto analytic/_reference/_velic.pyx, which is the kernels module, not the namespace. The ported solution is analytic/barr_houseman_96.py (class FaultedMedium), already exported from analytic/__init__.py and registered in the family, so the re-export has nothing left to do and is removed. function/_barr_houseman.py and tests/test_0210_barr_houseman_analytic.py go with it: test_1029_analytic_faulted_medium.py carries all eight of 0210's tests under the same names and ten more, including the contract, registry and tip-off-origin gates the AnalyticSolution base requires. Underworld development team with AI support from Claude Code --- .../analytic/_reference/_velic.pyx | 5 - src/underworld3/function/_barr_houseman.py | 314 ------------------ tests/test_0210_barr_houseman_analytic.py | 238 ------------- 3 files changed, 557 deletions(-) delete mode 100644 src/underworld3/function/_barr_houseman.py delete mode 100644 tests/test_0210_barr_houseman_analytic.py diff --git a/src/underworld3/analytic/_reference/_velic.pyx b/src/underworld3/analytic/_reference/_velic.pyx index 876613f4..93ddb1c3 100644 --- a/src/underworld3/analytic/_reference/_velic.pyx +++ b/src/underworld3/analytic/_reference/_velic.pyx @@ -2,11 +2,6 @@ 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/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py deleted file mode 100644 index 3482c551..00000000 --- a/src/underworld3/function/_barr_houseman.py +++ /dev/null @@ -1,314 +0,0 @@ -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): - # 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) - self._symbols = None - self._stress_fn = None - - @property - def _is_symbolic(self): - return any(isinstance(v, sympy.Basic) - for v in (self.U0, self.R0, self.eta)) - - # ------------------------------------------------------------------ sympy - @property - def symbols(self): - """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 - 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. - """ - velocity = self.evaluate_velocity(coords) - return velocity, self.evaluate_pressure(coords) - - 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)) - + 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)) - ) - 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) - ) - - 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 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. - - 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. - """ - 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 deleted file mode 100644 index 46be5082..00000000 --- a/tests/test_0210_barr_houseman_analytic.py +++ /dev/null @@ -1,238 +0,0 @@ -"""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 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. - -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. - - 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 - - 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(): - """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 - - 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 - - # 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. - - 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) - # 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_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"): - 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])