Barr & Houseman faulted-medium analytic solution in uw.analytic

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

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,