diff --git a/src/underworld3/function/_barr_houseman.py b/src/underworld3/function/_barr_houseman.py new file mode 100644 index 00000000..3482c551 --- /dev/null +++ b/src/underworld3/function/_barr_houseman.py @@ -0,0 +1,314 @@ +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/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..46be5082 --- /dev/null +++ b/tests/test_0210_barr_houseman_analytic.py @@ -0,0 +1,238 @@ +"""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])