diff --git a/docs/developer/subsystems/analytic-solutions.md b/docs/developer/subsystems/analytic-solutions.md index 9f7f079b..84143810 100644 --- a/docs/developer/subsystems/analytic-solutions.md +++ b/docs/developer/subsystems/analytic-solutions.md @@ -446,7 +446,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 @@ -565,6 +565,95 @@ 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_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. + +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), 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. + +**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 @@ -749,6 +838,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 diff --git a/src/underworld3/analytic/__init__.py b/src/underworld3/analytic/__init__.py index 6d7690c9..c44c7162 100644 --- a/src/underworld3/analytic/__init__.py +++ b/src/underworld3/analytic/__init__.py @@ -33,6 +33,7 @@ prescribed_velocity, ) +from .barr_houseman_96 import FaultedMedium from .inclusion import EllipticalInclusion from .kramer import CylindricalStokes from .richards import GardnerSteady, GardnerTransient @@ -47,6 +48,7 @@ "AdvectedFront", "EllipticalInclusion", "ErfcDiffusion", + "FaultedMedium", "CylindricalStokes", "GardnerSteady", "GardnerTransient", @@ -78,6 +80,7 @@ "CylindricalStokes": CylindricalStokes, "EllipticalInclusion": EllipticalInclusion, "ErfcDiffusion": ErfcDiffusion, + "FaultedMedium": FaultedMedium, "GardnerSteady": GardnerSteady, "GardnerTransient": GardnerTransient, "Poisson1D": Poisson1D, diff --git a/src/underworld3/analytic/barr_houseman_96.py b/src/underworld3/analytic/barr_houseman_96.py new file mode 100644 index 00000000..23e618de --- /dev/null +++ b/src/underworld3/analytic/barr_houseman_96.py @@ -0,0 +1,610 @@ +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: this solution needs a mesh Underworld cannot yet build. + + 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 + 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 + — 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), + :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). 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 — evaluate_velocity returns the " + f"limit, which is zero." + ) 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..a40e1a68 --- /dev/null +++ b/tests/test_1029_analytic_faulted_medium.py @@ -0,0 +1,526 @@ +"""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_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. + + 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.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 + 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" + )