diff --git a/docs/advanced/fault-mechanics-examples.md b/docs/advanced/fault-mechanics-examples.md new file mode 100644 index 000000000..f8bf4d053 --- /dev/null +++ b/docs/advanced/fault-mechanics-examples.md @@ -0,0 +1,354 @@ +--- +title: "Fault Mechanics Teaching Examples" +--- + +# Fault Mechanics: Teaching Examples + +Three worked examples built on [split-node faults](split-node-faults.md), +each a single short script (in `figures/fault-examples/`, sharing one +harness `common.py`). All run in minutes on a laptop: the meshes are +modest, because a zero-thickness fault needs no thin feature resolved. + +## The fault-strength ladder + +One fault, one far-field shear drive ($\tau_\infty = 1$ resolved on the +fault plane), and the whole constitutive ladder — the slip-rate profile +$V(s)$ tells each law's story: + +![The fault-strength ladder](figures/fault-examples/ladder.png) + +The right-hand column shows the shear-stress field $\sigma_{xy}$ for +each rung (colour range $0 \to 2\tau_\infty$, white at the far-field +value): a slipping fault shadows itself with a stress drop and +concentrates stress at its tips, in proportion to how much of +$\tau_\infty$ its law lets go; the stuck fault leaves the field +untouched. + +- **Frictionless** — the free crack: full stress drop, elliptical + profile (the dashed shape), slip vanishing at the unsplit tips. +- **Viscous, $\eta_f = \eta/a$** — the natural interface-viscosity + scale: the fault slips at roughly half its free rate. +- **Coulomb, weak** ($\mu\sigma_n < \tau_\infty$) — slides at the + reduced stress drop $\tau_\infty - \mu\sigma_n$, profile still + elliptical: the constant-strength crack. +- **Rate-and-state** — at steady state, a velocity-dependent strength + between the Coulomb end-members. +- **Coulomb, strong** ($\mu\sigma_n > \tau_\infty$) — stuck: creep at + the regularisation velocity, invisible at plot scale. + +Script: `figures/fault-examples/ladder.py`. + +## The Mohr circle, measured by faults + +All the stress-plane figures below use the **geological sign +convention**: compression positive, tension on the negative axis. (The +solver's tractions are tension-positive; only the plots flip.) There is +no confining pressure in these problems: the walls prescribe velocity +and the flow is incompressible, so pressure exists only up to the +nullspace constant and the solver's gauge centres the Mohr circle on +the origin — which is exactly why the tensile sector is reachable. + +A *welded* split fault (interface dashpot with large $\eta_f$) does not +perturb the flow — the welded limit recovers the uncut continuum — but +its machinery still reports the tractions across it: the no-opening +constraint's reaction gives $\sigma_n$, and the dashpot law reads the +signed shear traction off its own slip, $\tau = \eta_f V$. Each fault +orientation is therefore a passive stress probe, and sweeping the +orientation through $180°$ traces the full Mohr circle of the ambient +stress state: + +![Mohr circle from welded faults](figures/fault-examples/mohr-circle.png) + +The boundary condition is Dirichlet velocity on all four walls, +imposing the homogeneous flow $\mathbf v = (a(x-\tfrac12) + +\gamma(y-\tfrac12),\; -a(y-\tfrac12))$ with $a = 0.5$, $\gamma = 1$ +— a pure-shear (stretching) part plus a simple-shear part. Stress sees +only the symmetric gradient, so this is equivalent to pure shear of +magnitude $R = \eta\sqrt{4a^2+\gamma^2}$ with principal axes at +$22.5°$ to the box: deliberately not axis-aligned, so the circle's +orientation must be measured, not guessed. Because the linear flow is +an exact homogeneous Stokes solution, the stress is uniform to the +walls and every probe samples the same state; the pressure gauge (the +nullspace constant) sets the circle's centre. The probes (labelled by +fault angle) +land on that circle, centred at the (gauge-fixed) mean pressure. The +classical double-angle rule is the labels' spacing: $22.5°$ of fault +rotation moves a probe $45°$ around the circle, and $180°$ of rotation +closes it. + +Script: `figures/fault-examples/mohr_circle.py`. + +The measured points land on the circle too neatly to teach from a +static figure — the lesson is in the *construction*. The animated +version rotates the fault through $180°$ while its probe sweeps the +full circle at $2\theta$, one welded-fault solve per frame: + +![Building the Mohr circle](figures/fault-examples/mohr-circle-build.gif) + +Script: `figures/fault-examples/mohr_animate.py`. + +### The circle meets the friction envelope + +Give the rotating fault Coulomb friction (reaction-fed normal stress) +and the probe can no longer go everywhere the ambient stress points. +Three regimes appear in one sweep: + +![Coulomb fault vs the Mohr circle](figures/fault-examples/mohr-friction.png) + +- **Stuck** (circles): the ambient resolved stress lies inside the + envelope $|\tau| < \mu\sigma$ — the fault transmits it and the + probe sits on the Mohr circle; +- **Sliding** (squares): the ambient stress would exceed the envelope — + the fault slips, drops the shear traction to its strength, and the + probe is pinned to the yield line; +- **Held shut** ($\sigma < 0$, grey crosses): bare friction has no + strength under tension — a real fault would OPEN, and no static + solution exists. The bilateral no-opening constraint manufactures one + by gluing the surfaces together (a tensile reaction), so the solver + still converges; the probes ride the axis at $\tau = 0$, but they + are marked unphysical. The tell, in any real model, is the sign of + the recovered normal traction. + +The shear traction here is read from the Coulomb law at the measured +slip rate — which is exact in every regime, because the regularised +law *is* the traction the fault carries. The animated build shows the +switch happening as the fault rotates, with the slip sense drawn as +half-arrows when it unlocks: + +![A frictional fault against the Mohr circle](figures/fault-examples/mohr-friction-build.gif) + +Script: `figures/fault-examples/mohr_friction.py`. + +### Cohesion keeps more of the circle + +Add cohesion — the envelope becomes $\tau = \pm(C + \mu\sigma)$ — +and strength survives into MILD tension, declining along the envelope +and reaching zero at $\sigma = -C/\mu$. Stuck arcs now survive +around **both** principal poles, with envelope-pinned sliding between; +beyond the cutoff the fault is in the same held-shut (unphysical) +regime as the cohesionless case, just pushed further into tension. + +![Cohesive Mohr-Coulomb fault](figures/fault-examples/mohr-cohesion.png) + +The cohesive law is not a canned option — it is registered as a sympy +expression in the canonical interface symbols, which is the intended +extension path for any fault rheology: + +```python +V = fault_contact.slip_rate +S = fault_contact.normal_stress # reaction-fed, clamped >= 0 +law = fault_contact.SymbolicFaultLaw( + (C + mu * S) * (2 / sympy.pi) * sympy.atan(V / V0)) +``` + +![Cohesion keeps more of the circle](figures/fault-examples/mohr-cohesion-build.gif) + +Script: `figures/fault-examples/mohr_cohesion.py`. + +### The graded fault: depth-dependent stress along one surface + +Add a modest hydrostatic load (constant density, closed box — the flow +is untouched, pressure absorbs gravity exactly) and the per-node +character of the stress recovery becomes visible: every fault node sits +at its own depth, so a single welded fault contributes a **streak** of +probes, not a point. The streak is horizontal (only the pressure part +of $\sigma$ varies along the fault), longest for the vertical fault, +and collapses to a dot for the horizontal one: + +![The graded fault](figures/fault-examples/mohr-graded.png) + +Points are coloured by depth; the grey circles are the Mohr circles of +the shallowest, central, and deepest fault points. Nothing here is +averaged — $\sigma_n$ comes from the constraint reaction de-smeared +node by node, $\tau$ from the weld's own law at each node — which is +exactly the machinery that lets a friction law feel depth-dependent +strength along a single fault (the locking-depth structure of a +seismogenic zone). The closed-box pressure gauge is re-anchored in the +plot so $p = 0$ at the top surface; the shift is the analytically known +$\rho g H/2$. + +Script: `figures/fault-examples/mohr_graded.py`. + +## Orientation and slip: the circle's other face + +The same orientation sweep with *frictionless* faults under pure shear. +Now each fault drops whatever shear stress is resolved on its plane, so +the peak slip rate follows $|\cos 2\theta|$ — the slip-rate reading of +the Mohr circle. A fault at $45°$ to the shear plane feels no resolved +shear and barely slips; the aligned fault slips fully: + +![Orientation vs slip](figures/fault-examples/orientations.png) + +Script: `figures/fault-examples/orientations.py`. + +## Interacting faults: stress transfer on the Mohr plane + +Two faults via `add_fault`'s network form: the *source* slips freely, +the *receiver* is welded — the passive probe from the Mohr examples. +Differencing two solves on the same mesh gives the classical Coulomb +stress transfer, and the receiver's probes show it in the diagram +students already know: each node's point MOVES, and its colour is its +own $\Delta$CFF ($\mu' = 0.4$, King's value): + +![En echelon interaction](figures/fault-examples/interacting-faults.png) + +The receiver sits in the source's along-strike tip lobe, and the +friction bookkeeping is dressed so the crossing is visible: a declared +confining pressure $P_0 = 1$ puts the whole circle in compression, and +a cohesive envelope $\tau = \pm(C + \mu'\sigma)$ with $C = 0.75$ +passes just below the ambient cloud. Neither constant changes +$\Delta$CFF — they only place the failure line where the physics can +reach it. The loaded near-tip nodes then cross into the shaded failure +region while the far end barely moves. Rotating the regional +compression axis (the boundary velocities rotate with it) relocates +the ambient cloud around the circle, and what rotation controls is the +MARGIN to the envelope. (Fields here and below are P0 cell stress on +the split mesh's true connectivity, linear colour at $\pm 1$ — see +the California example for why.) + +![Rotating the regional stress](figures/fault-examples/interacting-rotation.png) + +Two gauge points, handled explicitly in the scripts because a closed +velocity-driven box fixes pressure only up to a constant: differenced +fields are anchored to zero in the far field (a slip event changes +nothing far away), and each solve's absolute probe pressures are +anchored to the analytic ambient state. Both constants are printed, +never silently absorbed. + +Script: `figures/fault-examples/interacting_faults.py`. + +### A schematic southern California + +The San Andreas as ONE continuous dextral trace with a smooth tanh +S-bend — the Big Bend, the smoothed version of a left-stepping +stepover, restraining under right-lateral slip. The Garlock, three +East California Shear Zone strands and a San Jacinto-like fault are +welded as probes — all inboard, in the continental crust. The drive is +right-lateral simple shear parallel to the plate boundary (~N40W), +which resolves dextral on the NW-striking faults and sinistral on the +Garlock, exactly the real senses; the slip arrows on the map are drawn +from the *measured* jump, not assumed: + +![Schematic southern California](figures/fault-examples/california.png) + +The curved trace is sampled as a polyline, and the fault carries the +smooth curve's **analytic normal** +(`add_fault_bc(..., normal=...)`). This matters: with the default +facet-averaged normals a sampled curve develops slip notches and +stress sawteeth at every sampling kink — roughness that *grows* under +refinement, because the constraint direction itself zig-zags (see the +practical notes in the split-node fault guide). With the analytic +normal the trace slips smoothly through the bend. + +One trunk-spanning event, and the map teaches by geography: the +restraining bend fills with strong compression — $\Delta$CFF deeply +negative in a bowtie exactly where the Transverse Ranges belong. The +neighbours give three verdicts: the parallel San Jacinto is deeply +relaxed ($\Delta$CFF $-0.4$, its cloud retreating into the safe +wedge); the conjugate Garlock and the distant ECSZ strands are mildly +loaded ($+0.05$ to $+0.06$). (Schematic geometry, not to scale; same +$P_0 = 1$, $C = 0.75$ envelope dressing.) + +Script: `figures/fault-examples/california.py`. + +### A branching rupture on an intersecting network + +Real trace sets cross and branch, and the split-node pipeline refuses +shared vertices — so intersecting geometry goes through +`uw.meshing.prepare_fault_network` first, which converts every +junction (X crossing, T abutment, near-miss) to the offset-junction +form: the traces stop short of the intersection, leaving a ligament of +intact material about a cell across, with angle-corrected pull-backs +and a printed report of every action. Here the raw input genuinely +intersects — a dextral trunk, a splay branching off its midpoint, and +a conjugate crossing it — and the trunk + splay rupture together while +the conjugate sits welded: + +![Branching rupture](figures/fault-examples/branching.png) + +The brightest features are the **junction plugs**: the intact +ligaments take the concentrated load shed by their slipped neighbours +on both sides, which is exactly where a through-going event would +break next. The ligament approximation is measured, not asserted: at a +T-junction the branch's slip profile is *identical* at two mesh +resolutions for a fixed plug size, and its peak slip changes by only a +few percent between ligaments of one and two cells (the study is +recorded in the developer notes). Keep ligaments at one to two local +cell sizes. + +Script: `figures/fault-examples/branching.py`. + +And the question the offset form begs — what does interrupting a fault +cost? — has a measured answer. The same network prepared two ways +(`through=["Trunk"]` keeps the trunk continuous, so the conjugate +yields on both sides of the crossing; the default cuts both): + +![Continuous vs interrupted](figures/fault-examples/branching-compare.png) + +The through-going trunk slips slightly *more* than it would alone (the +slipping splay feeds it); cutting it at the crossing pins it locally +and costs about a quarter of its peak slip — but the profiles converge +away from the junction, so an interruption's reach is the segment +scale, not the system scale. T abutments never cut the through-going +trace; declare master faults with `through=` where a crossing must not +either. + +Script: `figures/fault-examples/branching_compare.py`. + +### A true Y-branch, and which continuity to declare + +A genuine three-arm branch can be decomposed two ways: trunk +continuous with the splay abutting (A), or the *bent* west-arm-plus- +splay fault continuous with the east arm abutting (B). Running both, +with the gap swept in A: + +![The true branch, bracketed](figures/fault-examples/true-branch.png) + +Two conclusions. The **gap is second-order**: A at 1h, 2h and 4h +ligaments overlays on every arm except the splay's near-junction toe, +which simply moves with the trimmed tip. The **decomposition is +first-order**: B is a different mechanical system, because the +continuous bent fault LOCKS at its 33-degree kink — slip through a +kink is geometrically incompatible with the no-opening constraint +(the same mechanism that makes sampled curves rough, here as genuine +physics) — so both trunk arms lose a third of their slip and pinch to +zero at the branch point. The rule for real networks: **declare the +straightest path continuous; never route continuity around a +corner.** The energetically sensible decomposition (A) is admissible +for the true junction and differs from it only in the splay's +ligament-scale toe. + +The stress fields make the distinction physical. With the trunk +continuous, the stress shadow runs straight through the junction and +the branch point carries no feature at all — the junction is +mechanically invisible. With continuity routed around the bend, the +locked kink acts as a barrier: both arms' shadows terminate against it +and load the corner, the classic signature of a rupture-arresting +bend: + +![The branch point in stress](figures/fault-examples/true-branch-stress.png) + +Script: `figures/fault-examples/true_branch.py`. + +## What these models are, and are not + +For the static problem, incompressible linear elasticity and Stokes +are the same mathematics (displacement ↔ velocity, shear modulus ↔ +viscosity), so every stress pattern here is exactly the elastic +$\nu = 1/2$ pattern with slip rate standing in for coseismic slip. +What is deliberately absent: depth and half-space geometry (these are +plane-strain boxes with velocity-driven walls — no free-surface +amplification), topography, gravity except where stated (the graded +fault), compressibility ($\nu \ne 1/2$ shifts normal-stress lobes +modestly), and postseismic processes of any kind. The slipping faults +are completely weak during their event (frictionless: the full +ambient resolved shear drops); because the problem is linear, +$\Delta$CFF scales essentially linearly with the dropped stress, so +the transfer maps read as stress change per unit stress drop. A +partial-drop event is the Coulomb rung applied to the source fault. + +## Where next + +Remaining extensions on the same harness: tip-to-tip vs overlapped +en echelon arrangements side by side, a denser fault network, and the +graded (gravity-loaded) versions of the interaction cases. diff --git a/docs/advanced/figures/fault-examples/.gitignore b/docs/advanced/figures/fault-examples/.gitignore new file mode 100644 index 000000000..04a1f04a5 --- /dev/null +++ b/docs/advanced/figures/fault-examples/.gitignore @@ -0,0 +1,2 @@ +_*.png +*.log diff --git a/docs/advanced/figures/fault-examples/_california_probes.npz b/docs/advanced/figures/fault-examples/_california_probes.npz new file mode 100644 index 000000000..6c0bef08f Binary files /dev/null and b/docs/advanced/figures/fault-examples/_california_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/_interacting_probes.npz b/docs/advanced/figures/fault-examples/_interacting_probes.npz new file mode 100644 index 000000000..e0a3b4c70 Binary files /dev/null and b/docs/advanced/figures/fault-examples/_interacting_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/_mohr_cohesion_probes.npz b/docs/advanced/figures/fault-examples/_mohr_cohesion_probes.npz new file mode 100644 index 000000000..85e285682 Binary files /dev/null and b/docs/advanced/figures/fault-examples/_mohr_cohesion_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/_mohr_friction_probes.npz b/docs/advanced/figures/fault-examples/_mohr_friction_probes.npz new file mode 100644 index 000000000..d6fd0213a Binary files /dev/null and b/docs/advanced/figures/fault-examples/_mohr_friction_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/_mohr_graded_probes.npz b/docs/advanced/figures/fault-examples/_mohr_graded_probes.npz new file mode 100644 index 000000000..95ff368ed Binary files /dev/null and b/docs/advanced/figures/fault-examples/_mohr_graded_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/_mohr_probes.npz b/docs/advanced/figures/fault-examples/_mohr_probes.npz new file mode 100644 index 000000000..129bd50ff Binary files /dev/null and b/docs/advanced/figures/fault-examples/_mohr_probes.npz differ diff --git a/docs/advanced/figures/fault-examples/branching-compare.png b/docs/advanced/figures/fault-examples/branching-compare.png new file mode 100644 index 000000000..58deb39fa Binary files /dev/null and b/docs/advanced/figures/fault-examples/branching-compare.png differ diff --git a/docs/advanced/figures/fault-examples/branching.png b/docs/advanced/figures/fault-examples/branching.png new file mode 100644 index 000000000..67c4ca9c3 Binary files /dev/null and b/docs/advanced/figures/fault-examples/branching.png differ diff --git a/docs/advanced/figures/fault-examples/branching.py b/docs/advanced/figures/fault-examples/branching.py new file mode 100644 index 000000000..f55ac9405 --- /dev/null +++ b/docs/advanced/figures/fault-examples/branching.py @@ -0,0 +1,178 @@ +"""A branching rupture: intersecting fault traces and their Delta CFF. + +The raw geometry INTERSECTS: a dextral trunk, a splay branching off its +midpoint at ~30 degrees (a T junction), and a conjugate fault crossing +it outright (an X junction). ``prepare_fault_network`` converts both to +offset-junction form automatically (angle-corrected ligaments, loudly), +which is what makes the set splittable at all. The trunk and the splay +then rupture TOGETHER (frictionless) while the conjugate is welded as a +receiver, and the map shows Delta CFF on trunk-parallel planes — with +zooms at the two junctions, where the ligament-scale stress transfer +lives. Measured ligament sensitivity for this representation: +~/+Simulations/fault_junction_ligament/ (the branch response is +converged for a given plug size; ligaments of 1-2 local cells). +""" +import os +import time + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv + +import underworld3 as uw +from underworld3.meshing.surfaces import prepare_fault_network +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +MU_P = 0.4 +H = 0.012 +LIG_F = 1.5 + +# raw, INTERSECTING traces — the preparer makes them legal +TREND = np.degrees(np.arctan2(0.10, 0.70)) # trunk trend ~8 deg +M_RAW = ("Trunk", np.array([[0.15, 0.45], [0.85, 0.55]])) +S_RAW = ("Splay", np.array([[0.50, 0.50], [0.80, 0.76]])) # T at trunk +C_RAW = ("Conj", np.array([[0.30, 0.33], [0.42, 0.64]])) # X crossing + +prepared, report = prepare_fault_network( + [M_RAW, S_RAW, C_RAW], spacing=H, ligament=LIG_F, verbose=True) +names = [n for n, _ in prepared] +SLIP = [n for n in names if n.startswith(("Trunk", "Splay"))] +WELD = [n for n in names if n.startswith("Conj")] +ETA_WELD = 200.0 * common.ETA / 0.2 + +T_J = np.array([0.50, 0.50]) # the T junction +X_J = M_RAW[1][0] + (M_RAW[1][1] - M_RAW[1][0]) * 0.293 # approx X point + + +def solve_state(child, free): + stokes = common.stokes_on(child, + common.boundary_simple_shear(child, TREND)) + for n in names: + stokes.add_fault_bc(0.0 if (free and n in SLIP) else ETA_WELD, + boundary=n) + fault_contact.solve_with_fault(stokes, picard=2) + x, y = child.X + v, p = stokes.Unknowns.u, stokes.Unknowns.p + comps = {} + for cname, expr in ( + ("sxx", -p.sym[0] + 2 * common.ETA * v.sym[0].diff(x)), + ("syy", -p.sym[0] + 2 * common.ETA * v.sym[1].diff(y)), + ("sxy", common.ETA * (v.sym[0].diff(y) + v.sym[1].diff(x)))): + s_var = uw.discretisation.MeshVariable( + f"{cname}_{'a' if free else 'b'}", child, 1, degree=0, + continuous=False) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = expr + proj.smoothing = 0.0 + proj.solve() + row = common.split_mesh_cell_rows(child, s_var) + comps[cname] = np.asarray(s_var.data[:, 0])[row].copy() + return stokes, comps + + +t0 = time.perf_counter() +child = common.base_mesh(H).add_fault(prepared) +s1, c1 = solve_state(child, free=True) +print(f"[timing] slipping solve: {time.perf_counter() - t0:.1f} s", + flush=True) +for n in SLIP: + coords, jumps, normals = fault_contact.fault_pair_jumps( + s1, n, s1._rotated_freeslip_info) + tang = np.column_stack([-normals[:, 1], normals[:, 0]]) + V = np.einsum("ij,ij->i", jumps, tang) + print(f" {n:8s} peak slip {np.abs(V).max():.4f}", flush=True) +t0 = time.perf_counter() +_s0, c0 = solve_state(child, free=False) +print(f"[timing] welded solve: {time.perf_counter() - t0:.1f} s", + flush=True) + +# Delta CFF on trunk-parallel receiver planes +beta = np.radians(TREND) +nx, ny = -np.sin(beta), np.cos(beta) +tx, ty = np.cos(beta), np.sin(beta) + + +def resolve(c): + s_nn = c["sxx"] * nx * nx + 2 * c["sxy"] * nx * ny + c["syy"] * ny * ny + s_t = (c["sxx"] * tx * nx + c["sxy"] * (tx * ny + ty * nx) + + c["syy"] * ty * ny) + return s_nn, s_t + + +nn0, t_0 = resolve(c0) +nn1, t_1 = resolve(c1) +tau_dir = np.sign(np.median(t_0)) +dcff = tau_dir * (t_1 - t_0) + MU_P * (nn1 - nn0) + +pts, faces = common.split_mesh_cell_render(child) +fc = np.asarray(faces).reshape(-1, 4)[:, 1:] +cent = np.asarray(pts)[fc].mean(axis=1) +dcff, gauge = common.far_field_anchor( + cent, dcff, [p for _n, p in prepared], cut=0.18) +print(f"far-field gauge removed: {gauge:+.4f}", flush=True) + +# ---- renders: the map and the two junction zooms --------------------------- +COLOUR = {"Trunk": "black", "Splay": "black", "Conj": "#6a1b9a"} + + +def render(png, scale, focal): + pvm = pv.PolyData(np.asarray(pts, dtype=float), + faces=np.asarray(faces, dtype=np.int64)) + pvm.cell_data["dcff"] = dcff + pl = pv.Plotter(off_screen=True, window_size=(900, 850)) + pl.set_background("white") + pl.add_mesh(pvm, scalars="dcff", cmap="RdBu_r", clim=(-1.0, 1.0), + lighting=False, show_scalar_bar=False) + for n, p in prepared: + line = pv.lines_from_points( + np.column_stack([p, np.full(len(p), 1e-3)])) + pl.add_mesh(line, color=COLOUR[n.split("_")[0]], + line_width=5 if n.startswith("Trunk") else 4, + lighting=False) + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.parallel_scale = scale + pl.camera.focal_point = (focal[0], focal[1], 0.0) + pl.screenshot(png) + pl.close() + return png + + +map_png = render(os.path.join(D, "_branching_map.png"), 0.42, (0.5, 0.52)) +tz_png = render(os.path.join(D, "_branching_tzoom.png"), 0.075, T_J) +xz_png = render(os.path.join(D, "_branching_xzoom.png"), 0.075, X_J) + +fig = plt.figure(figsize=(12.6, 6.4)) +gs = fig.add_gridspec(2, 3, width_ratios=[2.1, 1.0, 0.05]) +axm = fig.add_subplot(gs[:, 0]) +axm.imshow(plt.imread(map_png)) +axm.set_xticks([]) +axm.set_yticks([]) +axm.set_title(r"$\Delta$CFF on trunk-parallel planes ($\mu' = 0.4$): " + "trunk + splay rupture together,\nconjugate welded; " + "junctions are auto-converted offset ligaments", fontsize=9.5) +for row, (png, label) in enumerate( + ((tz_png, "the branch point (T): splay fed through the ligament"), + (xz_png, "the crossing (X): four tips, one intact plug"))): + ax = fig.add_subplot(gs[row, 1]) + ax.imshow(plt.imread(png)) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_title(label, fontsize=8.5) + +from matplotlib import cm, colors as mcolors +sm = cm.ScalarMappable(norm=mcolors.Normalize(-1, 1), cmap="RdBu_r") +cax = fig.add_subplot(gs[:, 2]) +fig.colorbar(sm, cax=cax, label=r"$\Delta$CFF (per unit stress drop)") +fig.suptitle("A branching rupture on an intersecting fault network", + fontsize=11.5) +fig.tight_layout() +out = os.path.join(D, "branching.png") +fig.savefig(out, dpi=200) +print("wrote", out, flush=True) diff --git a/docs/advanced/figures/fault-examples/branching_compare.py b/docs/advanced/figures/fault-examples/branching_compare.py new file mode 100644 index 000000000..6f59e7d72 --- /dev/null +++ b/docs/advanced/figures/fault-examples/branching_compare.py @@ -0,0 +1,205 @@ +"""What does interrupting the trunk cost? Continuous vs cut at the +junction, against the isolated trunk. + +Same network as ``branching.py`` (trunk + splay slipping, conjugate +welded), prepared two ways: + +- ``through=["Trunk"]`` — the trunk is the MASTER: never cut; the + conjugate yields on both sides of the crossing and the splay is + trimmed where it abuts (T junctions never cut the through-going + trace). +- default — the X crossing cuts both traces, so the trunk is + interrupted mid-length by the junction plug. + +The isolated continuous trunk (no other faults) is the reference +profile. The comparison answers: how much slip does the interruption +forfeit, and what does it do to the Delta CFF field? +""" +import os +import time + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv + +import underworld3 as uw +from underworld3.meshing.surfaces import prepare_fault_network +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +MU_P = 0.4 +H = 0.012 +LIG_F = 1.5 +TREND = np.degrees(np.arctan2(0.10, 0.70)) +M_RAW = ("Trunk", np.array([[0.15, 0.45], [0.85, 0.55]])) +S_RAW = ("Splay", np.array([[0.50, 0.50], [0.80, 0.76]])) +C_RAW = ("Conj", np.array([[0.30, 0.33], [0.42, 0.64]])) +T_HAT = (M_RAW[1][1] - M_RAW[1][0]) +T_HAT = T_HAT / np.linalg.norm(T_HAT) +ETA_WELD = 200.0 * common.ETA / 0.2 + + +def run_case(tag, faults, weld_names, want_field): + child = common.base_mesh(H).add_fault(faults) + names = [n for n, _ in faults] + + def one_state(free): + stokes = common.stokes_on( + child, common.boundary_simple_shear(child, TREND)) + for n in names: + stokes.add_fault_bc( + 0.0 if (free and n not in weld_names) else ETA_WELD, + boundary=n) + fault_contact.solve_with_fault(stokes, picard=2) + return stokes + + t0 = time.perf_counter() + s1 = one_state(True) + # trunk slip vs GLOBAL trunk arc length, across however many pieces + ss, vv = [], [] + for n in names: + if not n.startswith("Trunk"): + continue + coords, jumps, normals = fault_contact.fault_pair_jumps( + s1, n, s1._rotated_freeslip_info) + tang = np.column_stack([-normals[:, 1], normals[:, 0]]) + ss.append((coords - M_RAW[1][0]) @ T_HAT) + vv.append(np.abs(np.einsum("ij,ij->i", jumps, tang))) + s_all = np.concatenate(ss) + v_all = np.concatenate(vv) + order = np.argsort(s_all) + profile = (s_all[order], v_all[order]) + + dcff_pack = None + if want_field: + def stress(stokes, tagc): + x, y = child.X + v, p = stokes.Unknowns.u, stokes.Unknowns.p + comps = {} + for cname, expr in ( + ("sxx", -p.sym[0] + 2 * common.ETA * v.sym[0].diff(x)), + ("syy", -p.sym[0] + 2 * common.ETA * v.sym[1].diff(y)), + ("sxy", common.ETA * (v.sym[0].diff(y) + + v.sym[1].diff(x)))): + s_var = uw.discretisation.MeshVariable( + f"{cname}_{tagc}", child, 1, degree=0, + continuous=False) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = expr + proj.smoothing = 0.0 + proj.solve() + row = common.split_mesh_cell_rows(child, s_var) + comps[cname] = np.asarray(s_var.data[:, 0])[row].copy() + return comps + + c1 = stress(s1, "a") + s0 = one_state(False) + c0 = stress(s0, "b") + beta = np.radians(TREND) + nx, ny = -np.sin(beta), np.cos(beta) + tx, ty = np.cos(beta), np.sin(beta) + + def resolve(c): + s_nn = (c["sxx"] * nx * nx + 2 * c["sxy"] * nx * ny + + c["syy"] * ny * ny) + s_t = (c["sxx"] * tx * nx + c["sxy"] * (tx * ny + ty * nx) + + c["syy"] * ty * ny) + return s_nn, s_t + + nn0, t_0 = resolve(c0) + nn1, t_1 = resolve(c1) + dcff = np.sign(np.median(t_0)) * (t_1 - t_0) + MU_P * (nn1 - nn0) + pts, faces = common.split_mesh_cell_render(child) + fc = np.asarray(faces).reshape(-1, 4)[:, 1:] + cent = np.asarray(pts)[fc].mean(axis=1) + dcff, gauge = common.far_field_anchor( + cent, dcff, [p for _n, p in faults], cut=0.18) + print(f"[{tag}] gauge {gauge:+.4f}", flush=True) + dcff_pack = (pts, faces, dcff) + + print(f"[{tag}] trunk peak slip {profile[1].max():.4f} " + f"({time.perf_counter() - t0:.1f} s)", flush=True) + return profile, dcff_pack, faults + + +# case i: trunk through-going (master) +prep_i, _ = prepare_fault_network( + [(n, p.copy()) for n, p in (M_RAW, S_RAW, C_RAW)], spacing=H, + ligament=LIG_F, through=["Trunk"], verbose=True) +weld_i = [n for n, _p in prep_i if n.startswith("Conj")] +prof_i, field_i, faults_i = run_case("through", prep_i, weld_i, True) + +# case ii: default — the crossing cuts the trunk +prep_ii, _ = prepare_fault_network( + [(n, p.copy()) for n, p in (M_RAW, S_RAW, C_RAW)], spacing=H, + ligament=LIG_F, verbose=True) +weld_ii = [n for n, _p in prep_ii if n.startswith("Conj")] +prof_ii, field_ii, faults_ii = run_case("cut", prep_ii, weld_ii, True) + +# reference: the trunk alone +prof_ref, _, _ = run_case("isolated", [("Trunk", M_RAW[1].copy())], + [], False) + + +def render(pack, faults, png): + pts, faces, dcff = pack + pvm = pv.PolyData(np.asarray(pts, dtype=float), + faces=np.asarray(faces, dtype=np.int64)) + pvm.cell_data["dcff"] = dcff + pl = pv.Plotter(off_screen=True, window_size=(900, 850)) + pl.set_background("white") + pl.add_mesh(pvm, scalars="dcff", cmap="RdBu_r", clim=(-1.0, 1.0), + lighting=False, show_scalar_bar=False) + for n, p in faults: + line = pv.lines_from_points( + np.column_stack([p, np.full(len(p), 1e-3)])) + pl.add_mesh(line, color=("#6a1b9a" if n.startswith("Conj") + else "black"), + line_width=5 if n.startswith("Trunk") else 4, + lighting=False) + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.parallel_scale = 0.42 + pl.camera.focal_point = (0.5, 0.52, 0.0) + pl.screenshot(png) + pl.close() + return png + + +png_i = render(field_i, faults_i, os.path.join(D, "_bcmp_through.png")) +png_ii = render(field_ii, faults_ii, os.path.join(D, "_bcmp_cut.png")) + +fig = plt.figure(figsize=(13.0, 5.4)) +gs = fig.add_gridspec(1, 4, width_ratios=[1.35, 1.0, 1.0, 0.05]) +axp = fig.add_subplot(gs[0, 0]) +axp.plot(*prof_ref, "k--", lw=1.3, label="trunk alone (continuous)") +axp.plot(*prof_i, "o-", ms=2.5, lw=1.0, color="#1565c0", + label="network, trunk through-going") +axp.plot(*prof_ii, "o-", ms=2.5, lw=1.0, color="#e65100", + label="network, trunk cut at the crossing") +axp.set_xlabel("arc length along the trunk") +axp.set_ylabel("|slip|") +axp.set_title("what the interruption costs the trunk", fontsize=10) +axp.legend(fontsize=8) +for png, col, title in ((png_i, 1, "trunk through-going"), + (png_ii, 2, "trunk cut at the crossing")): + ax = fig.add_subplot(gs[0, col]) + ax.imshow(plt.imread(png)) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_title(title, fontsize=10) +from matplotlib import cm, colors as mcolors +sm = cm.ScalarMappable(norm=mcolors.Normalize(-1, 1), cmap="RdBu_r") +cax = fig.add_subplot(gs[0, 3]) +fig.colorbar(sm, cax=cax, label=r"$\Delta$CFF") +fig.suptitle("Continuous vs interrupted: the same network, two " + "junction policies", fontsize=11.5) +fig.tight_layout() +out = os.path.join(D, "branching-compare.png") +fig.savefig(out, dpi=200) +print("wrote", out, flush=True) diff --git a/docs/advanced/figures/fault-examples/california.png b/docs/advanced/figures/fault-examples/california.png new file mode 100644 index 000000000..cec34eec6 Binary files /dev/null and b/docs/advanced/figures/fault-examples/california.png differ diff --git a/docs/advanced/figures/fault-examples/california.py b/docs/advanced/figures/fault-examples/california.py new file mode 100644 index 000000000..953714bcd --- /dev/null +++ b/docs/advanced/figures/fault-examples/california.py @@ -0,0 +1,334 @@ +"""A schematic southern California: the San Andreas with a smooth +restraining bend (the Big Bend), read by its neighbours. + +Geography (schematic, not to scale, x = east / y = north): the SAN +ANDREAS is ONE continuous dextral trace with a smooth tanh S-bend — the +smoothed version of a left-stepping stepover. Under right-lateral shear +the bend is RESTRAINING: the bend zone is where the TRANSVERSE RANGES +belong, right beside the Garlock. The GARLOCK (resolving sinistral — +the real sense, from the kinematics), three EAST CALIFORNIA SHEAR ZONE +strands and a SAN JACINTO-like fault sit inboard as welded probes. + +The curved trace is sampled as a polyline, and the fault's constraint +frame uses the curve's ANALYTIC normal (``add_fault_bc(normal=...)``). +Without it the per-node normal averages the adjacent facet normals and +zig-zags at the sampling kinks; the no-opening constraint then forbids +smooth slip past each kink — sawtooth tractions that GROW under mesh +refinement (measured: ~/+Simulations/curved_fault_roughness/). With the +analytic normal the polyline behaves as the smooth fault it represents. + +The San Andreas slips freely under right-lateral simple shear parallel +to the plate-boundary trend; the slip half-arrows on the map come from +the MEASURED jump. Field: Delta CFF on boundary-parallel planes, P0 +(cell) stress — continuous-P1 projection of the rough near-fault stress +rings at the node scale (measured) — rendered on the split mesh's true +connectivity. +""" +import os +import time + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import sympy +import pyvista as pv + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +MU_P = 0.4 +TAU0 = 1.0 +TREND = 132.0 # plate-boundary trend (~N42W), degrees CCW of E +ETA_WELD = 200.0 * common.ETA / 0.2 +P0 = 1.0 +COH = 0.75 + +# One continuous dextral trace with a smooth restraining bend: the +# centreline runs along the trend from A, offset by a tanh step of +# BEND_W toward the SW (the CCW normal) — an S in map view, stepping +# LEFT, so right-lateral slip must converge through the bend. Max +# deviation from the trend: atan(BEND_W / (2 LAM)) ~ 30 degrees, about +# the real Big Bend. +_A = np.array([0.88, 0.06]) +_t = np.array([np.cos(np.radians(TREND)), np.sin(np.radians(TREND))]) +_n = np.array([-_t[1], _t[0]]) # CCW normal (points SW) +S_END = 0.94 +BEND_W = 0.07 # total SW offset +BEND_S0 = 0.45 # bend centre, arc parameter +LAM = 0.06 # bend half-width + + +def _w(s): + return 0.5 * BEND_W * (1.0 + np.tanh((s - BEND_S0) / LAM)) + + +def saf_trace(n_seg=47): + """The smooth trace sampled as a polyline (kinks land on mesh + vertices; the analytic normal makes them harmless).""" + s = np.linspace(0.0, S_END, n_seg + 1) + return _A + np.outer(s, _t) + np.outer(_w(s), _n) + + +def saf_normal(child): + """The EXACT unit-normal direction of the smooth trace as a sympy + row matrix in mesh coordinates: with X(s) = A + s t + w(s) n and + s(X) = (X - A).t, the tangent is t + w'(s) n and the normal its + quarter turn, n - w'(s) t (normalisation is the caller's).""" + x, y = child.X + s = (x - _A[0]) * _t[0] + (y - _A[1]) * _t[1] + wp = (0.5 * BEND_W / LAM) * (1 - sympy.tanh((s - BEND_S0) / LAM) ** 2) + return sympy.Matrix([[_n[0] - wp * _t[0], _n[1] - wp * _t[1]]]) + + +SAF_PTS = saf_trace() + +MINORS = { + "Garlock": np.array([[0.52, 0.44], [0.80, 0.53]]), + "E1": np.array([[0.62, 0.64], [0.56, 0.78]]), + "E2": np.array([[0.72, 0.62], [0.66, 0.76]]), + "E3": np.array([[0.81, 0.68], [0.75, 0.82]]), + "SJF": np.array([[0.84, 0.24], [0.70, 0.38]]), +} +GROUPS = (("Garlock (sinistral)", ("Garlock",), "#6a1b9a"), + ("ECSZ, 3 strands (dextral)", ("E1", "E2", "E3"), "#1a6b1a"), + ("San Jacinto (dextral)", ("SJF",), "#e65100")) +COLOUR = {"Garlock": "#6a1b9a", "E1": "#1a6b1a", "E2": "#1a6b1a", + "E3": "#1a6b1a", "SJF": "#e65100"} + + +def build_and_solve(trunk_free): + faults = [("SAF", SAF_PTS)] + [(k, v) for k, v in MINORS.items()] + child = common.base_mesh(0.012).add_fault(faults) + stokes = common.stokes_on(child, + common.boundary_simple_shear(child, TREND, + TAU0)) + # the analytic normal in BOTH states, so the slipping and welded + # solves share one constraint frame and difference cleanly + stokes.add_fault_bc(0 if trunk_free else ETA_WELD, boundary="SAF", + normal=saf_normal(child)) + for k in MINORS: + stokes.add_fault_bc(ETA_WELD, boundary=k) + fault_contact.solve_with_fault(stokes, picard=2) + probes = {} + for k, pts in MINORS.items(): + _s, _xy, sig, tau = common.probe_nodes(stokes, k, pts[1] - pts[0], + ETA_WELD) + probes[k] = (sig, tau) + return child, stokes, probes + + +def stress_components(child, stokes, tag): + x, y = child.X + v = stokes.Unknowns.u + p = stokes.Unknowns.p + exprs = dict( + sxx=-p.sym[0] + 2 * common.ETA * v.sym[0].diff(x), + syy=-p.sym[0] + 2 * common.ETA * v.sym[1].diff(y), + sxy=common.ETA * (v.sym[0].diff(y) + v.sym[1].diff(x))) + out = {} + for name, expr in exprs.items(): + # P0 (cell) stress: projecting the rough near-fault stress onto + # CONTINUOUS P1 rings at the node scale (measured: residual rms + # 0.26 at half-wavelength h/2); cellwise averages are honest and + # pixel-scale at this resolution. + s_var = uw.discretisation.MeshVariable(f"{name}_{tag}", child, 1, + degree=0, continuous=False) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = expr + proj.smoothing = 0.0 + proj.solve() + row = common.split_mesh_cell_rows(child, s_var) + out[name] = np.asarray(s_var.data[:, 0])[row].copy() + return out + + +cache = os.path.join(D, "_california_probes.npz") +if os.path.exists(cache): + data = dict(np.load(cache, allow_pickle=True)) + print("loaded cached run") +else: + t_wall = time.perf_counter() + child, s1, probes1 = build_and_solve(trunk_free=True) + print(f"[timing] slipping solve + probes: " + f"{time.perf_counter() - t_wall:.1f} s") + t_wall = time.perf_counter() + s_saf, V_saf = common.slip_vs_position( + s1, _t, centre=_A + 0.47 * _t, name="SAF") + comp1 = stress_components(child, s1, "a") + + s0 = common.stokes_on(child, + common.boundary_simple_shear(child, TREND, + TAU0)) + s0.add_fault_bc(ETA_WELD, boundary="SAF", normal=saf_normal(child)) + for k in MINORS: + s0.add_fault_bc(ETA_WELD, boundary=k) + fault_contact.solve_with_fault(s0, picard=2) + comp0 = stress_components(child, s0, "b") + print(f"[timing] welded solve + all projections: " + f"{time.perf_counter() - t_wall:.1f} s") + probes0 = {} + for k, pts in MINORS.items(): + _s, _xy, sig, tau = common.probe_nodes(s0, k, pts[1] - pts[0], + ETA_WELD) + probes0[k] = (sig, tau) + + # Delta CFF on boundary-parallel planes (per cell) + beta = np.radians(TREND) + nx, ny = -np.sin(beta), np.cos(beta) + tx, ty = np.cos(beta), np.sin(beta) + + def resolve(c): + s_nn = (c["sxx"] * nx * nx + 2 * c["sxy"] * nx * ny + + c["syy"] * ny * ny) + s_t = (c["sxx"] * tx * nx + c["sxy"] * (tx * ny + ty * nx) + + c["syy"] * ty * ny) + return s_nn, s_t + + nn0, t0 = resolve(comp0) + nn1, t1 = resolve(comp1) + tau_dir = np.sign(np.median(t0)) + _pts, _faces = common.split_mesh_cell_render(child) + # cell centroids, for the far-field gauge anchor + fc = np.asarray(_faces).reshape(-1, 4)[:, 1:] + _cent = np.asarray(_pts)[fc].mean(axis=1) + data = dict(field_dcff=tau_dir * (t1 - t0) + MU_P * (nn1 - nn0), + field_points=_pts, field_faces=_faces, + field_centroids=_cent, saf_v=V_saf) + for k in MINORS: + data[f"{k}_sig0"], data[f"{k}_tau0"] = probes0[k] + data[f"{k}_sig1"], data[f"{k}_tau1"] = probes1[k] + np.savez(cache, **data) + data = dict(np.load(cache, allow_pickle=True)) + +v_med = float(np.median(data["saf_v"])) +# dextral: with the tangent pointing NW and the split's Plus side on +# its LEFT (SW, the Pacific side), a POSITIVE jump (v+ - v-) along +t +# means the Pacific side moves NW relative to North America. +sense = "right-lateral" if v_med > 0 else "LEFT-LATERAL?!" +print(f"SAF slip: median tangential jump {v_med:+.3f} ({sense})") + +# ---- the field render ------------------------------------------------------ +dcff_field, GAUGE_C = common.far_field_anchor( + data["field_centroids"], data["field_dcff"], + [SAF_PTS[:26], SAF_PTS[24:]] + list(MINORS.values()), cut=0.18) +print(f"far-field gauge constant removed: {GAUGE_C:+.4f}") +pvm = pv.PolyData(np.asarray(data["field_points"], dtype=float), + faces=np.asarray(data["field_faces"], dtype=np.int64)) +pvm.cell_data["dcff"] = dcff_field +pl = pv.Plotter(off_screen=True, window_size=(1000, 950)) +pl.set_background("white") +pl.add_mesh(pvm, scalars="dcff", cmap="RdBu_r", clim=(-1.0, 1.0), + show_edges=False, lighting=False, + scalar_bar_args=dict(title="dCFF", color="black")) + + +def polyline(pts): + return pv.lines_from_points( + np.column_stack([pts, np.full(len(pts), 0.001)])) + + +pl.add_mesh(polyline(SAF_PTS), color="black", line_width=5.0, + lighting=False) +for k, pts in MINORS.items(): + pl.add_mesh(polyline(pts), color=COLOUR[k], line_width=4.0, + lighting=False) + +# measured slip sense, half-arrows either side of the southern leg +mid = _A + 0.20 * _t +sense_sign = np.sign(v_med) +for pm in (+1.0, -1.0): + base = mid + pm * 0.05 * _n - pm * sense_sign * 0.06 * _t + pl.add_arrows(np.array([np.append(base, 0.002)]), + np.array([np.append(pm * sense_sign * _t, 0.0)]), + mag=0.11, color="black") + +pl.add_point_labels( + np.array([[0.30, 0.44, 0.002], [0.86, 0.12, 0.002], + [0.68, 0.50, 0.002], [0.60, 0.82, 0.002], + [0.84, 0.20, 0.002], [0.40, 0.24, 0.002]]), + ["San Andreas (N)", "San Andreas (S)", "Garlock", "ECSZ", + "San Jacinto", "Transverse Ranges\n(restraining bend)"], + font_size=20, text_color="black", shape=None, always_visible=True, + show_points=False) + +pl.view_xy() +pl.camera.parallel_projection = True +pl.camera.parallel_scale = 0.48 +pl.camera.focal_point = (0.48, 0.44, 0.0) +field_png = os.path.join(D, "_california_field.png") +pl.screenshot(field_png) +pl.close() + +# ---- the figure: field + one Mohr panel per fault group -------------------- +fig = plt.figure(figsize=(13.2, 6.6)) +gs = fig.add_gridspec(3, 3, width_ratios=[2.6, 1.15, 0.06]) + +axf = fig.add_subplot(gs[:, 0]) +axf.imshow(plt.imread(field_png)) +axf.set_xticks([]) +axf.set_yticks([]) +axf.set_title(r"$\Delta$CFF on boundary-parallel planes ($\mu' = 0.4$);" + "\nthe San Andreas slips (right-lateral) through its " + "restraining bend, neighbours welded as probes\n" + "(schematic geometry, not to scale)", + fontsize=9) + +for row, (label, members, col) in enumerate(GROUPS): + ax = fig.add_subplot(gs[row, 1]) + ss = np.linspace(-0.4, P0 + 1.9, 80) + strength = np.maximum(COH + MU_P * ss, 0.0) + for sgn in (+1, -1): + ax.plot(ss, sgn * strength, "-", color="0.4", lw=0.9) + ax.fill_between(ss, strength, 2.6, color="#c62828", alpha=0.06, lw=0) + ax.fill_between(ss, -strength, -2.6, color="#c62828", alpha=0.06, + lw=0) + tt = np.linspace(0, 2 * np.pi, 150) + ax.plot(P0 + TAU0 * np.cos(tt), TAU0 * np.sin(tt), "-", color="0.88", + lw=0.7) + medians = [] + for k in members: + t_hat = MINORS[k][1] - MINORS[k][0] + c0 = float(np.median(data[f"{k}_sig0"]) + - common.ambient_sigma_n_simple(TREND, t_hat, TAU0)) + sig0 = data[f"{k}_sig0"] - c0 + tau0 = data[f"{k}_tau0"] + sig1 = data[f"{k}_sig1"] - c0 - GAUGE_C / MU_P + tau1 = data[f"{k}_tau1"] + tau_dir = np.sign(np.median(tau0)) + dcff = tau_dir * (tau1 - tau0) + MU_P * (sig1 - sig0) + medians.append(np.median(dcff)) + sc0, sc1 = P0 - sig0, P0 - sig1 + ax.scatter(sc0, tau0, s=9, facecolors="none", edgecolors="0.6", + linewidths=0.7) + for j in range(0, len(sc0), 3): + ax.annotate("", xytext=(sc0[j], tau0[j]), + xy=(sc1[j], tau1[j]), + arrowprops=dict(arrowstyle="->", lw=0.45, + color="0.6")) + pts = ax.scatter(sc1, tau1, c=dcff, cmap="RdBu_r", s=16, + vmin=-0.3, vmax=0.3, zorder=5, + edgecolors="0.3", linewidths=0.2) + ax.axhline(0, color="0.92", lw=0.5) + ax.set_aspect("equal") + ax.set_xlim(-0.4, P0 + 1.9) + ax.set_ylim(-1.7, 1.7) + ax.tick_params(labelsize=7) + for spine in ax.spines.values(): + spine.set_color(col) + spine.set_linewidth(1.8) + ax.set_title(f"{label}: median $\\Delta$CFF " + f"{np.median(medians):+.2f}", fontsize=8.5, color=col) + +cax = fig.add_subplot(gs[:, 2]) +fig.colorbar(pts, cax=cax, label=r"node $\Delta$CFF") +fig.suptitle("A San Andreas slip event (Big Bend), read by its " + "neighbours", fontsize=11.5) +fig.tight_layout() +out = os.path.join(D, "california.png") +fig.savefig(out, dpi=200) +print("wrote", out) diff --git a/docs/advanced/figures/fault-examples/common.py b/docs/advanced/figures/fault-examples/common.py new file mode 100644 index 000000000..d55abe74d --- /dev/null +++ b/docs/advanced/figures/fault-examples/common.py @@ -0,0 +1,326 @@ +"""Shared harness for the split-node fault teaching examples. + +One base mesh, re-faulted per case (the non-cumulative add_fault +pattern); P2 velocity / P0 discontinuous pressure per the fault +pressure-space ruling; slip and traction read through the DOF pairing. +Run inside the fault-split-node worktree env with its bin on PATH. +""" +import numpy as np + +import underworld3 as uw +from underworld3.utilities import fault_contact + +ETA = 1.0 +CENTRE = np.array([0.5, 0.5]) + + +def base_mesh(cell_size=0.04): + return uw.meshing.UnstructuredSimplexBox(cellSize=cell_size) + + +def fault_segment(angle_deg, half_length=0.2, centre=CENTRE): + t = np.array([np.cos(np.radians(angle_deg)), + np.sin(np.radians(angle_deg))]) + return np.array([centre - half_length * t, centre + half_length * t]) + + +def split_with_fault(mesh, points, name="Fault"): + return mesh.add_fault((name, points)) + + +def stokes_on(child, drive, name="Fault"): + """Stokes with the wall drive on all four walls; fault BC added by + the caller. Variable names are derived from the mesh instance so + repeated calls on fresh children never collide.""" + tag = f"f{stokes_on.counter}" + stokes_on.counter += 1 + v = uw.discretisation.MeshVariable(f"V_{tag}", child, child.dim, + degree=2) + p = uw.discretisation.MeshVariable(f"P_{tag}", child, 1, degree=0, + continuous=False) + stokes = uw.systems.Stokes(child, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = ETA + stokes.bodyforce = [0.0, 0.0] + for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc(drive, wall) + # all-Dirichlet walls leave the pressure level to the solver unless + # the constant nullspace is declared — without it, two solves with + # different fault laws land on DIFFERENT gauges and any differenced + # quantity (Delta CFF above all) inherits a spurious constant + stokes.petsc_use_pressure_nullspace = True + stokes.tolerance = 1e-6 + return stokes + + +stokes_on.counter = 0 + + +def simple_shear(child, rate=1.0): + """v = (rate (y - 1/2), 0): sigma_xy = ETA * rate everywhere.""" + x, y = child.X + return (rate * (y - 0.5), 0.0) + + +def shear_plus_stretch(child, a=0.5, gamma=1.0): + """v = (a(x-c) + gamma(y-c), -a(y-c)): deviatoric stress + [[2 eta a, eta gamma], [eta gamma, -2 eta a]], Mohr radius + eta sqrt(4 a^2 + gamma^2).""" + x, y = child.X + return (a * (x - 0.5) + gamma * (y - 0.5), -a * (y - 0.5)) + + +def slip_profile(stokes, name="Fault"): + return fault_contact.fault_slip(stokes, name, + stokes._rotated_freeslip_info) + + +def normal_traction(stokes, name="Fault"): + return fault_contact.fault_normal_traction( + stokes, name, stokes._rotated_freeslip_info) + + +def inner(s, fraction=0.6): + """Mask selecting the central `fraction` of the along-fault range — + tip fields carry the crack singularity, the middle is the gauge.""" + lo = s.min() + (1 - fraction) / 2 * (s.max() - s.min()) + hi = s.max() - (1 - fraction) / 2 * (s.max() - s.min()) + return (s >= lo) & (s <= hi) + + +def slip_vs_position(stokes, tangent, centre=CENTRE, name="Fault"): + """(s, V) read through the pairing with EXACT positions: s is the + signed along-fault coordinate of each pair about the fault centre, + V the signed tangential jump. `fault_slip`'s arc-length origin sits + at the first PAIR, not the tip, so profiles built from it are offset + by half a node — this is the plotting-grade version.""" + coords, jumps, _normals = fault_contact.fault_pair_jumps( + stokes, name, stokes._rotated_freeslip_info) + t = np.asarray(tangent, dtype=float) + t = t / np.linalg.norm(t) + s = (coords - centre) @ t + V = jumps @ t + order = np.argsort(s) + return s[order], V[order] + + +def pure_shear_drive(child, phi_deg, tau0=1.0): + """Uniform pure shear with the COMPRESSION axis at ``phi_deg`` to x: + sigma' = tau0 (e_perp e_perp - e_phi e_phi), i.e. sigma'_xx = + -tau0 cos 2phi and sigma'_xy = -tau0 sin 2phi. (Check: compression + at phi = 45 deg gives sigma_xy = -tau0; a plane of strike + phi - 45 deg carries the full resolved shear.) Irrotational + velocity field, imposed as Dirichlet on all walls — rotating phi + rotates the whole regional stress field.""" + x, y = child.X + two_phi = 2.0 * np.radians(phi_deg) + exx = -tau0 * np.cos(two_phi) / (2.0 * ETA) + exy = -tau0 * np.sin(two_phi) / (2.0 * ETA) + return (exx * (x - 0.5) + exy * (y - 0.5), + exy * (x - 0.5) - exx * (y - 0.5)) + + +def boundary_simple_shear(child, trend_deg, rate=1.0): + """RIGHT-LATERAL simple shear parallel to a plate-boundary trend: + v = rate ((X - c) . n) t with t the boundary-parallel direction and + n its CCW normal — material on the +n side moves along +t, which an + observer on the fault sees moving to the RIGHT (dextral). Resolved + shear tau0 = ETA * rate on boundary-parallel planes.""" + x, y = child.X + t = np.array([np.cos(np.radians(trend_deg)), + np.sin(np.radians(trend_deg))]) + n = np.array([-t[1], t[0]]) + s = (x - 0.5) * n[0] + (y - 0.5) * n[1] + return (rate * s * t[0], rate * s * t[1]) + + +def ambient_sigma_n_simple(trend_deg, tangent, tau0=1.0): + """Analytic ambient normal stress (tension-positive) on a plane of + direction ``tangent`` under boundary_simple_shear(trend): + sigma_nn = -tau0 sin 2(theta - trend).""" + theta = np.degrees(np.arctan2(tangent[1], tangent[0])) + return -tau0 * np.sin(2.0 * np.radians(theta - trend_deg)) + + +def probe_nodes(stokes, name, tangent, eta_weld): + """Per-node stress probes on a WELDED fault: along-fault coordinate, + position, signed normal traction (tension-positive, as measured) and + signed shear traction from the weld's own law tau = eta_f V. + + Works with SEVERAL law-carrying faults on one mesh: the interface + assembler holds every fault's nodes, so this fault's are selected + through its OWN pairing (the plus points), never by position.""" + from underworld3.utilities.rotated_bc import _point_coord + + assembler = fault_contact._InterfaceAssembler(stokes, include=(name,)) + sig_all = assembler.nodal_normal_traction( + stokes, stokes._rotated_freeslip_info["reaction"]) + plus = set(stokes.mesh._fault_point_pairs[name].values()) + + dm = stokes.dm + dim = stokes.mesh.dim + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + rows = [(np.asarray(_point_coord(dm, dim, cvec, csec, v0, v1, q)), + sig_all[k]) for q, k in assembler._points.items() + if q in plus] + xy_sig = np.array([r[0] for r in rows]) + sig = np.array([r[1] for r in rows]) + + t = np.asarray(tangent, dtype=float) + t = t / np.linalg.norm(t) + order_sig = np.argsort(xy_sig @ t) + xy_sig, sig = xy_sig[order_sig], sig[order_sig] + + coords, jumps, _normals = fault_contact.fault_pair_jumps( + stokes, name, stokes._rotated_freeslip_info) + s = coords @ t + order = np.argsort(s) + V = (jumps @ t)[order] + assert len(sig) == len(V), "pair-node sets disagree" + assert np.allclose(coords[order], xy_sig), "node ordering disagrees" + return (s[order], coords[order], sig, eta_weld * V) + + +def mohr_probe(theta, a_rate=0.5, gamma=1.0, eta_weld=None, + half_length=0.2, cell_size=0.04): + """One welded-fault stress probe: (sigma_n, tau_signed) at fault + angle `theta` under the shear_plus_stretch drive. The weld's own + law reads the shear traction (tau = eta_f V); the no-opening + reaction reads the normal traction.""" + if eta_weld is None: + eta_weld = 200.0 * ETA / half_length + child = split_with_fault(base_mesh(cell_size), + fault_segment(theta, half_length)) + stokes = stokes_on(child, shear_plus_stretch(child, a_rate, gamma)) + stokes.add_fault_bc(eta_weld, boundary="Fault") + fault_contact.solve_with_fault(stokes, picard=2) + s, V, _leak = slip_profile(stokes) + s_n, sig = normal_traction(stokes) + tau = eta_weld * float(np.median(V[inner(s)])) + sigma_n = float(np.median(sig[inner(s_n)])) + return sigma_n, tau + + +def ambient_sigma_n(phi_deg, tangent, tau0=1.0): + """The analytic ambient normal stress (tension-positive) on a plane + of tangent direction ``tangent`` under pure_shear_drive(phi): + sigma_nn = tau0 cos 2(phi - theta). Used to anchor a probe set's + absolute pressure gauge, which a closed-box solve does not fix.""" + theta = np.degrees(np.arctan2(tangent[1], tangent[0])) + return tau0 * np.cos(2.0 * np.radians(phi_deg - theta)) + + +def far_field_anchor(points, dcff, segments, cut=0.3): + """Gauge the DIFFERENCED stress to the physics: a slip event changes + nothing far from the faults, so the far-field median of Delta CFF is + the spurious pressure-gauge constant between the two solves. Returns + (anchored dcff, the constant) — apply c / mu' to the probes' sigma.""" + p = np.asarray(points, dtype=float)[:, :2] + far = np.ones(len(p), dtype=bool) + for seg in segments: + a, b = np.asarray(seg[0]), np.asarray(seg[-1]) + t = b - a + L = np.linalg.norm(t) + t = t / L + s = np.clip((p - a) @ t, 0.0, L) + far &= np.linalg.norm(p - (a + s[:, None] * t), axis=1) > cut + c = float(np.median(np.asarray(dcff)[far])) + return np.asarray(dcff) - c, c + + +def signed_log(x, linthresh=0.02): + """Symmetric log transform for signed stress-change fields: linear + inside |x| < linthresh, logarithmic beyond — the near-fault values + saturate any linear colour scale and hide the far-field lobes.""" + x = np.asarray(x, dtype=float) + return np.sign(x) * np.log10(1.0 + np.abs(x) / linthresh) + + +def signed_log_annotations(values, linthresh=0.02): + """Scalar-bar tick positions/labels in transformed units.""" + return {float(np.sign(v) * np.log10(1.0 + abs(v) / linthresh)): + f"{v:+.2f}" if v else "0" for v in values} + + +def dedupe_average(points, values, decimals=9): + """Merge geometrically coincident points, averaging their values — + for RENDERING split-mesh fields only. The split mesh carries + coincident duplicated nodes along each fault; a Delaunay of the raw + cloud picks sides arbitrarily and paints alternating plus/minus + values along the trace (beading). Averaging collapses the jump in a + one-node strip that the drawn fault line covers.""" + pts = np.asarray(points, dtype=float) + vals = np.asarray(values, dtype=float) + key = np.round(pts[:, :2], decimals) + _uniq, inverse, counts = np.unique(key, axis=0, return_inverse=True, + return_counts=True) + sums = np.zeros(len(counts)) + np.add.at(sums, inverse, vals) + first = np.full(len(counts), -1, dtype=int) + for i, g in enumerate(inverse): + if first[g] < 0: + first[g] = i + return pts[first], sums / counts + + +def split_mesh_cell_render(child): + """Points (mesh vertices) + faces (true cell connectivity) for + rendering CELL data on a split mesh. Node-order mapping is trivial + (v - vS) because the points ARE the vertices; pair cell values via + split_mesh_cell_rows.""" + dm = child.dm + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + faces = [] + for c in range(cS, cE): + closure, _ = dm.getTransitiveClosure(c) + tri = [int(p) - vS for p in closure if vS <= int(p) < vE] + faces.extend([3, *tri]) + return (np.column_stack([X, np.zeros(len(X))]), + np.asarray(faces, dtype=np.int64)) + + +def split_mesh_cell_rows(child, var): + """Data-row index of a P0 variable for each cell, in plex cell + order (rank of the combined-section offsets).""" + dm = child.dm + sec = dm.getLocalSection() + cS, cE = dm.getHeightStratum(0) + offs = np.array([sec.getFieldOffset(c, var.field_id) + for c in range(cS, cE)]) + row = np.empty(len(offs), dtype=np.int64) + row[np.argsort(offs)] = np.arange(len(offs)) + return row + + +def split_mesh_faces(child, var): + """The mesh's own triangulation as a pyvista faces array, indexed in + ``var``'s DOF order (P1 vertex variable). Rendering a split-mesh + field through delaunay_2d of the DOF cloud is WRONG twice over: the + coincident fault pairs get arbitrary side-picking, and unconstrained + Delaunay edges hop across a curved trace — both paint beading along + the fault. The true connectivity keeps the two sheets separate and + renders the jump sharply.""" + dm = child.dm + sec = dm.getLocalSection() + vS, vE = dm.getDepthStratum(0) + cS, cE = dm.getHeightStratum(0) + fid = var.field_id + # getFieldOffset indexes the COMBINED local vector; the variable's + # own rows are the RANKS of those offsets (dof order = offset order) + offs = np.array([sec.getFieldOffset(v, fid) for v in range(vS, vE)]) + row = np.empty(len(offs), dtype=np.int64) + row[np.argsort(offs)] = np.arange(len(offs)) + X = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, 2) + assert np.allclose(np.asarray(var.coords)[row], X), \ + "variable DOF order does not match the offset ranking" + faces = [] + for c in range(cS, cE): + closure, _ = dm.getTransitiveClosure(c) + tri = [int(row[int(p) - vS]) for p in closure + if vS <= int(p) < vE] + faces.extend([3, *tri]) + return np.asarray(faces, dtype=np.int64) diff --git a/docs/advanced/figures/fault-examples/interacting-faults.png b/docs/advanced/figures/fault-examples/interacting-faults.png new file mode 100644 index 000000000..0720487f7 Binary files /dev/null and b/docs/advanced/figures/fault-examples/interacting-faults.png differ diff --git a/docs/advanced/figures/fault-examples/interacting-rotation.png b/docs/advanced/figures/fault-examples/interacting-rotation.png new file mode 100644 index 000000000..0498ad47b Binary files /dev/null and b/docs/advanced/figures/fault-examples/interacting-rotation.png differ diff --git a/docs/advanced/figures/fault-examples/interacting_faults.py b/docs/advanced/figures/fault-examples/interacting_faults.py new file mode 100644 index 000000000..7e4bc647e --- /dev/null +++ b/docs/advanced/figures/fault-examples/interacting_faults.py @@ -0,0 +1,251 @@ +"""Interacting faults, King-style: stress transfer read on the Mohr plane. + +Two en echelon faults (right-stepping, overlapping). The SOURCE slips +freely; the RECEIVER is welded — a passive per-node stress probe, the +instrument built in the Mohr examples. Two solves on the SAME mesh: + + (0) both faults welded -> the ambient regional state + (1) source slips, receiver weld -> the perturbed state + +The difference is the classical Coulomb stress transfer: + +- the FIELD: Delta CFF = tau_dir * d sigma_xy + mu' * d sigma_yy on + receiver-parallel planes (King's mu' = 0.4), the familiar lobes; +- ON the receiver: each node's probe MOVES in the (sigma, tau) plane, + and the displacement decomposes the CFF change into its shear and + normal(unclamping) parts — Coulomb stress transfer in the diagram + students already know how to read. + +Rotating the regional compression axis phi (the boundary velocities +rotate with it) changes the interaction pattern; the phi sweep shows +the receiver cloud pushed toward failure under one orientation and +away under another. +""" +import os +import time + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv + +import underworld3 as uw +import underworld3.visualisation as vis +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +MU_P = 0.4 # King's effective friction +TAU0 = 1.0 +ETA_WELD = 200.0 * common.ETA / 0.18 +PHIS = (20.0, 45.0, 70.0) # regional compression axis + +SOURCE = np.array([[0.22, 0.46], [0.58, 0.46]]) +# The receiver sits in the source's ALONG-STRIKE tip lobe, where the +# transfer is strongly positive at its near end and fades along it. +RECEIVER = np.array([[0.65, 0.44], [0.86, 0.44]]) +T_HAT = np.array([1.0, 0.0]) + +# The friction bookkeeping is dressed with a declared confining +# pressure and cohesion — NEITHER changes Delta CFF (both are +# constants under differencing), they only place the failure envelope +# where the teaching needs it: P0 shifts the whole circle into +# compression (no tensile regime anywhere), and the cohesive envelope +# tau = +-(C + mu' sigma) arcs over it with the phi = 45 receiver +# ambient sitting ~0.15 below the line, so the loaded nodes CROSS. +P0 = 1.0 +COH = 0.75 + + +def solve_pair(phi, source_free): + child = common.base_mesh(0.016).add_fault( + [("Source", SOURCE), ("Receiver", RECEIVER)]) + stokes = common.stokes_on(child, common.pure_shear_drive(child, phi, + TAU0)) + stokes.add_fault_bc(0 if source_free else ETA_WELD, boundary="Source") + stokes.add_fault_bc(ETA_WELD, boundary="Receiver") + fault_contact.solve_with_fault(stokes, picard=2) + s, xy, sig, tau = common.probe_nodes(stokes, "Receiver", T_HAT, + ETA_WELD) + return child, stokes, (s, xy, sig, tau) + + +def stress_fields(child, stokes): + """sigma_xy and sigma_yy (full stress, pressure included) projected + to continuous P1 — the components CFF on horizontal planes needs.""" + x, y = child.X + v = stokes.Unknowns.u + p = stokes.Unknowns.p + out = [] + for name, expr in (("sxy", common.ETA * (v.sym[0].diff(y) + + v.sym[1].diff(x))), + ("syy", -p.sym[0] + 2 * common.ETA + * v.sym[1].diff(y))): + # P0 cell stress: continuous-P1 projection rings at node scale + # near the fault (see california.py) + s_var = uw.discretisation.MeshVariable( + f"{name}_{stress_fields.counter}", child, 1, degree=0, + continuous=False) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = expr + proj.smoothing = 0.0 + proj.solve() + row = common.split_mesh_cell_rows(child, s_var) + out.append(np.asarray(s_var.data[:, 0])[row].copy()) + stress_fields.counter += 1 + return out, s_var # arrays share the P1 layout + + +stress_fields.counter = 0 + +cache = os.path.join(D, "_interacting_probes.npz") +if os.path.exists(cache): + data = dict(np.load(cache, allow_pickle=True)) + print("loaded cached sweep") +else: + data = {} + for phi in PHIS: + t_wall = time.perf_counter() + _c0, s0, probe0 = solve_pair(phi, source_free=False) + child, s1, probe1 = solve_pair(phi, source_free=True) + data[f"probe0_{phi}"] = np.array(probe0[0]), probe0[1], probe0[2], \ + probe0[3] + data[f"p0_s_{phi}"], data[f"p0_xy_{phi}"] = probe0[0], probe0[1] + data[f"p0_sig_{phi}"], data[f"p0_tau_{phi}"] = probe0[2], probe0[3] + data[f"p1_sig_{phi}"], data[f"p1_tau_{phi}"] = probe1[2], probe1[3] + if phi == 45.0: + # the Delta CFF FIELD needs both solves' stress components on + # one mesh: redo the welded reference ON the slipping child + (sxy1, syy1), s_var = stress_fields(child, s1) + s0b = common.stokes_on(child, common.pure_shear_drive( + child, phi, TAU0)) + s0b.add_fault_bc(ETA_WELD, boundary="Source") + s0b.add_fault_bc(ETA_WELD, boundary="Receiver") + fault_contact.solve_with_fault(s0b, picard=2) + (sxy0, syy0), _ = stress_fields(child, s0b) + tau_dir = np.sign(np.median(data[f"p0_tau_{phi}"])) + dcff = tau_dir * (sxy1 - sxy0) + MU_P * (syy1 - syy0) + _pts, _faces = common.split_mesh_cell_render(child) + data["field_points"] = _pts + data["field_faces"] = _faces + data["field_dcff"] = dcff + print(f"phi {phi}: receiver ambient tau " + f"{np.median(data[f'p0_tau_{phi}']):+.3f}, sigma_n " + f"{np.median(data[f'p0_sig_{phi}']):+.3f} " + f"[{time.perf_counter() - t_wall:.1f} s]") + np.savez(cache, **{k: v for k, v in data.items() + if not k.startswith("probe0")}) + data = dict(np.load(cache, allow_pickle=True)) + +# ---- figure A: the field + the receiver's Mohr move (phi = 45) ----------- +# gauge the difference to the far field (see common.far_field_anchor) +_fc = np.asarray(data["field_faces"]).reshape(-1, 4)[:, 1:] +_cent = np.asarray(data["field_points"])[_fc].mean(axis=1) +dcff_field, GAUGE_C = common.far_field_anchor( + _cent, data["field_dcff"], (SOURCE, RECEIVER)) +print(f"far-field gauge constant removed: {GAUGE_C:+.4f}") +pvm = pv.PolyData(np.asarray(data["field_points"], dtype=float), + faces=np.asarray(data["field_faces"], dtype=np.int64)) +pvm.cell_data["dcff"] = dcff_field +pl = pv.Plotter(off_screen=True, window_size=(1050, 820)) +pl.set_background("white") +# linear colour scale with generous limits (the gate-study +# convention): pale far field, graded lobes +pl.add_mesh(pvm, scalars="dcff", cmap="RdBu_r", clim=(-1.0, 1.0), + show_edges=False, lighting=False, + scalar_bar_args=dict(title="dCFF", color="black")) +for pts, col, w in ((SOURCE, "black", 4.0), (RECEIVER, "#1a6b1a", 4.0)): + line = pv.Line(tuple(pts[0]) + (0.001,), tuple(pts[1]) + (0.001,)) + pl.add_mesh(line, color=col, line_width=w, lighting=False) +pl.view_xy() +pl.camera.parallel_projection = True +pl.camera.parallel_scale = 0.33 +pl.camera.focal_point = (0.5, 0.5, 0.0) +field_png = os.path.join(D, "_interacting_field.png") +pl.screenshot(field_png) +pl.close() + + +def mohr_panel(ax, phi, legend=False): + # anchor the welded probes to the analytic ambient for this drive, + # the after-probes to the far-field-anchored difference on top + c0 = float(np.median(data[f"p0_sig_{phi}"]) + - common.ambient_sigma_n(phi, T_HAT, TAU0)) + sig0 = data[f"p0_sig_{phi}"] - c0 + tau0 = data[f"p0_tau_{phi}"] + sig1 = data[f"p1_sig_{phi}"] - c0 - GAUGE_C / MU_P + tau1 = data[f"p1_tau_{phi}"] + tau_dir = np.sign(np.median(tau0)) + dcff = tau_dir * (tau1 - tau0) + MU_P * (sig1 - sig0) + # geo convention, with the declared confining pressure superposed + sc0, sc1 = P0 - sig0, P0 - sig1 + ss = np.linspace(-0.4, P0 + 1.8, 90) + strength = np.maximum(COH + MU_P * ss, 0.0) + for sgn in (+1, -1): + ax.plot(ss, sgn * strength, "-", color="0.35", lw=1.1, + label=(r"envelope $\tau = \pm(C + \mu'\sigma)$" + if sgn > 0 and legend else None)) + ax.fill_between(ss, strength, 2.4, color="#c62828", alpha=0.06, + lw=0) + ax.fill_between(ss, -strength, -2.4, color="#c62828", alpha=0.06, + lw=0) + tt = np.linspace(0, 2 * np.pi, 200) + ax.plot(P0 + TAU0 * np.cos(tt), TAU0 * np.sin(tt), "-", + color="0.85", lw=0.8, + label="regional circle" if legend else None) + ax.scatter(sc0, tau0, s=16, facecolors="none", edgecolors="0.55", + linewidths=1.0, + label="receiver nodes, before" if legend else None) + for k in range(0, len(sc0), 2): + ax.annotate("", xytext=(sc0[k], tau0[k]), xy=(sc1[k], tau1[k]), + arrowprops=dict(arrowstyle="->", lw=0.6, + color="0.65")) + pts = ax.scatter(sc1, tau1, c=dcff, cmap="RdBu_r", s=26, + vmin=-0.4, vmax=0.4, zorder=5, edgecolors="0.3", + linewidths=0.3, + label="after (colour: dCFF)" if legend else None) + ax.axhline(0, color="0.9", lw=0.5) + ax.set_xlim(-0.4, P0 + 1.8) + ax.set_ylim(-1.75, 1.75) + ax.set_aspect("equal") + ax.set_title(rf"$\phi = {phi:.0f}°$ " + rf"(median dCFF {np.median(dcff):+.2f})", fontsize=9) + ax.set_xlabel(r"$\sigma$ (compression +, $P_0 = 1$)", fontsize=8) + return pts + + +fig = plt.figure(figsize=(11.5, 4.9)) +axf = fig.add_subplot(1, 2, 1) +axf.imshow(plt.imread(field_png)) +axf.set_xticks([]) +axf.set_yticks([]) +axf.set_title(r"$\Delta$CFF field ($\mu' = 0.4$), source slips freely;" + "\nreceiver (green) welded as a probe", fontsize=9) +axm = fig.add_subplot(1, 2, 2) +pts = mohr_panel(axm, 45.0, legend=True) +axm.set_ylabel(r"$\tau$", fontsize=9) +axm.legend(fontsize=7, loc="lower right") +fig.colorbar(pts, ax=axm, label=r"node $\Delta$CFF", shrink=0.8) +fig.suptitle("En echelon interaction: the receiver's probes move in the " + "Mohr plane", fontsize=11) +fig.tight_layout() +out = os.path.join(D, "interacting-faults.png") +fig.savefig(out, dpi=200) +print("wrote", out) + +# ---- figure B: rotate the regional stress -------------------------------- +fig, axes = plt.subplots(1, len(PHIS), figsize=(11.5, 4.2), + sharey=True) +for ax, phi in zip(axes, PHIS): + pts = mohr_panel(ax, phi, legend=(phi == PHIS[0])) +axes[0].set_ylabel(r"$\tau$", fontsize=9) +axes[0].legend(fontsize=6.5, loc="lower right") +fig.colorbar(pts, ax=axes, label=r"node $\Delta$CFF", shrink=0.8) +fig.suptitle("Rotating the regional compression axis changes the " + "interaction", fontsize=11) +out = os.path.join(D, "interacting-rotation.png") +fig.savefig(out, dpi=200) +print("wrote", out) diff --git a/docs/advanced/figures/fault-examples/ladder.png b/docs/advanced/figures/fault-examples/ladder.png new file mode 100644 index 000000000..1b9292773 Binary files /dev/null and b/docs/advanced/figures/fault-examples/ladder.png differ diff --git a/docs/advanced/figures/fault-examples/ladder.py b/docs/advanced/figures/fault-examples/ladder.py new file mode 100644 index 000000000..8dd1d67b9 --- /dev/null +++ b/docs/advanced/figures/fault-examples/ladder.py @@ -0,0 +1,130 @@ +"""The fault-strength ladder: one fault, one drive, four laws. + +A horizontal fault under far-field simple shear (resolved shear stress +tau_infty = eta * rate = 1 on the fault plane), solved with each rung of +the constitutive ladder. Left: the slip profiles V(s), with the exact +zero-slip tips appended (positions read through the DOF pairing — the +generic fault_slip arc-length starts at the first pair, not the tip). +Right: the shear-stress field sigma_xy for each rung, stacked for +context — the stress drop shadows the fault where it slips, and the +welded/stuck rungs leave the far-field stress untouched. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import numpy as np +import pyvista as pv + +import underworld3 as uw +import underworld3.visualisation as vis +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +SIGMA_N = 2.0 # prescribed normal stress for friction + +segment = common.fault_segment(0.0, half_length=HALF) +runs = [] + + +def render_shear_stress(stokes, child, tag): + """sigma_xy projected to continuous P1 (smoothing 0), rendered as a + cropped planform panel with the fault trace overlaid.""" + x, y = child.X + v = stokes.Unknowns.u + s_var = uw.discretisation.MeshVariable(f"Sxy_{tag}", child, 1, degree=1) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = common.ETA * (v.sym[0].diff(y) + v.sym[1].diff(x)) + proj.smoothing = 0.0 + proj.solve() + + pv_m = vis.meshVariable_to_pv_mesh_object(s_var) + pv_m.point_data["S"] = np.asarray(s_var.data[:, 0]) + pl = pv.Plotter(off_screen=True, window_size=(560, 340)) + pl.set_background("white") + pl.add_mesh(pv_m, scalars="S", cmap="RdBu_r", clim=(0.0, 2.0), + show_edges=False, lighting=False, show_scalar_bar=False) + trace = pv.Line((0.5 - HALF, 0.5, 0.001), (0.5 + HALF, 0.5, 0.001)) + pl.add_mesh(trace, color="black", line_width=2.5, lighting=False) + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.parallel_scale = 0.26 + pl.camera.focal_point = (0.5, 0.5, 0.0) + out = os.path.join(D, f"_ladder_stress_{tag}.png") + pl.screenshot(out) + pl.close() + return out + + +def solve_case(label, register, tag): + child = common.split_with_fault(common.base_mesh(0.04), segment) + stokes = common.stokes_on(child, common.simple_shear(child)) + register(stokes) + fault_contact.solve_with_fault(stokes, picard=2) + s, V = common.slip_vs_position(stokes, tangent=(1.0, 0.0)) + # the unsplit tips carry exactly zero slip by construction + s = np.concatenate([[-HALF], s, [HALF]]) + V = np.concatenate([[0.0], V, [0.0]]) + panel = render_shear_stress(stokes, child, tag) + runs.append((label, s, np.abs(V), panel)) + print(f"{label:34s} peak slip {np.abs(V).max():.4f}") + + +solve_case("frictionless", + lambda st: st.add_fault_bc(0, boundary="Fault"), "free") +solve_case(r"viscous, $\eta_f = \eta/a$", + lambda st: st.add_fault_bc(common.ETA / HALF, boundary="Fault"), + "visc") +solve_case(r"Coulomb, $\mu\sigma_n = 0.6 < \tau_\infty$", + lambda st: fault_contact.add_coulomb_fault_bc( + st, 0.3, "Fault", sigma_n=SIGMA_N, V0=1e-4), "weak") +solve_case(r"rate-state (steady, $f_{ss}\sigma_n \approx 0.86$)", + lambda st: fault_contact.add_rate_state_fault_bc( + st, 0.42, "Fault", a=0.02, b=0.01, V0=1e-3, Dc=1e-2, + sigma_n=SIGMA_N), "rs") +solve_case(r"Coulomb, $\mu\sigma_n = 1.2 > \tau_\infty$ (stuck)", + lambda st: fault_contact.add_coulomb_fault_bc( + st, 0.6, "Fault", sigma_n=SIGMA_N, V0=1e-4), "stuck") + +colors = ["#c62828", "#e57373", "#d9960a", "#4a7bf7", "#555555"] + +fig = plt.figure(figsize=(9.8, 6.2)) +gs = gridspec.GridSpec(len(runs), 2, width_ratios=[1.7, 1.0], + wspace=0.12, hspace=0.15) + +ax = fig.add_subplot(gs[:, 0]) +for (label, s, V, _panel), col in zip(runs, colors): + ax.plot(s, V, ".-", ms=3.5, lw=1.0, color=col, label=label) +ss = np.linspace(-HALF, HALF, 200) +ax.plot(ss, np.sqrt(np.maximum(HALF**2 - ss**2, 0)) / HALF + * runs[0][2].max(), "k--", lw=0.8, + label="elliptical profile (shape)") +ax.set_xlabel("position along the fault $s$") +ax.set_ylabel("slip rate $|V(s)|$") +ax.set_title("The fault-strength ladder: one fault, four laws " + r"($\tau_\infty = 1$)") +ax.legend(fontsize=8, loc="upper right") +ax.set_xlim(-HALF * 1.15, HALF * 1.15) + +for k, ((label, _s, _V, panel), col) in enumerate(zip(runs, colors)): + axi = fig.add_subplot(gs[k, 1]) + axi.imshow(plt.imread(panel)) + axi.set_xticks([]) + axi.set_yticks([]) + for spine in axi.spines.values(): + spine.set_color(col) + spine.set_linewidth(2.2) + if k == 0: + axi.set_title(r"shear stress $\sigma_{xy}$ (0 → 2, " + r"white $= \tau_\infty$)", fontsize=8) + +fig.tight_layout() +out = os.path.join(D, "ladder.png") +fig.savefig(out, dpi=200) +print("wrote", out) diff --git a/docs/advanced/figures/fault-examples/mohr-circle-build.gif b/docs/advanced/figures/fault-examples/mohr-circle-build.gif new file mode 100644 index 000000000..1ae195223 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-circle-build.gif differ diff --git a/docs/advanced/figures/fault-examples/mohr-circle.png b/docs/advanced/figures/fault-examples/mohr-circle.png new file mode 100644 index 000000000..da35e6ea9 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-circle.png differ diff --git a/docs/advanced/figures/fault-examples/mohr-cohesion-build.gif b/docs/advanced/figures/fault-examples/mohr-cohesion-build.gif new file mode 100644 index 000000000..7f65e260a Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-cohesion-build.gif differ diff --git a/docs/advanced/figures/fault-examples/mohr-cohesion.png b/docs/advanced/figures/fault-examples/mohr-cohesion.png new file mode 100644 index 000000000..46354fe38 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-cohesion.png differ diff --git a/docs/advanced/figures/fault-examples/mohr-friction-build.gif b/docs/advanced/figures/fault-examples/mohr-friction-build.gif new file mode 100644 index 000000000..c1b8e3a22 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-friction-build.gif differ diff --git a/docs/advanced/figures/fault-examples/mohr-friction.png b/docs/advanced/figures/fault-examples/mohr-friction.png new file mode 100644 index 000000000..1dd769024 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-friction.png differ diff --git a/docs/advanced/figures/fault-examples/mohr-graded.png b/docs/advanced/figures/fault-examples/mohr-graded.png new file mode 100644 index 000000000..e6edcafc8 Binary files /dev/null and b/docs/advanced/figures/fault-examples/mohr-graded.png differ diff --git a/docs/advanced/figures/fault-examples/mohr_animate.py b/docs/advanced/figures/fault-examples/mohr_animate.py new file mode 100644 index 000000000..ba47c43de --- /dev/null +++ b/docs/advanced/figures/fault-examples/mohr_animate.py @@ -0,0 +1,129 @@ +"""The Mohr circle, built frame by frame as the fault rotates. + +The teaching version of mohr_circle.py: the left panel shows the fault +physically rotating in the box; the right panel shows its stress probe +sweeping around the Mohr circle at TWICE the rate — the double-angle +rule as motion, not as a formula. Each frame is a full welded-fault +solve (the probes are measured, not drawn), assembled into a GIF. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from PIL import Image + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +A_RATE, GAMMA = 0.5, 1.0 +R_ANALYTIC = common.ETA * np.sqrt(4 * A_RATE**2 + GAMMA**2) +STEP = 7.5 +angles = np.arange(0.0, 180.0 + 1e-9, STEP) + +# The sweep is 25 welded-fault solves; cache it so iterating on the +# animation's look does not re-measure the physics. +cache = os.path.join(D, "_mohr_probes.npz") +if os.path.exists(cache): + probes = np.load(cache)["probes"] + assert len(probes) == len(angles) + print(f"loaded {len(probes)} cached probes") +else: + probes = [] + for theta in angles: + sigma_n, tau = common.mohr_probe(theta, A_RATE, GAMMA, + half_length=HALF) + probes.append((theta, sigma_n, tau)) + print(f"theta {theta:6.1f}: sigma_n {sigma_n:8.4f} " + f"tau {tau:8.4f}") + probes = np.array(probes) + np.savez(cache, probes=probes) +centre = float(np.mean(probes[:, 1])) +# GEOLOGICAL convention on the stress plane: compression positive. +# The left panel's traction vector stays in physical (tension-positive) +# components; only the Mohr diagram flips. +cg = -centre +scg = -probes[:, 1] + +frames = [] +for k in range(len(probes)): + theta, sig_k, tau_k = probes[k] + fig, (axl, axr) = plt.subplots( + 1, 2, figsize=(9.6, 4.6), + gridspec_kw=dict(width_ratios=[1.0, 1.25])) + + # ---- left: the fault in the box ------------------------------------- + axl.add_patch(plt.Rectangle((0, 0), 1, 1, fill=False, lw=1.0, + edgecolor="0.4")) + t = np.array([np.cos(np.radians(theta)), np.sin(np.radians(theta))]) + n = np.array([-t[1], t[0]]) + c = common.CENTRE + axl.plot([c[0] - HALF * t[0], c[0] + HALF * t[0]], + [c[1] - HALF * t[1], c[1] + HALF * t[1]], + "-", color="#c62828", lw=2.5) + axl.annotate("", xytext=c, xy=c + 0.14 * n, + arrowprops=dict(arrowstyle="->", lw=1.0, color="0.35")) + axl.text(*(c + 0.18 * n), r"$\hat n$", fontsize=9, ha="center", + color="0.35") + # the traction vector on the fault plane, from the MEASURED probe: + # purely normal (aligned with n-hat) exactly at the principal + # orientations — the same instant the probe crosses tau = 0 + T = sig_k * n + tau_k * t + scale = 0.16 / R_ANALYTIC + axl.annotate("", xytext=c, xy=c + scale * T, + arrowprops=dict(arrowstyle="-|>", lw=2.2, + color="#4a7bf7")) + axl.text(0.06, 0.9, rf"$\theta = {theta:.1f}°$", fontsize=12) + axl.text(0.06, 0.135, r"$\sigma\cdot\hat n$: traction on the plane", + fontsize=9, color="#4a7bf7", transform=axl.transAxes) + if abs(tau_k) < 0.05 * R_ANALYTIC: + axl.text(0.5, 0.06, "principal orientation: traction ∥ normal", + fontsize=10, color="#c62828", ha="center", + transform=axl.transAxes) + axl.set_xlim(-0.06, 1.06) + axl.set_ylim(-0.06, 1.06) + axl.set_aspect("equal") + axl.set_xticks([]) + axl.set_yticks([]) + axl.set_title("the (welded) fault rotates ...", fontsize=10) + + # ---- right: the probe sweeps the circle, twice as fast -------------- + tt = np.linspace(0, 2 * np.pi, 300) + axr.plot(cg + R_ANALYTIC * np.cos(tt), R_ANALYTIC * np.sin(tt), + "-", color="0.75", lw=0.9) + axr.axhline(0, color="0.8", lw=0.6) + axr.axvline(cg, color="0.8", lw=0.6) + axr.plot([cg - R_ANALYTIC, cg + R_ANALYTIC], [0, 0], "D", + ms=5, color="0.3", zorder=4) + axr.text(cg + R_ANALYTIC, -0.16 * R_ANALYTIC, + "principal\nstresses", fontsize=7, ha="center", va="top", + color="0.3") + axr.plot(scg[:k + 1], probes[:k + 1, 2], "o", ms=5, + mfc="none", mec="#c62828", mew=1.2) + axr.plot([cg, -sig_k], [0.0, tau_k], "-", + color="#4a7bf7", lw=1.2) + axr.plot([-sig_k], [tau_k], "o", ms=9, color="#c62828", zorder=5) + axr.text(0.04, 0.94, r"... its stress probe sweeps at $2\theta$", + fontsize=10, transform=axr.transAxes) + axr.set_xlabel(r"normal stress $\sigma$ (compression positive)") + axr.set_ylabel(r"shear traction $\tau$") + axr.set_xlim(cg - 1.35 * R_ANALYTIC, cg + 1.35 * R_ANALYTIC) + axr.set_ylim(-1.35 * R_ANALYTIC, 1.35 * R_ANALYTIC) + axr.set_aspect("equal") + + fig.suptitle("Building the Mohr circle with a rotating fault", + fontsize=11) + fig.tight_layout() + frame = os.path.join(D, f"_mohr_frame_{k:03d}.png") + fig.savefig(frame, dpi=110) + plt.close(fig) + frames.append(frame) + +images = [Image.open(f) for f in frames] +out = os.path.join(D, "mohr-circle-build.gif") +images[0].save(out, save_all=True, append_images=images[1:] + + [images[-1]] * 6, # hold the completed circle + duration=280, loop=0) +print("wrote", out, f"({len(frames)} frames)") diff --git a/docs/advanced/figures/fault-examples/mohr_circle.py b/docs/advanced/figures/fault-examples/mohr_circle.py new file mode 100644 index 000000000..f670ea0f2 --- /dev/null +++ b/docs/advanced/figures/fault-examples/mohr_circle.py @@ -0,0 +1,87 @@ +"""The Mohr circle, measured by faults. + +A WELDED split fault transmits the full stress across itself, so it is +a passive stress probe: the no-opening constraint's reaction gives the +normal traction sigma_n, and the stiff interface dashpot reads the +shear traction off its own law, tau = eta_f V. Sweeping the fault +orientation theta and plotting (−sigma_n, |tau|) traces the Mohr +circle of the ambient stress state — measured by the machinery itself, +against the analytic circle of the imposed flow. + +Boundary conditions: Dirichlet velocity on all four walls imposing +the homogeneous flow v = (a(x−c) + gamma(y−c), −a(y−c)) — a PURE-SHEAR +(irrotational stretching) part a plus a SIMPLE-SHEAR part gamma. +Stress sees only the symmetric gradient, so the state is equivalent to +pure shear of magnitude R = eta sqrt(4 a^2 + gamma^2) with principal +axes at 22.5 degrees to the box (deliberately not axis-aligned: the +circle's orientation must be measured, not guessed). The linear flow +is an exact homogeneous Stokes solution — constant stress, no wall +boundary layers — and the centre sits at the (gauge-fixed) mean +pressure, which the fit reports rather than assumes. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +A_RATE, GAMMA = 0.5, 1.0 +ETA_WELD = 200.0 * common.ETA / HALF # welded: slip ~ tau/eta_f -> 0 +R_ANALYTIC = common.ETA * np.sqrt(4 * A_RATE**2 + GAMMA**2) + +angles = np.arange(0.0, 180.0 - 1e-9, 22.5) +points = [] +for theta in angles: + sigma_n, tau = common.mohr_probe(theta, A_RATE, GAMMA, + eta_weld=ETA_WELD, half_length=HALF) + points.append((theta, sigma_n, tau)) + print(f"theta {theta:5.1f}: sigma_n {sigma_n:8.4f} tau {tau:7.4f}") + +points = np.array(points) +# GEOLOGICAL sign convention for the stress plane: compression positive, +# tension on the negative axis (the solver's tractions are +# tension-positive; only the plot flips). +points[:, 1] *= -1.0 + +# circle fit: centre c on the sigma axis, radius R, from the probes +sig, tau = points[:, 1], points[:, 2] +c0 = float(np.mean(sig)) +for _ in range(60): # two-parameter Gauss fit + r = np.sqrt((sig - c0) ** 2 + tau ** 2) + R0 = r.mean() + c0 -= float(np.mean((c0 - sig) * (1 - R0 / np.maximum(r, 1e-30)))) +print(f"fit: centre {c0:.4f}, radius {R0:.4f} " + f"(analytic radius {R_ANALYTIC:.4f})") + +fig, ax = plt.subplots(figsize=(6.8, 4.8)) +tt = np.linspace(0, 2 * np.pi, 300) +ax.plot(c0 + R_ANALYTIC * np.cos(tt), R_ANALYTIC * np.sin(tt), "k-", + lw=1.0, label=f"analytic circle, $R = {R_ANALYTIC:.3f}$") +ax.plot(c0 + R0 * np.cos(tt), R0 * np.sin(tt), "--", lw=0.9, + color="#4a7bf7", label=f"fit through the probes, $R = {R0:.3f}$") +ax.plot(sig, tau, "o", ms=6, color="#c62828", zorder=5, + label="welded-fault probes") +ax.axvline(c0, color="0.6", lw=0.6) +for theta, sg, tu in points: + ax.annotate(f"{theta:.0f}°", (sg, tu), textcoords="offset points", + xytext=(6, 5), fontsize=8) +ax.axhline(0, color="0.6", lw=0.6) +ax.set_xlabel(r"normal stress $\sigma$ (compression positive)") +ax.set_ylabel(r"shear traction $\tau$") +ax.set_title("The Mohr circle, measured by welded split-node faults") +ax.set_aspect("equal") +ax.legend(fontsize=8, loc="center") +fig.tight_layout() +out = os.path.join(D, "mohr-circle.png") +fig.savefig(out, dpi=200) +print("wrote", out) + +assert abs(R0 - R_ANALYTIC) < 0.06 * R_ANALYTIC, "radius off by > 6%" diff --git a/docs/advanced/figures/fault-examples/mohr_cohesion.py b/docs/advanced/figures/fault-examples/mohr_cohesion.py new file mode 100644 index 000000000..720f9870d --- /dev/null +++ b/docs/advanced/figures/fault-examples/mohr_cohesion.py @@ -0,0 +1,211 @@ +"""Cohesive Mohr-Coulomb: strength survives into the tensile sector. + +The cohesive sequel to mohr_friction: the fault's yield envelope is +tau = +-(C + mu sigma) (compression positive), declining through mild +tension and reaching ZERO at sigma = -C/mu. Three physical regimes and +one UNPHYSICAL one appear as the fault rotates: + +- stuck (on the circle) — now around BOTH principal poles, because + cohesion holds shear through mild tension; +- sliding, pinned to the declining envelope; +- HELD SHUT — beyond sigma = -C/mu the strength is zero and the normal + stress is tensile: a real fault would OPEN, no static solution + exists, and the bilateral no-opening constraint manufactures one by + gluing the surfaces together (a tensile constraint reaction). The + solver converges; the physics has failed. Detect it from the SIGN of + the recovered normal traction. + +The law is not a canned option — it is registered as a sympy +expression in the canonical symbols, which is the whole design: a new +fault rheology is four lines. + +All Mohr figures on this page use the GEOLOGICAL sign convention: +compression positive, tension on the negative axis. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import sympy +from PIL import Image + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +A_RATE, GAMMA = 0.5, 1.0 +MU, C, V0 = 0.6, 0.6, 1e-4 +R_ANALYTIC = common.ETA * np.sqrt(4 * A_RATE**2 + GAMMA**2) +STEP = 7.5 +angles = np.arange(0.0, 180.0 + 1e-9, STEP) + + +def register_cohesive_law(stokes): + """Mohr-Coulomb with cohesion as a symbolic law. normal_stress is + the SIGNED effective normal stress (reaction-fed; negative in + tension), and the strength clamp is the law's own: Max(0, C + mu s) + declines through mild tension and vanishes at s = -C/mu.""" + V = fault_contact.slip_rate + S = fault_contact.normal_stress + law = fault_contact.SymbolicFaultLaw( + sympy.Max(C + MU * S, 0) * (2 / sympy.pi) * sympy.atan(V / V0)) + fault_contact.add_frictionless_fault_bc(stokes, "Fault") + fault_contact._register_law(stokes, "Fault", law) + + +cache = os.path.join(D, "_mohr_cohesion_probes.npz") +if os.path.exists(cache): + probes = np.load(cache)["probes"] + assert len(probes) == len(angles) + print(f"loaded {len(probes)} cached probes") +else: + rows = [] + for theta in angles: + child = common.split_with_fault( + common.base_mesh(0.04), common.fault_segment(theta, HALF)) + stokes = common.stokes_on( + child, common.shear_plus_stretch(child, A_RATE, GAMMA)) + register_cohesive_law(stokes) + fault_contact.solve_with_fault(stokes, picard=3) + s, V, _leak = common.slip_profile(stokes) + s_n, sig = common.normal_traction(stokes) + v_med = float(np.median(V[common.inner(s)])) + sigma_n = float(np.median(sig[common.inner(s_n)])) + strength = max(C + MU * (-sigma_n), 0.0) + tau = strength * (2 / np.pi) * np.arctan(v_med / V0) + rows.append((theta, sigma_n, tau, v_med)) + print(f"theta {theta:6.1f}: sigma_n {sigma_n:8.4f} " + f"tau {tau:8.4f} V {v_med:9.5f}") + probes = np.array(rows) + np.savez(cache, probes=probes) + +# GEO convention for plotting: compression positive +sc = -probes[:, 1] +tau = probes[:, 2] +# beyond sigma = -C/mu the strength is zero and the fault is in +# tension: it would open — the constraint holds it shut (unphysical) +held_shut = sc < -C / MU + 1e-6 +sliding = (np.abs(probes[:, 3]) > 5 * V0) & ~held_shut + + +def draw_stress_plane(ax): + tt = np.linspace(0, 2 * np.pi, 300) + ax.plot(R_ANALYTIC * np.cos(tt), R_ANALYTIC * np.sin(tt), + "-", color="0.8", lw=0.9, + label="ambient stress (welded circle)") + s_env = np.linspace(-C / MU, 1.5 * R_ANALYTIC, 60) + for sgn in (+1, -1): + ax.plot(s_env, sgn * (C + MU * s_env), "--", color="0.35", + lw=1.0, label=(r"envelope $\tau = \pm(C + \mu\sigma)$," + r" zero at $\sigma = -C/\mu$" + if sgn > 0 else None)) + ax.axvspan(-1.6 * R_ANALYTIC, -C / MU, color="0.92", zorder=0) + ax.text(-1.5 * R_ANALYTIC, 1.1 * R_ANALYTIC, + "fault would open:\nno static solution\n(held shut by the\n" + "no-opening constraint)", fontsize=7.5, va="top", + color="0.35") + ax.axhline(0, color="0.85", lw=0.6) + ax.axvline(0, color="0.85", lw=0.6) + ax.set_xlabel(r"normal stress $\sigma$ (compression positive)") + ax.set_ylabel(r"shear traction $\tau$") + ax.set_xlim(-1.6 * R_ANALYTIC, 1.6 * R_ANALYTIC) + ax.set_ylim(-1.35 * R_ANALYTIC, 1.35 * R_ANALYTIC) + ax.set_aspect("equal") + + +fig, ax = plt.subplots(figsize=(7.6, 5.4)) +draw_stress_plane(ax) +stuck = ~sliding & ~held_shut +ax.plot(sc[stuck], tau[stuck], "o", ms=7, color="#c62828", + label="stuck: on the circle", zorder=5) +ax.plot(sc[sliding], tau[sliding], "s", ms=6, color="#d9960a", + label="sliding: on the envelope", zorder=5) +ax.plot(sc[held_shut], tau[held_shut], "x", ms=8, mew=2.0, + color="0.45", label="held shut (unphysical)", zorder=5) +ax.legend(fontsize=8, loc="lower right") +ax.set_title(rf"Cohesive Mohr-Coulomb fault, $C = {C}$, $\mu = {MU}$") +fig.tight_layout() +out = os.path.join(D, "mohr-cohesion.png") +fig.savefig(out, dpi=200) +print("wrote", out) + +# ---- the animated build ------------------------------------------------------ +frames = [] +for k in range(len(probes)): + theta, sig_k, tau_k, v_k = probes[k] + shut_k = bool(held_shut[k]) + slide_k = (abs(v_k) > 5 * V0) and not shut_k + fig, (axl, axr) = plt.subplots( + 1, 2, figsize=(9.6, 4.6), + gridspec_kw=dict(width_ratios=[1.0, 1.25])) + + axl.add_patch(plt.Rectangle((0, 0), 1, 1, fill=False, lw=1.0, + edgecolor="0.4")) + t = np.array([np.cos(np.radians(theta)), np.sin(np.radians(theta))]) + n = np.array([-t[1], t[0]]) + c = common.CENTRE + axl.plot([c[0] - HALF * t[0], c[0] + HALF * t[0]], + [c[1] - HALF * t[1], c[1] + HALF * t[1]], + "-", color="#c62828", lw=2.5) + T = sig_k * n + tau_k * t + scale = 0.16 / R_ANALYTIC + axl.annotate("", xytext=c, xy=c + scale * T, + arrowprops=dict(arrowstyle="-|>", lw=2.2, + color="#4a7bf7")) + if slide_k: + off = 0.035 * n + sgn = np.sign(v_k) + for pm in (+1, -1): + axl.annotate("", xytext=c + pm * off - pm * sgn * 0.09 * t, + xy=c + pm * off + pm * sgn * 0.09 * t, + arrowprops=dict(arrowstyle="->", lw=1.4, + color="#d9960a")) + axl.text(0.06, 0.9, rf"$\theta = {theta:.1f}°$", fontsize=12) + status, scol = (("HELD SHUT (unphysical)", "0.45") if shut_k + else ("SLIDING", "#d9960a") if slide_k + else ("stuck", "#c62828")) + axl.text(0.06, 0.82, status, fontsize=10, color=scol) + axl.set_xlim(-0.06, 1.06) + axl.set_ylim(-0.06, 1.06) + axl.set_aspect("equal") + axl.set_xticks([]) + axl.set_yticks([]) + axl.set_title(rf"cohesive fault, $C = {C}$, $\mu = {MU}$ ...", + fontsize=10) + + draw_stress_plane(axr) + axr.legend(fontsize=7, loc="lower right") + stuck_prev = (~sliding & ~held_shut)[:k + 1] + slide_prev = sliding[:k + 1] + shut_prev = held_shut[:k + 1] + axr.plot(sc[:k + 1][stuck_prev], tau[:k + 1][stuck_prev], "o", ms=5, + mfc="none", mec="#c62828", mew=1.2) + axr.plot(sc[:k + 1][slide_prev], tau[:k + 1][slide_prev], "s", ms=5, + mfc="none", mec="#d9960a", mew=1.2) + axr.plot(sc[:k + 1][shut_prev], tau[:k + 1][shut_prev], "x", ms=6, + mew=1.6, color="0.45") + mark, mcol = (("x", "0.45") if shut_k else + ("s", "#d9960a") if slide_k else ("o", "#c62828")) + axr.plot([sc[k]], [tau[k]], mark, ms=9, mew=2.2, color=mcol, + zorder=6) + axr.set_title("... strength declines to zero, then the fault " + "would open", fontsize=10) + + fig.suptitle("Cohesion keeps more of the Mohr circle", fontsize=11) + fig.tight_layout() + frame = os.path.join(D, f"_mohrc_frame_{k:03d}.png") + fig.savefig(frame, dpi=110) + plt.close(fig) + frames.append(frame) + +images = [Image.open(f) for f in frames] +out = os.path.join(D, "mohr-cohesion-build.gif") +images[0].save(out, save_all=True, + append_images=images[1:] + [images[-1]] * 6, + duration=280, loop=0) +print("wrote", out, f"({len(frames)} frames)") diff --git a/docs/advanced/figures/fault-examples/mohr_friction.py b/docs/advanced/figures/fault-examples/mohr_friction.py new file mode 100644 index 000000000..a4731bbaf --- /dev/null +++ b/docs/advanced/figures/fault-examples/mohr_friction.py @@ -0,0 +1,204 @@ +"""The Mohr circle meets the friction envelope. + +The frictional sequel to mohr_circle/mohr_animate: the same rotating +fault, but now carrying Coulomb friction with reaction-fed normal +stress. Three regimes appear as the fault rotates: + +- STUCK — the ambient resolved stress lies inside the envelope + |tau| < mu |sigma_n| (compressive side): the fault transmits the + full stress and its probe sits ON the Mohr circle; +- SLIDING — the ambient stress would exceed the envelope: the fault + slips, drops the shear traction to its strength, and the probe is + pinned to the yield line tau = ±mu |sigma_n|; +- HELD SHUT — under tensile normal stress bare friction has no + strength: a real fault would OPEN, no static solution exists, and + the bilateral no-opening constraint manufactures one by gluing the + surfaces (tensile reaction). The solver converges; the physics has + failed. The probes ride the axis at tau = 0, marked as unphysical. + +sigma_n comes from the no-opening constraint's reaction; tau is read +from the Coulomb law at the measured slip rate — exact in both +regimes, because the regularised law IS the traction the fault +carries. Outputs: a static summary figure and the animated build. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from PIL import Image + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +A_RATE, GAMMA = 0.5, 1.0 +MU, V0 = 0.7, 1e-4 +R_ANALYTIC = common.ETA * np.sqrt(4 * A_RATE**2 + GAMMA**2) +STEP = 7.5 +angles = np.arange(0.0, 180.0 + 1e-9, STEP) + +cache = os.path.join(D, "_mohr_friction_probes.npz") +if os.path.exists(cache): + probes = np.load(cache)["probes"] + assert len(probes) == len(angles) + print(f"loaded {len(probes)} cached probes") +else: + rows = [] + for theta in angles: + child = common.split_with_fault( + common.base_mesh(0.04), common.fault_segment(theta, HALF)) + stokes = common.stokes_on( + child, common.shear_plus_stretch(child, A_RATE, GAMMA)) + fault_contact.add_coulomb_fault_bc(stokes, MU, "Fault", + sigma_n="reaction", V0=V0) + fault_contact.solve_with_fault(stokes, picard=3) + s, V, _leak = common.slip_profile(stokes) + s_n, sig = common.normal_traction(stokes) + v_med = float(np.median(V[common.inner(s)])) + sigma_n = float(np.median(sig[common.inner(s_n)])) + sigma_eff = max(-sigma_n, 0.0) + tau = MU * sigma_eff * (2 / np.pi) * np.arctan(v_med / V0) + rows.append((theta, sigma_n, tau, v_med)) + print(f"theta {theta:6.1f}: sigma_n {sigma_n:8.4f} " + f"tau {tau:8.4f} V {v_med:9.5f}") + probes = np.array(rows) + np.savez(cache, probes=probes) + +centre = 0.0 # the welded sweep's measured gauge +# GEOLOGICAL convention on the stress plane: compression positive. +scg = -probes[:, 1] +# bare friction: strength vanishes the moment the stress turns tensile +held_shut = scg < -1e-6 +sliding = (np.abs(probes[:, 3]) > 5 * V0) & ~held_shut + + +def draw_stress_plane(ax): + tt = np.linspace(0, 2 * np.pi, 300) + ax.plot(centre + R_ANALYTIC * np.cos(tt), R_ANALYTIC * np.sin(tt), + "-", color="0.8", lw=0.9, + label="ambient stress (welded circle)") + ss = np.linspace(0, 1.35 * R_ANALYTIC, 50) + for sgn in (+1, -1): + ax.plot(ss, sgn * MU * ss, "--", color="0.35", lw=1.0, + label=(r"friction envelope $\tau = \pm\mu\sigma$" + if sgn > 0 else None)) + ax.axhline(0, color="0.85", lw=0.6) + ax.axvline(centre, color="0.85", lw=0.6) + ax.axvspan(-1.6 * R_ANALYTIC, 0, color="0.92", zorder=0) + ax.text(-1.45 * R_ANALYTIC, 1.1 * R_ANALYTIC, + "fault would open:\nno static solution\n(held shut by the\n" + "no-opening constraint)", fontsize=7.5, va="top", + color="0.35") + ax.set_xlabel(r"normal stress $\sigma$ (compression positive)") + ax.set_ylabel(r"shear traction $\tau$") + ax.set_xlim(centre - 1.5 * R_ANALYTIC, centre + 1.5 * R_ANALYTIC) + ax.set_ylim(-1.35 * R_ANALYTIC, 1.35 * R_ANALYTIC) + ax.set_aspect("equal") + + +# ---- the static summary ----------------------------------------------------- +fig, ax = plt.subplots(figsize=(7.2, 5.4)) +draw_stress_plane(ax) +stuck = ~sliding & ~held_shut +ax.plot(scg[stuck], probes[stuck, 2], "o", ms=7, + color="#c62828", label="stuck: on the circle", zorder=5) +ax.plot(scg[sliding], probes[sliding, 2], "s", ms=6, + color="#d9960a", label="sliding: on the envelope", zorder=5) +ax.plot(scg[held_shut], probes[held_shut, 2], "x", ms=8, mew=2.0, + color="0.45", label="held shut (unphysical)", zorder=5) +ax.legend(fontsize=8, loc="upper left") +ax.set_title(rf"Coulomb fault probes, $\mu = {MU}$: " + "the stress switches to the yield envelope") +fig.tight_layout() +out = os.path.join(D, "mohr-friction.png") +fig.savefig(out, dpi=200) +print("wrote", out) + +# ---- the animated build ------------------------------------------------------ +frames = [] +for k in range(len(probes)): + theta, sig_k, tau_k, v_k = probes[k] + shut_k = bool(held_shut[k]) + slide_k = (abs(v_k) > 5 * V0) and not shut_k + fig, (axl, axr) = plt.subplots( + 1, 2, figsize=(9.6, 4.6), + gridspec_kw=dict(width_ratios=[1.0, 1.25])) + + axl.add_patch(plt.Rectangle((0, 0), 1, 1, fill=False, lw=1.0, + edgecolor="0.4")) + t = np.array([np.cos(np.radians(theta)), np.sin(np.radians(theta))]) + n = np.array([-t[1], t[0]]) + c = common.CENTRE + axl.plot([c[0] - HALF * t[0], c[0] + HALF * t[0]], + [c[1] - HALF * t[1], c[1] + HALF * t[1]], + "-", color="#c62828", lw=2.5) + axl.annotate("", xytext=c, xy=c + 0.14 * n, + arrowprops=dict(arrowstyle="->", lw=1.0, color="0.35")) + axl.text(*(c + 0.18 * n), r"$\hat n$", fontsize=9, ha="center", + color="0.35") + T = sig_k * n + tau_k * t + scale = 0.16 / R_ANALYTIC + axl.annotate("", xytext=c, xy=c + scale * T, + arrowprops=dict(arrowstyle="-|>", lw=2.2, + color="#4a7bf7")) + if slide_k: + # half-arrows for the slip sense + off = 0.035 * n + sgn = np.sign(v_k) + for pm in (+1, -1): + axl.annotate("", xytext=c + pm * off - pm * sgn * 0.09 * t, + xy=c + pm * off + pm * sgn * 0.09 * t, + arrowprops=dict(arrowstyle="->", lw=1.4, + color="#d9960a")) + axl.text(0.06, 0.9, rf"$\theta = {theta:.1f}°$", fontsize=12) + status, scol = (("HELD SHUT (unphysical)", "0.45") if shut_k + else ("SLIDING", "#d9960a") if slide_k + else ("stuck", "#c62828")) + axl.text(0.06, 0.82, status, fontsize=10, color=scol) + axl.text(0.06, 0.135, r"$\sigma\cdot\hat n$: traction on the plane", + fontsize=9, color="#4a7bf7", transform=axl.transAxes) + axl.set_xlim(-0.06, 1.06) + axl.set_ylim(-0.06, 1.06) + axl.set_aspect("equal") + axl.set_xticks([]) + axl.set_yticks([]) + axl.set_title(rf"Coulomb fault, $\mu = {MU}$, rotating ...", + fontsize=10) + + draw_stress_plane(axr) + axr.legend(fontsize=7, loc="upper left") + stuck_prev = (~sliding & ~held_shut)[:k + 1] + slide_prev = sliding[:k + 1] + shut_prev = held_shut[:k + 1] + axr.plot(scg[:k + 1][stuck_prev], probes[:k + 1][stuck_prev, 2], + "o", ms=5, mfc="none", mec="#c62828", mew=1.2) + axr.plot(scg[:k + 1][slide_prev], probes[:k + 1][slide_prev, 2], + "s", ms=5, mfc="none", mec="#d9960a", mew=1.2) + axr.plot(scg[:k + 1][shut_prev], probes[:k + 1][shut_prev, 2], + "x", ms=6, mew=1.6, color="0.45") + mark, mcol = (("x", "0.45") if shut_k else + ("s", "#d9960a") if slide_k else ("o", "#c62828")) + axr.plot([-sig_k], [tau_k], mark, ms=9, mew=2.2, color=mcol, + zorder=6) + axr.set_title("... its traction cannot leave the envelope", + fontsize=10) + + fig.suptitle("A frictional fault against the Mohr circle", + fontsize=11) + fig.tight_layout() + frame = os.path.join(D, f"_mohrf_frame_{k:03d}.png") + fig.savefig(frame, dpi=110) + plt.close(fig) + frames.append(frame) + +images = [Image.open(f) for f in frames] +out = os.path.join(D, "mohr-friction-build.gif") +images[0].save(out, save_all=True, + append_images=images[1:] + [images[-1]] * 6, + duration=280, loop=0) +print("wrote", out, f"({len(frames)} frames)") diff --git a/docs/advanced/figures/fault-examples/mohr_graded.py b/docs/advanced/figures/fault-examples/mohr_graded.py new file mode 100644 index 000000000..c2d1fde97 --- /dev/null +++ b/docs/advanced/figures/fault-examples/mohr_graded.py @@ -0,0 +1,130 @@ +"""The graded fault: depth-dependent strength sampled along one fault. + +Gravity joins the Mohr experiment. With constant density and closed +(velocity-Dirichlet) walls the flow is untouched — pressure absorbs the +body force exactly — but the WELDED fault's per-node probes now sample +the hydrostatic gradient: every node sits at its own depth, so each +fault orientation contributes a horizontal STREAK of points in the +(sigma, tau) plane rather than a single value. The streak is longest +for the vertical fault (largest depth range along the fault), a single +dot for the horizontal one, and the whole cloud scatters about the +family of Mohr circles between the shallowest and deepest fault points. + +The per-node sampling is native: sigma_n comes from the no-opening +reaction de-smeared NODE BY NODE, and tau from the weld's own law +tau = eta_f V(s) at each node — nothing is averaged. + +Gauge: a closed box fixes pressure only up to a constant (the solver +uses a mean-zero gauge). The plot re-anchors it so p = 0 at the top +surface — the shift is exactly rho g H / 2, known analytically, and is +applied to the plotted sigma only. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +import underworld3 as uw +from underworld3.utilities import fault_contact + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 +A_RATE, GAMMA = 0.5, 1.0 +RHO_G = 0.75 # modest: streaks visible, circle intact +ETA_WELD = 200.0 * common.ETA / HALF +R_ANALYTIC = common.ETA * np.sqrt(4 * A_RATE**2 + GAMMA**2) + +angles = np.arange(0.0, 180.0 - 1e-9, 22.5) + +cache = os.path.join(D, "_mohr_graded_probes.npz") +if os.path.exists(cache): + d = np.load(cache) + theta_all, y_all, sig_all, tau_all = (d["theta"], d["y"], d["sig"], + d["tau"]) + print(f"loaded {len(sig_all)} cached node probes") +else: + theta_all, y_all, sig_all, tau_all = [], [], [], [] + for theta in angles: + child = common.split_with_fault( + common.base_mesh(0.04), common.fault_segment(theta, HALF)) + stokes = common.stokes_on( + child, common.shear_plus_stretch(child, A_RATE, GAMMA)) + stokes.bodyforce = [0.0, -RHO_G] + stokes.add_fault_bc(ETA_WELD, boundary="Fault") + fault_contact.solve_with_fault(stokes, picard=2) + + # per-NODE probes, matched by the along-fault coordinate + coords, jumps, normals = fault_contact.fault_pair_jumps( + stokes, "Fault", stokes._rotated_freeslip_info) + t_hat = np.array([np.cos(np.radians(theta)), + np.sin(np.radians(theta))]) + s_v = (coords - common.CENTRE) @ t_hat + order_v = np.argsort(s_v) + V = (jumps @ t_hat)[order_v] + y = coords[order_v, 1] + s_n, sig = common.normal_traction(stokes) + assert len(sig) == len(V), "pair-node sets disagree" + mid = common.inner(s_v[order_v]) + theta_all.extend([theta] * int(mid.sum())) + y_all.extend(y[mid]) + sig_all.extend(sig[mid]) + tau_all.extend(ETA_WELD * V[mid]) + print(f"theta {theta:6.1f}: {int(mid.sum())} nodes, " + f"sigma_n [{sig[mid].min():.3f}, {sig[mid].max():.3f}]") + theta_all, y_all, sig_all, tau_all = (np.array(theta_all), + np.array(y_all), + np.array(sig_all), + np.array(tau_all)) + np.savez(cache, theta=theta_all, y=y_all, sig=sig_all, tau=tau_all) + +# GEO convention + the top-zero pressure anchor (mean-zero solver gauge +# shifted by the known rho g H / 2) +sc = -sig_all + RHO_G * 0.5 +depth = 1.0 - y_all + +fig, ax = plt.subplots(figsize=(7.8, 5.6)) +# the family of circles between the shallowest and deepest fault points: +# fault nodes span y in [0.3, 0.7] -> pressure rho g (1 - y) +tt = np.linspace(0, 2 * np.pi, 300) +for pshift, style in ((RHO_G * 0.3, ":"), (RHO_G * 0.5, "-"), + (RHO_G * 0.7, ":")): + ax.plot(pshift + R_ANALYTIC * np.cos(tt), R_ANALYTIC * np.sin(tt), + style, color="0.75", lw=0.9, + label=("Mohr circles: shallowest / centre / deepest" + if style == "-" else None)) +ax.axhline(0, color="0.85", lw=0.6) +ax.axvline(0, color="0.85", lw=0.6) + +pts = ax.scatter(sc, tau_all, c=depth, cmap="viridis", s=22, zorder=5, + edgecolors="none") +fig.colorbar(pts, ax=ax, label="depth below surface", shrink=0.8) + +# annotate one streak: the vertical fault has the largest depth range +vert = theta_all == 90.0 +if vert.any(): + ax.annotate("one fault (90°):\na streak, not a point", + xy=(sc[vert].min(), tau_all[vert].mean()), + xytext=(-1.30, -1.30), fontsize=8, ha="left", + arrowprops=dict(arrowstyle="->", lw=0.8, color="0.3")) +flat = theta_all == 0.0 +if flat.any(): + ax.annotate("horizontal fault (0°):\none depth, one point", + xy=(sc[flat].mean(), tau_all[flat].mean()), + xytext=(-1.15, 1.22), fontsize=8, ha="center", + arrowprops=dict(arrowstyle="->", lw=0.8, color="0.3")) + +ax.set_xlabel(r"normal stress $\sigma$ (compression positive, " + r"$p = 0$ at the surface)") +ax.set_ylabel(r"shear traction $\tau$") +ax.set_title(rf"The graded fault: hydrostatic load $\rho g = {RHO_G}$, " + "per-node probes") +ax.set_aspect("equal") +ax.legend(fontsize=8, loc="center") +fig.tight_layout() +out = os.path.join(D, "mohr-graded.png") +fig.savefig(out, dpi=200) +print("wrote", out) diff --git a/docs/advanced/figures/fault-examples/orientations.png b/docs/advanced/figures/fault-examples/orientations.png new file mode 100644 index 000000000..ae9ebd2a1 Binary files /dev/null and b/docs/advanced/figures/fault-examples/orientations.png differ diff --git a/docs/advanced/figures/fault-examples/orientations.py b/docs/advanced/figures/fault-examples/orientations.py new file mode 100644 index 000000000..4b2e9c148 --- /dev/null +++ b/docs/advanced/figures/fault-examples/orientations.py @@ -0,0 +1,62 @@ +"""Fault orientation vs slip: the Mohr circle's other face. + +The same orientation sweep as the Mohr demo, but with FRICTIONLESS +faults under pure simple shear: each fault drops the shear stress +resolved on its own plane, so the peak slip rate follows the resolved +shear tau(theta) = tau_infty cos(2 theta) — the slip-rate version of +reading the Mohr circle. Faults near 45 degrees to the shear plane +barely slip; the fault aligned with it slips fully. +""" +import os + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +import underworld3 as uw + +import common + +D = os.path.dirname(os.path.abspath(__file__)) +HALF = 0.2 + +angles = np.arange(0.0, 90.0 + 1e-9, 15.0) +peaks, profiles = [], [] +for theta in angles: + child = common.split_with_fault( + common.base_mesh(0.04), common.fault_segment(theta, HALF)) + stokes = common.stokes_on(child, common.simple_shear(child)) + stokes.add_fault_bc(0, boundary="Fault") + stokes.solve(verbose=False) + s, V, leak = common.slip_profile(stokes) + assert np.abs(leak).max() < 1e-10 + peaks.append(np.abs(V).max()) + profiles.append((theta, s - s.min() - HALF, np.abs(V))) + print(f"theta {theta:5.1f}: peak slip {peaks[-1]:.4f}") + +peaks = np.array(peaks) + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.4, 4.4)) +cmap = plt.cm.magma +for (theta, s, V) in profiles: + ax1.plot(s, V, ".-", ms=3, lw=0.9, color=cmap(theta / 120.0), + label=f"{theta:.0f}°") +ax1.set_xlabel("position along the fault $s$") +ax1.set_ylabel("slip rate $|V(s)|$") +ax1.set_title("frictionless slip profiles by orientation") +ax1.legend(fontsize=8, title="fault angle", ncol=2) + +tt = np.linspace(0, 90, 200) +ax2.plot(tt, np.abs(np.cos(np.radians(2 * tt))) * peaks[0], "k-", + lw=1.0, label=r"$|\cos 2\theta|$ (resolved shear)") +ax2.plot(angles, peaks, "o", ms=6, color="#c62828", + label="measured peak slip") +ax2.set_xlabel(r"fault angle $\theta$ to the shear plane") +ax2.set_ylabel("peak slip rate") +ax2.set_title("peak slip follows the resolved shear stress") +ax2.legend(fontsize=8) +fig.tight_layout() +out = os.path.join(D, "orientations.png") +fig.savefig(out, dpi=200) +print("wrote", out) diff --git a/docs/advanced/figures/fault-examples/true-branch-stress.png b/docs/advanced/figures/fault-examples/true-branch-stress.png new file mode 100644 index 000000000..d33cdb1f2 Binary files /dev/null and b/docs/advanced/figures/fault-examples/true-branch-stress.png differ diff --git a/docs/advanced/figures/fault-examples/true-branch.png b/docs/advanced/figures/fault-examples/true-branch.png new file mode 100644 index 000000000..e8e918d00 Binary files /dev/null and b/docs/advanced/figures/fault-examples/true-branch.png differ diff --git a/docs/advanced/figures/fault-examples/true_branch.py b/docs/advanced/figures/fault-examples/true_branch.py new file mode 100644 index 000000000..8036821e0 --- /dev/null +++ b/docs/advanced/figures/fault-examples/true_branch.py @@ -0,0 +1,240 @@ +"""A true Y-branch, bracketed: does the near-miss tributary capture it? + +A genuine branch (three arms meeting at a point) cannot be split — but +it can be DECOMPOSED two ways, and the pair brackets the truth: + +- A: the TRUNK is continuous (west + east arms as one fault); the + splay abuts, stopping a ligament short. +- B: the BENT fault is continuous (west arm + splay as one deliberately + kinked polyline — the kink response is the physics, so no smoothed + normal); the east arm abuts. + +Each decomposition welds a different pair of arms through the junction +and feeds the third across a ligament. If both give the same slip on +every arm at a small gap, the offset representation reproduces the true +branch; sweeping the gap in decomposition A shows how close "close by" +has to be. All three arms slip freely under the same drive. +""" +import os +import time + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv + +import underworld3 as uw +from underworld3.meshing.surfaces import prepare_fault_network +from underworld3.utilities import fault_contact + +import common + +pv.OFF_SCREEN = True +D = os.path.dirname(os.path.abspath(__file__)) +H = 0.012 +MU_P = 0.4 +TREND = np.degrees(np.arctan2(0.10, 0.70)) + +J = np.array([0.50, 0.50]) # the branch point +W_END = np.array([0.15, 0.45]) +E_END = np.array([0.85, 0.55]) +S_END = np.array([0.80, 0.76]) +ARMS = {"west": (W_END - J) / np.linalg.norm(W_END - J), + "east": (E_END - J) / np.linalg.norm(E_END - J), + "splay": (S_END - J) / np.linalg.norm(S_END - J)} + + +def decomposition(kind, lig): + if kind == "A": # trunk continuous + faults = [("Trunk", np.array([W_END, E_END])), + ("Splay", np.array([J, S_END]))] + through = ["Trunk"] + else: # bent fault continuous + faults = [("Bent", np.array([W_END, J, S_END])), + ("East", np.array([J, E_END]))] + through = ["Bent"] + return prepare_fault_network(faults, spacing=H, ligament=lig, + through=through, verbose=False) + + +ETA_WELD = 200.0 * common.ETA / 0.2 + + +def _stress_p0(child, stokes, tagc): + x, y = child.X + v, p = stokes.Unknowns.u, stokes.Unknowns.p + comps = {} + for cname, expr in ( + ("sxx", -p.sym[0] + 2 * common.ETA * v.sym[0].diff(x)), + ("syy", -p.sym[0] + 2 * common.ETA * v.sym[1].diff(y)), + ("sxy", common.ETA * (v.sym[0].diff(y) + v.sym[1].diff(x)))): + s_var = uw.discretisation.MeshVariable(f"{cname}_{tagc}", child, 1, + degree=0, continuous=False) + proj = uw.systems.Projection(child, s_var) + proj.uw_function = expr + proj.smoothing = 0.0 + proj.solve() + row = common.split_mesh_cell_rows(child, s_var) + comps[cname] = np.asarray(s_var.data[:, 0])[row].copy() + return comps + + +def run(tag, prepared, want_field=False): + child = common.base_mesh(H).add_fault(prepared) + stokes = common.stokes_on( + child, common.boundary_simple_shear(child, TREND)) + for n, _p in prepared: + stokes.add_fault_bc(0, boundary=n) + t0 = time.perf_counter() + fault_contact.solve_with_fault(stokes, picard=2) + + field = None + if want_field: + c1 = _stress_p0(child, stokes, f"a{run.n}") + s0 = common.stokes_on( + child, common.boundary_simple_shear(child, TREND)) + for n, _p in prepared: + s0.add_fault_bc(ETA_WELD, boundary=n) + fault_contact.solve_with_fault(s0, picard=2) + c0 = _stress_p0(child, s0, f"b{run.n}") + run.n += 1 + beta = np.radians(TREND) + nx, ny = -np.sin(beta), np.cos(beta) + tx, ty = np.cos(beta), np.sin(beta) + + def resolve(c): + s_nn = (c["sxx"] * nx * nx + 2 * c["sxy"] * nx * ny + + c["syy"] * ny * ny) + s_t = (c["sxx"] * tx * nx + c["sxy"] * (tx * ny + ty * nx) + + c["syy"] * ty * ny) + return s_nn, s_t + + nn0, t_0 = resolve(c0) + nn1, t_1 = resolve(c1) + dcff = np.sign(np.median(t_0)) * (t_1 - t_0) + MU_P * (nn1 - nn0) + pts, faces = common.split_mesh_cell_render(child) + fc = np.asarray(faces).reshape(-1, 4)[:, 1:] + cent = np.asarray(pts)[fc].mean(axis=1) + dcff, gauge = common.far_field_anchor( + cent, dcff, [p for _n, p in prepared], cut=0.18) + print(f"[{tag}] gauge {gauge:+.4f}", flush=True) + field = (pts, faces, dcff) + arms = {k: ([], []) for k in ARMS} + for n, _p in prepared: + coords, jumps, normals = fault_contact.fault_pair_jumps( + stokes, n, stokes._rotated_freeslip_info) + if not len(coords): + continue + tang = np.column_stack([-normals[:, 1], normals[:, 0]]) + V = np.abs(np.einsum("ij,ij->i", jumps, tang)) + R = coords - J + # classify each pair node by the arm it lies along + proj = {k: R @ t for k, t in ARMS.items()} + dist2 = {k: np.einsum("ij,ij->i", R - np.outer( + np.clip(proj[k], 0, None), ARMS[k]), + R - np.outer(np.clip(proj[k], 0, None), ARMS[k])) + for k in ARMS} + best = np.argmin(np.vstack([dist2[k] for k in ARMS]), axis=0) + for i, k in enumerate(ARMS): + sel = best == i + arms[k][0].extend(proj[k][sel]) + arms[k][1].extend(V[sel]) + out = {} + for k in ARMS: + s = np.asarray(arms[k][0]) + v = np.asarray(arms[k][1]) + order = np.argsort(s) + out[k] = (s[order], v[order]) + print(f"[{tag}] peaks " + " ".join( + f"{k} {out[k][1].max():.4f}" for k in ARMS) + + f" ({time.perf_counter() - t0:.1f} s)", flush=True) + return out, field, prepared + + +run.n = 0 + +cases = [("A, gap 1h", "A", 1.0, "#1565c0", "-"), + ("A, gap 2h", "A", 2.0, "#64b5f6", "--"), + ("A, gap 4h", "A", 4.0, "#b3d7f7", ":"), + ("B, gap 1h", "B", 1.0, "#c62828", "-")] +profiles, fields = {}, {} +for tag, kind, lig, _c, _ls in cases: + prepared, _rep = decomposition(kind, lig) + want = tag in ("A, gap 1h", "B, gap 1h") + profiles[tag], field, prep = run(tag, prepared, want_field=want) + if field is not None: + fields[tag] = (field, prep) + +fig, axes = plt.subplots(1, 3, figsize=(13.2, 4.3), sharey=True) +for ax, arm in zip(axes, ("west", "east", "splay")): + for tag, _k, _l, col, ls in cases: + s, v = profiles[tag][arm] + ax.plot(s, v, ls, lw=1.4, color=col, label=tag) + ax.set_title(f"the {arm} arm", fontsize=10.5) + ax.set_xlabel("distance from the branch point") +axes[0].set_ylabel("|slip|") +axes[0].legend(fontsize=8, title="decomposition, ligament") +fig.suptitle( + "One true Y-branch, two decompositions: continuous-trunk (A) vs " + "continuous-bend (B).\nAgreement between A and B at small gap = the " + "near-miss tributary reproduces the true branch", fontsize=10.5) +fig.tight_layout() +out = os.path.join(D, "true-branch.png") +fig.savefig(out, dpi=200) +print("wrote", out, flush=True) + +# ---- the stress maps: the kink-lock made visible --------------------------- + + +def render(field, prepared, png, scale, focal): + pts, faces, dcff = field + pvm = pv.PolyData(np.asarray(pts, dtype=float), + faces=np.asarray(faces, dtype=np.int64)) + pvm.cell_data["dcff"] = dcff + pl = pv.Plotter(off_screen=True, window_size=(850, 800)) + pl.set_background("white") + pl.add_mesh(pvm, scalars="dcff", cmap="RdBu_r", clim=(-1.0, 1.0), + lighting=False, show_scalar_bar=False) + for n, p in prepared: + line = pv.lines_from_points( + np.column_stack([p, np.full(len(p), 1e-3)])) + pl.add_mesh(line, color="black", line_width=4, lighting=False) + pl.view_xy() + pl.camera.parallel_projection = True + pl.camera.parallel_scale = scale + pl.camera.focal_point = (focal[0], focal[1], 0.0) + pl.screenshot(png) + pl.close() + return png + + +panels = [] +for tag, sub in (("A, gap 1h", "a"), ("B, gap 1h", "b")): + field, prep = fields[tag] + panels.append((render(field, prep, os.path.join( + D, f"_tb_map_{sub}.png"), 0.40, (0.5, 0.55)), + f"{tag.split(',')[0]}: " + + ("trunk continuous" if sub == "a" else "bend continuous"))) + panels.append((render(field, prep, os.path.join( + D, f"_tb_zoom_{sub}.png"), 0.09, J), + "the branch point")) + +fig = plt.figure(figsize=(12.8, 11.6)) +gs = fig.add_gridspec(2, 3, width_ratios=[1.25, 1.0, 0.05]) +for k, (png, title) in enumerate(panels): + ax = fig.add_subplot(gs[k // 2, k % 2]) + ax.imshow(plt.imread(png)) + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_title(title, fontsize=10.5) +from matplotlib import cm, colors as mcolors +sm = cm.ScalarMappable(norm=mcolors.Normalize(-1, 1), cmap="RdBu_r") +cax = fig.add_subplot(gs[:, 2]) +fig.colorbar(sm, cax=cax, label=r"$\Delta$CFF (per unit stress drop)") +fig.suptitle("The same Y-branch, two continuity choices: " + r"$\Delta$CFF on trunk-parallel planes", fontsize=12) +fig.tight_layout() +out2 = os.path.join(D, "true-branch-stress.png") +fig.savefig(out2, dpi=200) +print("wrote", out2, flush=True) diff --git a/docs/advanced/index.md b/docs/advanced/index.md index eaf2d99b3..5ede20451 100644 --- a/docs/advanced/index.md +++ b/docs/advanced/index.md @@ -103,6 +103,8 @@ multigrid-preconditioning solver-iteration-callbacks complex-rheologies vep-transverse-isotropy-faults +split-node-faults +fault-mechanics-examples custom-meshes curved-boundary-conditions mesh-adaptation diff --git a/docs/advanced/split-node-faults.md b/docs/advanced/split-node-faults.md new file mode 100644 index 000000000..4eadbeba2 --- /dev/null +++ b/docs/advanced/split-node-faults.md @@ -0,0 +1,300 @@ +--- +title: "Split-Node Faults: Zero-Thickness Frictional Surfaces" +--- + +# Split-Node Faults + +A fault, mechanically, is a surface across which the velocity field jumps: +slip, not strain. Underworld3 can represent a fault exactly that way — the +mesh nodes along a conforming fault surface are duplicated so the two +sides share no degrees of freedom, the two coincident copies are tied by a +strong no-opening constraint, and the tangential jump (the slip rate) +either emerges freely or is governed by a friction law. There is no weak +layer, no viscosity contrast, and no band width to resolve: the fault has +genuinely zero thickness, and the surrounding mesh can be uniform. + +This is one of two fault representations in Underworld3, and they are +complementary rather than competing: + +- **Split-node surface** (this page) — the fault is a sharp interface. + The right tool when slip is the quantity of interest: earthquake-cycle + models, stress transfer, plate-boundary-style discontinuities. +- **[Finite-width TI weak zone](vep-transverse-isotropy-faults.md)** — + the fault is a thin volume with anisotropic rheology. The right tool + when the fault-zone width is itself constitutive (rate-and-state + linked to a physical gouge width), or when the same fault must appear + as a damage zone for another equation (e.g. porous flow). + +## Quick start (2-D) + +A fault is a polyline with both tips strictly inside the domain. One call +places its points onto mesh vertices, cuts the mesh so the fault becomes +a conforming facet chain, and splits it: + +```python +import numpy as np +import underworld3 as uw + +mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.05) +fault_points = np.array([[0.3, 0.4], [0.7, 0.6]]) + +child = mesh.add_fault(("Fault", fault_points)) +``` + +`child` is a new, standalone mesh carrying boundaries `FaultPlus` and +`FaultMinus` — two geometrically coincident copies of the fault — plus +the record of which degree of freedom pairs with which. The source mesh +is untouched: when the fault moves, call `add_fault` again on the base +mesh at the new position (nothing is cumulative). + +A frictionless (perfectly slippery) fault is then a one-line boundary +condition, and the solve is an ordinary `solve()`: + +```python +v = uw.discretisation.MeshVariable("V", child, child.dim, degree=2) +p = uw.discretisation.MeshVariable("P", child, 1, degree=0, + continuous=False) +stokes = uw.systems.Stokes(child, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + +x, y = child.X +for wall in ("Bottom", "Top", "Left", "Right"): + stokes.add_dirichlet_bc((y - 0.5, 0.0), wall) # far-field shear + +stokes.add_fault_bc(0, boundary="Fault") # frictionless +stokes.solve() +``` + +The no-opening constraint $[\mathbf v]\cdot\hat n = 0$ is always imposed +strongly (the measured normal jump is machine zero), and with `conds = 0` +the shear traction on the fault is exactly zero, so slip develops freely +— the stress-driven crack. Read the slip through the pairing: + +```python +from underworld3.utilities import fault_contact + +s, V, leak = fault_contact.fault_slip( + stokes, "Fault", stokes._rotated_freeslip_info) +# s: along-fault coordinate; V: slip rate; leak: normal jump (~1e-17) +``` + +```{warning} +Always read fault quantities through the pairing-based helpers +(`fault_slip`, `fault_pair_jumps`, `fault_normal_traction`). The two +sides of the fault are *geometrically coincident*, so any query by +coordinate — `uw.function.evaluate` included — sees one side only and +silently averages or picks arbitrarily. +``` + +## Quick start (3-D) + +In 3-D the fault is a triangulated patch whose rim stays strictly inside +the domain. `uw.meshing.BoxInternalPatch` embeds a planar polygon patch +conformingly in a simplex box (a disc is a many-sided polygon), and +`split_fault` does the rest: + +The preferred input is a `uw.meshing.FaultSurface` — the object names +the patch, its rim polygon becomes the conforming embed, and it rides +onto the split mesh so the constraint frame can come from the surface's +own face normals: + +```python +from underworld3.utilities.fault_split import split_fault + +fault = uw.meshing.FaultSurface("Fault", points) # or .from_vtk(...) +fault.triangulate() +mesh = uw.meshing.BoxInternalPatch(cellSize=0.05, patch_points=fault, + patch_cellSize=0.02) +child = split_fault(mesh, "Fault") +... +stokes.add_fault_bc(0, boundary="Fault", normal="surface") +``` + +The surface must be planar for now (`rim_polygon` refuses otherwise — +a genuinely curved sheet waits on the discrete-entity embed). A raw +`(N, 3)` polygon via `patch_points=` still works for quick setups: + +```python +patch = np.array([[0.5, 0.3, 0.3], [0.5, 0.7, 0.3], + [0.5, 0.7, 0.7], [0.5, 0.3, 0.7]]) +mesh = uw.meshing.BoxInternalPatch(cellSize=0.05, patch_points=patch, + patch_name="Fault") +child = split_fault(mesh, "Fault") +``` + +Everything downstream is identical to 2-D — `add_fault_bc`, `solve()`, +and the law functions below all work unchanged. In 3-D the slip is an +in-plane vector; read it with `fault_pair_jumps`: + +```python +coords, jumps, normals = fault_contact.fault_pair_jumps( + stokes, "Fault", stokes._rotated_freeslip_info) +leak = np.einsum("ij,ij->i", jumps, normals) # ~ machine zero +slip_vec = jumps - leak[:, None] * normals # the slip vectors +``` + +## The fault laws + +The tangential condition on the fault is an interface constitutive law +$\tau(V, \sigma_n, \theta)$ relating the shear traction to the slip rate +$V$, the effective normal stress $\sigma_n$, and (for rate-and-state) the +state variable $\theta$. Laws are sympy expressions; their consistent +Newton tangents are derived symbolically, so every law converges with +full Newton — none of them is Picard-limited. + +| law | call | condition | +|-----|------|-----------| +| frictionless | `stokes.add_fault_bc(0, boundary="Fault")` | $\tau = 0$ | +| viscous | `stokes.add_fault_bc(eta_f, boundary="Fault")` | $\tau = \eta_f V$ | +| Coulomb | `fault_contact.add_coulomb_fault_bc(stokes, mu, "Fault", sigma_n="reaction", V0=1e-5)` | $\tau = \mu\,\sigma_n\,\tfrac{2}{\pi}\arctan(V/V_0)$ | +| rate-and-state | `fault_contact.add_rate_state_fault_bc(stokes, f0, "Fault", a=..., b=..., V0=..., Dc=...)` | regularised arcsinh form, ageing law | + +Notes on each: + +- **Viscous**: `eta_f` is a viscosity per unit length — the + zero-thickness limit of a shear band is $\eta_f = \eta_{\rm band}/w$ + with the *band's own* (weak-zone) viscosity over its width. Small + `eta_f` approaches the frictionless crack; large `eta_f` *welds* the + fault — and welding recovers the uncut continuum exactly, because the + law penalises only the jump. The natural scale is $\eta/a$ (bulk + viscosity over fault half-length), at which the fault slips at roughly + half its free rate. +- **Coulomb**: `sigma_n` may be a prescribed number or `"reaction"`, in + which case the effective normal stress is recovered from the + no-opening constraint's own reaction — no auxiliary projection, and + compressive stress generated by the flow feeds straight into the + frictional strength. The reaction-fed value is SIGNED: where it turns + tensile the strength is zero, and a real fault would OPEN — the + bilateral no-opening constraint holds it shut instead (a tensile + reaction), so check the sign of `fault_normal_traction` before + trusting results on faults that see tension. `V0` is the regularisation velocity: choose it + well below the slip rates the flow produces. Below it the fault sticks + (creep $\sim V_0$); above it the traction saturates at + $\mu\,\sigma_n$ and the fault slides at constant stress drop. +- **Rate-and-state**: the state $\theta$ advances *between* solves by + the ageing law $\dot\theta = 1 - V\theta/D_c$, integrated exactly over + the interval. Nonlinear fault solves go through `solve_with_fault`, + and a time loop alternates solve and state update: + +```python +fault_contact.add_rate_state_fault_bc( + stokes, 0.6, "Fault", a=0.015, b=0.010, V0=1e-6, Dc=1e-2) + +for step in range(n_steps): + result = fault_contact.solve_with_fault(stokes, picard=2) + monitor = fault_contact.update_fault_state( + stokes, "Fault", dt, solve_result=result) + # monitor -> 1.0 as the fault approaches steady state (theta V = Dc) +``` + +The recovered normal traction along the fault (also the quantity a +Coulomb law reads when `sigma_n="reaction"`): + +```python +s, sigma_nn = fault_contact.fault_normal_traction( + stokes, "Fault", stokes._rotated_freeslip_info) +# negative in compression; in 3-D the first return is pair coordinates +``` + +## Prescribed (kinematic) slip + +Because the two sides are ordinary named boundaries, a *prescribed* slip +distribution is just Dirichlet data — no law involved. The slip must +taper to zero at the tips (2-D) or rim (3-D), which stay unsplit: + +```python +import sympy + +x, y = child.X +u = ((x - cx) * tx + (y - cy) * ty) / half_length # along-fault coord +taper = sympy.sqrt(sympy.Max(1 - u**2, 0)) # elliptical profile + +stokes.add_dirichlet_bc((vx_bg + s0 / 2 * taper * tx, + vy_bg + s0 / 2 * taper * ty), "FaultPlus") +stokes.add_dirichlet_bc((vx_bg - s0 / 2 * taper * tx, + vy_bg - s0 / 2 * taper * ty), "FaultMinus") +``` + +The elliptical taper is the constant-stress-drop crack profile; because +it vanishes at the tips, the shared tip vertex receives the same value +from both sides and no special treatment is needed. + +## Practical notes + +- **Pressure space**: use P2 velocity with P0 *discontinuous* pressure + on split meshes. A continuous pressure space smears the pressure jump + across the fault and measurably pollutes the near-fault stress. +- **Solvers**: split meshes carry no geometric-multigrid tail (the + coarse levels never contain the fault), so the velocity block takes + its algebraic-multigrid default. Set accuracy with + `stokes.tolerance = ...` — never raw `ksp_rtol`, which the solver + configuration overrides silently. +- **Conditioning**: the contact formulation is dramatically better + conditioned than a thin weak inclusion of the same fault — there is + no thin feature and no viscosity contrast for the Schur complement + to fight (measured: 10 outer iterations vs 147 for a $10^{-4}$ + contrast band on the same mesh). +- **Curved faults need their analytic normal**: a curved trace is + sampled as a polyline, and by default each fault node's normal is the + average of its adjacent facet normals — exact on a straight fault, but + zig-zagging at the sampling kinks of a curve. The no-opening + constraint then forbids smooth slip past each kink, producing slip + notches and normal-traction sawteeth that *grow* under mesh + refinement. Pass the smooth curve's normal instead: + + ```python + x, y = mesh.X + # e.g. a circular arc about (cx, cy): the radial direction + stokes.add_fault_bc(0, boundary="Arc", + normal=sympy.Matrix([[x - cx, y - cy]])) + ``` + + Same conventions as `add_rotated_freeslip_bc`: a sympy `1×dim` matrix + in `mesh.X` (need not be unit length; it is normalised per node and + sign-aligned to the split's Plus→Minus orientation), or a constant + array. For a **digitized trace with no analytic formula** — real + mapped faults — pass `normal="trace"`: the smoothed normal is built + from the fault's own control polyline (central-difference tangents at + the control points, tangent angle interpolated along each segment), + which `add_fault` stores on the mesh. Every fault-law variant + (`add_coulomb_fault_bc`, `add_rate_state_fault_bc`, ...) accepts the + same `normal=` argument, and the slip/traction diagnostics read the + same frame automatically. Measured on a sampled circular arc: the + smooth normal reduces the kink sawtooth by an order of magnitude and + restores convergence under refinement. On straight faults it changes + nothing — omit it. And a deliberately *kinked* fault should not be + smoothed: there the kink response is the physics. +- **Moving faults**: re-derive, don't update. Cut and split again from + the static base mesh at the new fault position; transfer fields with + the standard re-adaptation machinery. +- **Networks**: pass a list of faults to `add_fault`. Segments must not + share vertices — represent a branch or crossing as offset segments + separated by a ligament of one or two cell sizes. +- **Parallel**: `add_fault` / `split_fault` REDISTRIBUTE first, in + both 2-D and 3-D: the fault's cell star (plus one growth layer — a + thin skin, not the refined band) is gathered onto the rank that + already owns most of it, everything else stays with the + load-balanced partition, and the split then runs with serial + topology at any rank count. In 2-D the move happens only when the + cut chain actually touches the partition seam, and a network passed + to one `add_fault` call is redistributed ONCE, keyed on all its + faults together — which is also what keeps every fault's pairing + valid, since a pairing cannot yet migrate through a redistribution. + In 3-D the move is unconditional. The cost is a bounded imbalance on + the fault-owning rank (measured ~1.8x at np = 8 on a graded box). + One current exception follows from the pairing rule: adding several + faults ONE SPLIT AT A TIME in parallel is refused when a later fault + needs redistribution — pass the whole network to one `add_fault` + call (2-D), or split 3-D networks in serial. All refusals are + collective and name the actual problem. + +## Current limitations + +Refused loudly rather than mishandled: closed-loop faults (rings, +spheres); faults that reach the domain boundary (daylighting); junctions +sharing vertices; 3-D multi-fault networks in parallel (the pairing does +not yet migrate through the redistribution). The design +documents in `docs/developer/design/` (`SPLIT_NODE_FAULT_METHOD_2026-08` +and `FAULT_CONTACT_DEPLOYMENT_2026-08`) record the method, the +validation benchmarks, and the roadmap for these extensions. diff --git a/docs/developer/ai-notes/fault-teaching-examples-handoff-2026-08.md b/docs/developer/ai-notes/fault-teaching-examples-handoff-2026-08.md new file mode 100644 index 000000000..d7c810ff5 --- /dev/null +++ b/docs/developer/ai-notes/fault-teaching-examples-handoff-2026-08.md @@ -0,0 +1,224 @@ +# Fault-mechanics teaching examples: curation handoff + +Written 2026-08-05 for the session curating the teaching materials. +Everything described here is committed and pushed on +`origin/feature/fault-split-node` (through the commit adding this +file). Read this document first; every path below is repo-relative on +that branch. + +## The one critical dependency + +**These examples run ONLY on `feature/fault-split-node`.** They are +built on the split-node fault capability (`fault_split.py`, +`fault_contact.py`, the pair blocks in `rotated_bc.py`), which is not +on `development` yet. Curate the *placement and prose* freely, but the +examples cannot execute from any other branch, and they should not be +merged into user-facing docs ahead of the capability itself. The +figures, however, are committed PNGs/GIFs — prose pages that embed +them render fine anywhere. + +## Where everything is + +- **Page**: `docs/advanced/fault-mechanics-examples.md` (in the + `docs/advanced` toctree). Nine figures + three animations, each + with a short teaching narrative and a closing "what these models + are, and are not" caveat section. +- **Scripts, figures, caches**: `docs/advanced/figures/fault-examples/` + — one Python script per example, committed output PNG/GIF alongside, + probe caches as `_*.npz` (committed: they let anyone re-style every + figure and both animations without re-running ~100 solves). + `_*.png` frame intermediates and `*.log` files are gitignored. +- **Shared harness**: `figures/fault-examples/common.py` — every + example builds on it (see contracts below). +- **Companion user documentation**: `docs/advanced/split-node-faults.md` + (the API guide the examples assume) and the method/benchmark + write-up `docs/developer/design/SPLIT_NODE_FAULT_METHOD_2026-08.md` + (with its own figure set under + `docs/developer/design/figures/split-node-faults/` — cetz sources, + generators, and PNGs for the split anatomy, pair transform, mesh + stack, and grid hierarchy; those serve the method paper / blog + posts rather than the teaching page). + +## The example inventory (teaching order) + +1. **`ladder.py` → `ladder.png`** — *the fault-strength ladder.* One + fault, one shear drive, five solves down the constitutive ladder + (frictionless / viscous eta/a / Coulomb-weak / rate-state / + Coulomb-stuck), slip profiles against the elliptical crack shape, + plus a stacked column of per-rung shear-stress panels. Teaches: what + an interface law IS. ~5 solves. +2. **`mohr_circle.py` → `mohr-circle.png`** — welded faults as passive + stress probes trace the Mohr circle (fitted radius 1.411 vs analytic + sqrt(2)). Teaches: the probe instrument + the circle exists. +3. **`mohr_animate.py` → `mohr-circle-build.gif`** — the circle built + frame by frame as the fault rotates; the measured traction vector + snaps to the fault normal at the principal orientations, flagged on + both panels. Teaches: the 2-theta rule as motion. 25 solves, cached. +4. **`mohr_friction.py` → `mohr-friction.png` + `-build.gif`** — give + the rotating fault Coulomb friction: stuck probes reproduce the + circle, sliding probes pin to the envelope, tensile probes are + HELD SHUT (unphysical — the bilateral constraint glues an opening + fault; the tell is the sign of the normal traction). Teaches: the + yield envelope truncates the circle. +5. **`mohr_cohesion.py` → `mohr-cohesion.png` + `-build.gif`** — + cohesion: strength declines through mild tension to zero at + sigma = -C/mu; stuck arcs at both poles; the law is registered as a + four-line sympy expression (the extension path for any rheology). +6. **`mohr_graded.py` → `mohr-graded.png`** — hydrostatic load: each + welded fault becomes a depth-coloured STREAK spanning the family of + Mohr circles between its shallowest and deepest points. Teaches: + per-node stress recovery = depth-dependent strength along one fault. +7. **`orientations.py` → `orientations.png`** — frictionless faults at + swept orientations: peak slip follows |cos 2 theta| with an exact + zero at 45 degrees. Teaches: resolved stress controls slip. +8. **`interacting_faults.py` → `interacting-faults.png` + + `interacting-rotation.png`** — King-style stress transfer: the + source slips, the receiver is welded in the source's along-strike + tip lobe. Delta CFF as a field (linear RdBu_r at +-1, P0 cells) + AND as the receiver's probe cloud moving in the Mohr plane, with + the loaded near-tip nodes CROSSING the cohesive envelope into the + shaded failure zone while the far end retreats. The friction + dressing (P0 = 1 confining pressure, C = 0.75 cohesion) enters the + ENVELOPE only, never Delta CFF — both docs say so explicitly. The + rotation sweep's message: the regional stress orientation controls + the MARGIN to the envelope (grazing at phi = 20, crossing at 45, + parked safe at 70). 6 solves at h = 0.016, 13-18 s each. +9. **`california.py` → `california.png`** — schematic southern + California: the San Andreas as ONE continuous dextral trace with a + smooth tanh S-bend (the Big Bend — the smoothed stepover; left + step = restraining). The curved trace is sampled as a polyline and + carries the smooth curve's ANALYTIC NORMAL via + `add_fault_bc(..., normal=...)` — the capability that made curved + traces viable (see the rendering/roughness rules below). The + restraining bend fills with strong compression: a dCFF bowtie + exactly where the Transverse Ranges belong — the model puts the + mountains where the mountains are. Garlock (resolves sinistral + from the kinematics), three ECSZ strands and a San Jacinto + analogue welded as probes, all inboard. Verdicts: SJF relaxed + (-0.41), Garlock and ECSZ mildly loaded (+0.05/+0.06); slip + arrows from the MEASURED jump (+0.242 right-lateral — a + continuous trace slips more than the two offset segments did). + ~50 s of compute at h = 0.012. The capstone. + + RENDERING RULES the interaction fields obey: see the revised + field-rendering entry in the harness contracts below (P0 and P1 + are both legitimate now — choose by physics; never Delaunay a + split-mesh DOF cloud; RdBu_r +-1 for dCFF, Oranges for strain + rate). + +## The harness contracts (`common.py`) + +- `base_mesh(h)` + `mesh.add_fault([...])`: one static base, re-faulted + per case (the non-cumulative pattern). Networks = list form; + segments must not share vertices (ligament >= 2h). +- `stokes_on(child, drive)`: P2 velocity / P0 *discontinuous* pressure + (the fault pressure-space ruling), all-wall Dirichlet drive, + `petsc_use_pressure_nullspace = True`, `stokes.tolerance` (never raw + ksp_rtol). +- Drives: `simple_shear`, `shear_plus_stretch(a, gamma)` (Mohr radius + eta sqrt(4a^2+gamma^2)), `pure_shear_drive(phi)` (phi = COMPRESSION + axis — an earlier sign error made it the extension axis; fixed), + `boundary_simple_shear(trend)` (right-lateral along a plate-boundary + trend; dextral verified from the measured jump). +- **Read fault quantities through the DOF pairing only** + (`slip_vs_position`, `probe_nodes`, `fault_pair_jumps`): the two + sides are geometrically coincident — coordinate queries see one side. + `probe_nodes` selects the named fault's nodes via its own pairing + (several law-carrying faults share one assembler). +- **Pressure-gauge discipline** (closed velocity-driven boxes fix + pressure only up to a per-solve constant; observed drifts up to + ~300 stress units on large slip events): differenced fields are + anchored with `far_field_anchor` (a slip event changes nothing far + away), absolute probe pressures with `ambient_sigma_n` / + `ambient_sigma_n_simple` (analytic ambient states). Every removed + constant is printed, never silently absorbed. +- Field rendering (the standing rules, measured not guessed — REVISED + 2026-08-07): render on the split mesh's TRUE connectivity via + `split_mesh_cell_render`/`split_mesh_cell_rows` (P0 cells) or + `split_mesh_faces` (P1 nodes). Never Delaunay the DOF point cloud of + a split mesh (coincident fault pairs + trace-crossing edges paint + false beading). BOTH projections are now legitimate: the old "never + P1" rule was calibrated against the kinked-normal sawtooth fields we + later cured at the source, and the split mesh's P1 space is + naturally discontinuous across faults (doubled nodes), so the jump + is representable. Choose by physics: P1 nodal for smooth fields + (dCFF away from tips renders smoothly, no triangle facets); P0 cells + where the field carries genuine cell-scale structure (yield-zone + boundaries, material interfaces — P1 blurs them by a node spacing; + measured side by side in + ~/+Simulations/fault_junction_rheology/p1-vs-p0.png). Recover + COMPONENTS and form invariants in numpy, never project an invariant + directly. Colour: linear RdBu_r at +-1 for dCFF; one-sided Oranges + on white for strain rate (maintainer convention — no dark + backgrounds, keep the mesh edges visible). (`signed_log` remains in + the harness but the maintainer rejected it for these figures.) +- Curved faults carry their ANALYTIC NORMAL + (`add_fault_bc(..., normal=sympy 1×dim Matrix in mesh.X)`). The kink + roughness of a sampled curve was diagnosed (2026-08-05, + ~/+Simulations/curved_fault_roughness/): NOT the meshing, NOT the + integration — the default per-node normal AVERAGES the adjacent + facet normals and zig-zags at the sampling kinks, so the no-opening + constraint forbids smooth slip past each kink (sawteeth that GROW + under refinement). The analytic normal on the same kinked mesh cuts + the sawtooth 7-17x and restores h-convergence. Straight faults need + nothing (the average is exact there); a deliberately kinked fault + should NOT be given a smooth normal — the kink response is then the + physics. + +## Conventions (consistent across every figure) + +- Mohr planes in the GEOLOGICAL sign convention: compression positive, + tension on the negative axis (solver tractions are tension-positive; + only plots flip). +- The interaction figures dress friction with a declared confining + pressure P0 = 1 and cohesion C = 0.75 (envelope + tau = +-(C + mu' sigma)); **neither enters Delta CFF** (constants + under differencing) — they place the failure line where teaching + needs it, and the pages say so. +- mu' = 0.4 (King's value) for all Delta CFF. +- Fields render in pyvista (RdBu_r, white background, lighting off); + line/scatter plots in matplotlib. Method-paper figures (cetz) use + sans labels (Noto Sans -> Helvetica fallback) at small sizes. + +## Regeneration + +Everything runs inside the worktree env, with its bin on PATH (the JIT +needs mpicc): + + cd /docs/advanced/figures/fault-examples + PATH=/.pixi/envs/runtime/bin:$PATH \ + /.pixi/envs/runtime/bin/python -u