Split-node faults: zero-thickness fault contacts in 2-D and 3-D, parallel, with interface constitutive laws - #502
Open
lmoresi wants to merge 64 commits into
Open
Split-node faults: zero-thickness fault contacts in 2-D and 3-D, parallel, with interface constitutive laws#502lmoresi wants to merge 64 commits into
lmoresi wants to merge 64 commits into
Conversation
…genuine discontinuity split_along_label grows the point chart the way rebuild_without_vertices compacts it: replica vertices for the chain interior, doubled fault facets (<name>Plus keeps the originals' labels, <name>Minus clones them via the clone map), cells kept in source order, edges re-derived so plus-side facets survive with their labels, coordinates repeated, and the star-forest renumbered by the same single broadcast. Tips stay unsplit, so a slip datum tapers to zero there. split_fault wraps the result as a standalone Mesh (no MG tail: the coarse levels do not carry the fault, so an essential condition there would leave the custom-P coarse operator singular) and records the Minus->Plus point pairing that any interface condition needs -- the sides are geometrically coincident, so no coordinate query can ever recover it. Refused loudly rather than mishandled: junctions, loops, single-facet chains, boundary-touching faults, and any fault whose cell fans touch the partition seam. The seam verdict outranks the chain-fragment symptoms each rank sees locally and is allgathered so every rank raises together. Serial tests: exact chart arithmetic (Euler characteristic 1 -> 0), conserved geometry with per-cell data alignment, DOF independence proven through a solve, label pairing into UW_Boundaries, five refusals, and re-application at a second fault position from the same base (the fault moves). Parallel (np=2/3): a rank-interior fault splits with zero star-forest coordinate drift; a seam-crossing fault is refused on every rank. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
Essential-BC DOFs are absent from the global vector, so the rotated path's per-field globalToLocal scatter left them at ZERO in the output fields wherever the Dirichlet datum was non-zero. The solve itself was always right -- the residual path inserts the values internally -- but every field-based diagnostic (projection, Integral, evaluate, renders) read a garbage boundary strip. Measured: far-field sigma_xy of 0.79 for a true 1.00, and a projected stress row off by 3x one fault length out. Invisible to every earlier rotated test because their Dirichlet walls were homogeneous: zero was accidentally the datum. Fix: petsc_dm_insert_boundary_values in cython/petsc_discretisation.pyx -- wrapping DMPlexInsertBoundaryValues, which petsc4py does not expose; the same insertion the consistent-boundary-flux paths already use (issues #407/#411) -- called per field after the scatter in _finalize_rotated_solution. The regression test drives an inhomogeneous lid with rotated free-slip sides and asserts the FIELD carries the wall values. The tell that found it: the divergence theorem, Integral(2 e_xy) against its own boundary term. The rotated boundary's prescribed datum and the sigma_nn reaction recovery were never affected (both live in the global vector). Fixes #497. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…inery For each coincident DOF pair of a split fault, build_rotation now writes an orthogonal 2*dim mean/jump block rotated into the fault (n,t) frame: mean rows on the Plus point, jump rows on the Minus point. The jump-normal row joins the strongly-constrained set (no opening, datum 0) and the slip row stays FREE, which is the zero-shear-traction condition -- its conjugate reaction is the shear traction and an unconstrained row carries none. Both points of a pair are rank-local by the split's seam refusal, so every block sits inside one rank's diagonal portion of Q exactly like the single-node wall blocks, and the whole rotated Newton loop, feasibility projection and line search carry over unchanged. A future datum on the slip row is jump-only prescribed slip; a nonlinear relation on it is friction. fault_contact.py holds the registration (against the pairing recorded by split_fault -- the sides are coincident, so the clone map is the only route to the pairs), the solve driver, and the fault_slip diagnostic that reads the jump through the pairing from the global vector. Validated against the analytic mode-II crack at the gate study's resolution: leak max |[v].n| = 5e-18; emergent slip elliptical to 1.35% RMS shape error with the residual confined to the under-resolved sqrt(r) tip zone; peak 92% of the infinite-medium value Delta-tau a / eta (finite box); Coulomb stress lobes in the same decay family as the inclusion and kinematic-slip representations, with amplitudes ordered exactly as their slips; converged in one Newton increment. Study and figures: ~/+Simulations/fault_split_gate/. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…esign add_viscous_fault_bc(solver, conds, boundary): tau = eta_f * V on the split fault, entering the rotated system as 2*eta_f*M on the slip rows (M the consistent 1-D P2 trace mass per fault facet; the tips drop out topologically because the jump space vanishes where the sides share a point, so the interface law needs no boundary condition of its own). The solve driver gains the two hooks every interface law uses: K.(Qu) added to the rotated residual, K added to a COPY of the ptap-refreshed operator (never injected -- the refresh owns its structure; the null space is re-attached on each rebuilt copy). Measured at the study resolution: the family bridges welded to free monotonically over three decades, following the crack compliance V/V_free = 1/(1 + 0.91 eta_f a / eta) with half-slip at eta_f = eta/a, and machine-zero opening at every member. The welded limit REMOVES the fault (recovers the uncut continuum) rather than stiffening it: only the jump is penalised, the mean velocity is never touched. A friction law replaces the constant with 2*(d tau/d V)*M in the same hook, and the measured compliance curve is its acceptance oracle at the secant eta_f = tau/V. docs/developer/design/FAULT_CONTACT_DEPLOYMENT_2026-08.md records the agreed deployment architecture: the three layers (persistent replicated manifold / ephemeral split mesh with the trace mapping / pair-transform constitutive layer), the two entry paths over one backend, and the parallel seam-vertex crossing rule (pin a vertex at each fault-seam crossing; one shared replica pair per crossing, keyed by root point and side). Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…-sector Three levels, each a legitimate model on its own. J0: disjoint segments with a ligament of one or two local h -- available today under the existing refusals, validated by the King two-fault interaction pattern; a crossing is two offset abutments. J1: true tip-on-fault abutment -- split the through-going fault first, then the abutter with its unsplit tip on the master's slit; the only code change is the boundary-touching refusal distinguishing, for tips only, a domain boundary from a prior fault's slit. J2: the exact degree-d sector split, where branch compatibility is a telescoping identity of sector differences (no cycle constraint, no multiplier) and the acceptance test is machine-precision closure of the branch slips at the junction. Junctions stay off partition seams at every level. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…plain solve() The user-facing layer over the split-node backend, per the deployment design. Mesh.add_fault(faults) runs the whole pipeline in one call — tips placed onto vertices, conforming cut, split with the coincident pairing recorded — accepting a Surface, a (name, points) pair, or a sequence of either. A sequence is a NETWORK: every fault is cut first, then every fault is split, which makes the offset-junction (J0) pattern a one-liner. Splitting renumbers the whole chart, so prior faults' pairings are carried through the new split's point_map rather than copied verbatim — verbatim ids silently index the wrong points, which is exactly how the network case failed first. solver.add_fault_bc(conds, boundary) is the value-first BC: conds = 0 the frictionless contact, conds > 0 the viscous interface tau = eta_f V. With a fault registered, an ORDINARY solve() dispatches to the rotated fault-contact path; guard() and estimate_difficulty() refuse loudly there, exactly as for rotated free-slip. End-to-end tests: one-call fault + plain solve reproducing the crack behaviour, the viscous law through the solver method, Surface-object input, and a two-segment offset network slipping on both faults with machine-zero opening. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…named The blanket refusal (any fault-fan point shared) was far stricter than the invariant requires. What the rank-local split actually needs is that every DUPLICATED vertex owns its whole cell fan — and an UNSHARED vertex always does (any cell touching it from another rank would make it shared), which also keeps every re-homed spoke rank-local. So the rules become: a support-1 facet in a fault vertex's star is the domain boundary only if it is UNSHARED (on a seam every facet looks one-sided locally); and the split refuses only a CHAIN VERTEX on the seam — a seam crossing — with a message naming the crossing milestone. Faults may now run arbitrarily close to a partition seam. np=2: a fault reaching toward the box centre (refused by the old rule) splits cleanly with zero star-forest drift; np=3 lands a chain vertex on a seam and gets the collective crossing refusal. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
The seam-vertex crossing rule of the deployment design, end to end. A crossing is a shared chain ENDPOINT whose global fault-facet count is 2 (an SF-sum of owned facets — shared facets would double-count); it gets a replica via an OPEN-fan walk (anticlockwise from the one local fault neighbour to the fan's dead end sweeps exactly the left of the directed chain, with the other rank's fault direction inside the gap), under a globally consistent chain orientation passed down from the fault polyline — fragments cannot use the lexicographic-tip rule, each rank would pick its own. The star-forest rebuild gains a second broadcast carrying each point's replica/replacement index, so a leaf rank associates ITS copy of a duplicated shared point with the owner's: one keyed entry per crossing vertex, plus the re-homed seam spokes, with hard errors on any one-sided decision. The existing compacted-chart path is untouched (reconnect parallel suite stays green). Mesh.add_fault manufactures the clean crossing the rule requires: an organic cut STRADDLES the seam with a facet held by both ranks (refused, with its own message), so add_fault intersects the polyline with the seam edges and pulls a SHARED vertex onto each crossing point — preferring vertices shared by exactly two ranks (corners are refused downstream), capped at the local cell scale so the pull cannot invert cells that the cut would then reject. Refused, collectively, with named messages: junctions on the seam, crossings on three-rank corners, faults touching the seam without crossing, and crossing splits without a global orientation. Interface trace-mass assembly now allows off-rank rows (a crossing pair is owned across the seam; PETSc communicates stashed entries). np=2: all four parallel tests including a frictionless solve slipping across the ranks with machine-zero opening; np=4: crossing passes; np=3: this partition's geometry draws the honest collective refusals. Serial suites and the reconnect parallel suite unchanged. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
The np=3/4/5 robustness sweep (crossing_sweep.py in the study directory) found one genuine rule: EVERY polyline control point needs a vertex pulled onto it — an interior kink is the same problem as a tip, and without its vertex the cut leaves the chain fragmented at the turn. add_fault now pulls all control points. Sweep results: 17/27 geometries split cleanly across np=3/4/5 with every success validated (zero star-forest drift, global Euler 0, conformity, machine-zero solve leak) and every refusal collective and categorised (three-rank corners; along-seam faults). No hangs, no invalid meshes. The deployment design gains the representation policy: the fault object is the physical entity and each equation system chooses its adapter — surface contact (this branch), volume weak/damage zone or TI director (the existing meshing/faults.py machinery on the same manifold). Default composition is different meshes per continuity class (Stokes on the split child, Darcy/thermal on the continuous parent); a shared split mesh needs a scalar weld constraint, deferred until tight coupling demands it. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
…angent The interface hook graduates from a constant operator to a law layer that keeps sympy working for the fault exactly as it does for the bulk. A law is a sympy expression in fault_contact.slip_rate; its consistent tangent d tau/dV comes from sympy.diff — never hand-coded — and both are lambdified once at registration (the interface analogue of the JIT step). The _InterfaceAssembler caches the fault-trace geometry and, at EVERY Newton iterate, adds the interface force integral(tau(V) dV) to the rotated residual and the consistent tangent 2 (dtau/dV) M to the rotated operator: full Newton in V, the stiff direction. There is no Picard limitation from the hand-rolled coupling — the manual loop makes the tangent exactly dF/du; only a future reaction-fed sigma_n argument would be lagged, by choice. Tip nodes enter the quadrature with V = 0 and carry no DOF; a crossing pair's off-rank rows go through the PETSc stash. ViscousFaultLaw becomes the linear member of the same symbolic family (the welded-to-free suite passes unchanged through the general assembler); CoulombFaultLaw is the first nonlinear one, tau = mu sigma_n (2/pi) arctan(V/V0) with prescribed sigma_n. Tested against both analytic regimes on one mesh: strength below the resolved driving shear slides at the constant-stress-drop scaling of the frictionless crack; strength above it sticks, slip collapsing to the regularisation creep ~V0 — converged Newton in both. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
sigma_n = "reaction" feeds a Coulomb law from the no-opening constraint's own reaction instead of a prescribed value. The recovery is the pair analogue of the wall sigma_nn machinery: the Cartesian residual stashed at each iteration START is the pure transmitted traction (the interface force is added after the stash), so the nodal normal load is (n.r+ - n.r-)/2, de-smeared by the lumped 1-D P2 trace mass (positive row sums, unlike P2 triangles; completed over the star-forest so a seam-crossing node sees both ranks' facets). Compression enters the law positive; tension clamps the strength to zero — an opening-tending fault has no friction, and the no-opening constraint holds it shut regardless. The lag is placed INSIDE rotated_residual, between the reaction stash and the interface force, gated on keep_cartesian: sigma updates exactly once per Newton iteration, and line-search trials evaluate against the frozen sigma — a consistent comparison, the Picard-lag done properly. Iteration zero sees a rest state, so a reaction-fed law starts frictionless and tightens as the stress field forms: a built-in homotopy. Laws become two-argument symbolic expressions (slip_rate, normal_stress), tangents still by sympy.diff. New public diagnostic fault_normal_traction returns signed sigma_nn along the fault from the same geometry. Tested: the recovered mean matches the resolved background compression 2 n_x n_y on the frictionless fault, and the reaction-fed Coulomb slides at the constant-stress-drop ratio its recovered strength predicts. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
RateStateFaultLaw is the arcsinh-regularised form (never the raw logarithmic law, singular at V = 0): tau = a sigma_n asinh[(V/2V0) exp((f0 + b ln(V0 theta/Dc))/a)] — odd in V, smooth, with the state entering per node through the new state_variable symbol exactly as the normal stress does, and the consistent tangent in V still from sympy.diff. At steady state it reduces to sigma(f0 + (a-b) ln V/V0). The state lives per fault on the solver (by plus-point id — stable across assembler rebuilds; manifold MeshVariables replace this store when the 1-D mesh type lands, with the law layer unchanged) and advances between solves by update_fault_state: the ageing law integrated EXACTLY for piecewise-constant V, written with expm1 so V -> 0 degrades to pure ageing. The function returns the steady-state monitor median theta*V/Dc. One numerical lesson: theta varies by orders of magnitude along a fault, and quadratic shape functions undershoot — plain interpolation sent a positive-only state negative at a quadrature point (fractional powers of it, from b/a exponents, then poisoned the law). Theta is interpolated in the LOG, which is both positive-preserving and the natural interpolant for a variable the law only reads through ln(theta). Tested: a velocity-strengthening fault under steady loading cycles solve <-> state-update to the fixed point — slip settles and the monitor reaches Dc/V — with sigma_n reaction-fed throughout. Underworld development team with AI support from Claude Code Claude-Session: https://claude.ai/code/session_01WDKP73LD3dCFQSNgEXbbw9
… rule The 3-D counterpart of split_along_label: the fault is a triangulated patch of interior faces; everything strictly inside its rim — vertices, edges, and the faces themselves — is doubled, and the rim stays unsplit (slip must taper to zero there, exactly as at the 2-D tips). Where the 2-D rebuild hand-wires edge cones, the 3-D one builds the substituted cells UNINTERPOLATED (vertex 4-tuples carry no orientation) and lets DMPlexInterpolate derive faces and edges: cell and vertex numbering are preserved (asserted), and every old face or edge is recovered in the new chart by joining its mapped vertex tuple — a tuple that no longer joins was re-homed onto replicas, the 2-D fresh-facet rule one level deeper. Sides come from an orientability propagation over the patch plus a signed-volume test per fault face, and each interior vertex's cell star must flood into exactly two half-balls. Refusals: non-manifold, pinched, non-orientable or disconnected patches; a patch touching the domain boundary; a patch so coarse that a face has no interior vertex (its two copies would be one triple); and ANY partition-seam contact — a 3-D seam crossing is a curve, not the 2-D crossing vertex, and is deferred by design. uw.meshing.BoxInternalPatch supplies the conforming geometry for the science-first path: a planar polygon patch gmsh-embedded in a simplex box, physical name = fault name, rim strictly interior. split_fault dispatches on dimension; the pairing dict now carries vertex, edge AND face twins (P2 needs the edge midpoints). Tests (test_0848): exact chart deltas against an independent patch census; the Euler invariant chi 1 -> 2 (the slit adds the patch interior's characteristic); conformity and one-sided side labels; volume conservation and bit-identical replica coordinates; the load-bearing P2 DOF-independence solve (Plus datum must not leak onto coincident Minus DOFs, rim reads it from both sides); refusals. Underworld development team with AI support from Claude Code
The pair blocks become 2x3 mean/jump in the (n, t1, t2) frame — the jump-normal row constrained as always, TWO free slip rows — with the frame drawn from a single authority (_tangent_frame) shared between rotated_bc's rotation and the interface assembler's slip components. _InterfaceAssembler generalises to (dim-1)-component slip. The 3-D trace element is the P2 triangle (6-point degree-4 Dunavant rule); laws act on the slip magnitude with collinear traction tau_vec = tau(|V|) V-hat, and the consistent tangent is the 2x2 block (dtau/dV) VV^T + (tau/|V|)(I - VV^T), regularised at |V| -> 0 by the laws' own smoothness (the zero-slip block is dtau/dV(0) I). The sigma_n de-smear uses the P1 SUB-LUMPED triangle mass (A/12 per vertex, A/4 per midpoint) because the lumped P2-triangle vertex rows vanish — the known trap. The 2-D path keeps its signed-V evaluation and lumped line mass; the full 2-D suite is unregressed. Validation: penny-shaped frictionless crack against the classical oracle dv(r) = (8/3pi)(dtau/eta) sqrt(a^2 - r^2) (nu = 1/2) — leak at machine zero, slip parallel to the resolved shear, amplitude below the infinite-medium value and converging toward it from below under both h-refinement (peak ratio 0.930 -> 0.947) and box growth (-> 0.941), as the finite Dirichlet box demands. Study in ~/+Simulations/fault_split_3d/penny_crack. New level_2 tests: the viscous family (monotone to the weld) and Coulomb slide/stick with reaction-fed sigma_n, both in 3-D. Also from liaison with the placed-surface session: eta_f = eta_layer/w is the LAYER'S weak-zone viscosity, not the background's (both docstrings reworded), and the manifold contract is recorded as a protocol with Surface (2-D) and FaultSurface (3-D) as concrete types. Underworld development team with AI support from Claude Code
A bare technical outline for publication drafting: formulation (split space, pair transform, symbolic interface laws, consistent tangent, reaction-fed normal stress, exact ageing), the mesh operation in 2-D and 3-D with the chart arithmetic and refusal table, the trace assembly (including the vanishing lumped P2-triangle mass and its P1 sub-lumping), and the validation numbers: the 10-vs-147 conditioning gate, the 2-D crack/compliance/friction/rate-state results, and the 3-D penny-crack convergence table. Figures follow the cetz-figures pattern (generator + JSON + typ + committed PNG) in docs/developer/design/figures/split-node-faults/: the 2-D split anatomy, the pair-transform row diagram, and a pyvista render of the doubled patch and shared rim from an actual split mesh. Underworld development team with AI support from Claude Code
Five panels from deterministic Delaunay point sets: the static base with the fault as a geometric manifold, two adapt-on-top refinement layers graded toward it, the cut that turns it into a conforming facet chain, and the split (exploded), with the non-cumulative re-derive loop annotated. The chain's conformity in the generated data is verified programmatically, not assumed. Underworld development team with AI support from Claude Code
The band points were laid out in the fault's own frame, so layer 2 manufactured the conforming line before the cut had anything to do — and clearing the corridor around it left visibly lopsided flanks. Refinement points now come from an axis-aligned lattice with deterministic jitter, restricted to the distance band: the layers refine TOWARD the manifold without containing it, and the conforming chain in the cut panel is unambiguously the cut's contribution. Underworld development team with AI support from Claude Code
…hout The stack-progression's split panel now shows the mesh in its TRUE (unexploded) geometry — the copies are coincident — with the doubled surface drawn as the conventional double line: two thin parallel strokes offset along the fault normal, white between, tip rings at the shared ends. All three cetz figures move to finer line weights and compile at 300 ppi. Underworld development team with AI support from Claude Code
The vertical companion to the stack progression: the static base is the finest level of a geometric (FMG) hierarchy — two coarser levels drawn below it — and the adapted children extend that hierarchy upward, one level per doubling of resolution. The cut + split working mesh sits on top of the stack but outside the hierarchy, separated by the dashed rule "a cut is not a multigrid level". Typography across the whole set: labels move to a sans stack (Noto Sans, Helvetica, Arial — the first available wins) at half the previous size; math stays in the math font. Underworld development team with AI support from Claude Code
The coarse levels are now genuine structured refinements of one another — the same diagonal-split generator at h = 1.0, 0.5, 0.25 on a 2.0 x 1.0 domain, so midpoint refinement of each level IS the next, and the generator VERIFIES the nesting (coarse vertices a subset of fine) rather than assuming it. The single FMG bracket splits in two: the grey bracket spans the exactly-nested structured levels (lossless transfers, perfect FMG); an amber bracket and amber "adapt, non-nested" arrows mark the adapted children, where transfers become approximate. The stack-progression panels move to the same domain so the two figures share one base mesh. Underworld development team with AI support from Claude Code
The parallel sweep caught a mis-diagnosis: a rank holding only a ghost FRAGMENT of the patch sees the fragment's edge as rim, its faces as support-1 and its interior as empty, so the fragment-sensitive checks (boundary contact, too-coarse) fired with misleading messages where the true verdict was the seam. The seam rule now runs FIRST — its patch-vertex union is fragment-correct even though the rim/interior classification is not — and the topology checks only run on a rank whose patch is seam-free. ptest_0848: at np = 2 an elongated box pins the balance cut across the long axis (probed), so a patch near one end is genuinely rank-interior: it splits, the star forest rebuilds with zero coordinate drift, and a P2 solve shows the Plus datum decoupled from the coincident Minus DOFs. A patch on the cut refuses with the seam RuntimeError on every rank. np = 3 and 4 partition through the patch — a legal layout — and now refuse with the correct collective seam verdict. Underworld development team with AI support from Claude Code
The user-facing page for the capability: quick starts in 2-D (mesh.add_fault -> add_fault_bc -> solve -> fault_slip) and 3-D (BoxInternalPatch -> split_fault), the four fault laws with their call signatures and the physics of each parameter, prescribed kinematic slip with the tip-taper rule, and the practical notes (P2/P0-disc pressure, tolerance not ksp_rtol, moving faults re-derive from the base, the parallel contract). Cross-linked with the TI weak-zone page as the two complementary fault representations, and with the warning that fault quantities must be read through the DOF pairing — coordinate queries see one side of a coincident pair. Linked into the advanced toctree; sphinx build verified clean for the new page. Underworld development team with AI support from Claude Code
…ation vs slip Three short scripts on one shared harness (base mesh re-faulted per case — the non-cumulative add_fault pattern), each producing a committed figure, hosted on a new docs/advanced page: - The LADDER: one fault, one shear drive, five solves down the constitutive ladder. Measured: frictionless peak 0.174, the eta/a dashpot at ~half, Coulomb-weak at exactly the reduced stress drop (1 - mu sigma_n / tau) x free, rate-state between, Coulomb-strong stuck at the V0 creep floor. - The MOHR CIRCLE: welded faults as passive stress probes — sigma_n from the constraint reaction, SIGNED tau from the dashpot's own law — swept through 180 degrees of orientation. Eight probes trace the full circle; fitted radius 1.411 vs the analytic sqrt(2), centre at the zero-pressure gauge. The double-angle rule is the label spacing. - ORIENTATION vs SLIP: the same sweep, frictionless: peak slip follows |cos 2 theta| to three digits, with an EXACT zero at 45 degrees. Underworld development team with AI support from Claude Code
The ladder's two fixes: slip positions are now read EXACTLY through the DOF pairing (fault_slip's arc-length origin is the first pair, not the tip — the profiles were offset by half a node and stopped short of the fault ends) with the zero-slip tips appended, so the profiles span the same fault the elliptical shape is drawn for — the frictionless curve now sits ON the ellipse. And each rung gains a stacked shear-stress panel (sigma_xy to continuous P1, RdBu_r about the far-field value): the stress-drop shadow and tip concentrations fade down the ladder to a featureless field for the stuck fault. The Mohr circle gains its teaching form: an animated build (mohr-circle-build.gif, 25 welded-fault solves, one per frame) with the fault physically rotating in the left panel while its measured stress probe sweeps the circle at 2*theta on the right — the double-angle rule as motion. The static figure stays for print; the probe loop is hoisted into common.mohr_probe as the single authority. Underworld development team with AI support from Claude Code
Dirichlet velocity on all four walls imposing the homogeneous flow v = (a(x-c) + gamma(y-c), -a(y-c)): pure shear (stretching, a) plus simple shear (gamma). Stress sees only the symmetric gradient, so the state is equivalent to pure shear of magnitude eta sqrt(4a^2+gamma^2) with principal axes at 22.5 degrees — deliberately not axis-aligned. The linear flow is an exact homogeneous Stokes solution, so the stress is uniform to the walls and every probe samples the same state; the pressure gauge sets the circle centre. Underworld development team with AI support from Claude Code
The left panel now draws the MEASURED traction vector sigma.n at the fault centre: it swings with the rotation, lies along the fault plane near the maximum-shear orientations, and snaps (anti-)parallel to the normal exactly at the principal orientations — at the same instant the probe crosses tau = 0 on the right, where the two principal stresses are now marked. A caption flashes at the principal frames. The 25-solve probe sweep is cached (_mohr_probes.npz, untracked) so animation cosmetics no longer re-measure the physics. Underworld development team with AI support from Claude Code
…elope The rotating fault now carries Coulomb friction (mu = 0.7, reaction-fed sigma_n). Three regimes appear in one measured sweep, exactly as the theory demands: stuck probes reproduce the welded circle values to four digits (the ambient stress lies inside the envelope); sliding probes are pinned to the yield rays tau = +-mu|sigma_n| (measured tau = 0.491 at sigma_n = -0.702 against mu|sigma_n| = 0.491); the tensile sector carries no strength and collapses to the sigma axis. Both principal orientations end up not slipping — the compressive pole because it is inside the envelope, the tensile pole because zero resolved shear leaves nothing to drive. tau 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. Static summary + animated build (slip sense as half-arrows when the fault unlocks); probe sweep cached like the welded one. Underworld development team with AI support from Claude Code
Cohesive Mohr-Coulomb (C = 0.6, mu = 0.6), registered as a sympy law in the canonical interface symbols — the intended extension path, four lines. Measured: stuck arcs survive around BOTH principal poles (compressive because C + mu sigma = 1.2 exceeds the ambient shear; tensile because C = 0.6 does), sliding probes pin to the rising Coulomb rays on the compressive side (tau = 1.0215 vs C + mu sigma = 1.023) and to the FLAT cohesion lines under tension (tau = 0.5994 vs C = 0.6) — cohesion under tension is shear strength only, the no-opening constraint already excludes opening. Static + animated. All four Mohr figures now use the geological sign convention — compression positive, tension on the negative axis — with the flip applied at plot time only (solver tractions stay tension-positive, and the animations' physical-space traction vectors are untouched). Regenerated from the cached probe sweeps. Underworld development team with AI support from Claude Code
… own The assembler clamped the reaction-fed normal stress at zero before the laws saw it — so a cohesive law could not know the fault was in tension and its strength sat flat at C forever, a physics defect caught by the maintainer on the Mohr figure. Where strength vanishes is CONSTITUTIVE: bare friction loses it at sigma = 0, cohesion at sigma = -C/mu. update_normal_stress now feeds the signed stress and each law clamps itself (sympy Max; the kink never enters the Newton direction because sigma is Picard-lagged). Canned Coulomb and rate-state laws gain Max(sigma, 0) — behaviour unchanged in compression, all 15 law tests pass. The Mohr examples now tell the truth about tension: the cohesive envelope declines through mild tension and reaches zero at sigma = -C/mu (measured: tau = 0.375 at sigma = -0.375, 0.172 at -0.714 — C + mu sigma to three digits); beyond the apex, and everywhere tensile for bare friction, the fault WOULD OPEN — no static solution exists, and the bilateral no-opening constraint manufactures one by carrying a tensile reaction. Those probes are now marked HELD SHUT (unphysical) in a shaded no-solution region, and the docs state the diagnostic: the sign of the recovered normal traction. Also answered in the docs: there is no confining pressure in these velocity-driven incompressible problems — the pressure gauge centres the circle on the origin, which is why tension is reachable. Underworld development team with AI support from Claude Code
A modest gravity load (rho g = 0.75, constant density, closed box — the flow is untouched, pressure absorbs the body force exactly) makes the per-node stress recovery visible on the Mohr diagram: each welded fault contributes a horizontal STREAK of depth-coloured probes rather than a single value, spanning the family of Mohr circles between its shallowest and deepest points. Streak length follows the fault's vertical extent exactly (a dot for the horizontal fault, rho g dy for the vertical one), and nothing is averaged: sigma_n per node from the constraint reaction, tau per node from the weld's own law. The closed-box pressure gauge is re-anchored in the plot so p = 0 at the surface (the analytically known rho g H/2 shift, stated openly). This is the machinery that lets a friction law feel depth-dependent strength along one fault — the locking-depth structure. Underworld development team with AI support from Claude Code
Member
Author
Adversarial reviewWe reviewed this branch as adversaries before asking anyone else to. Six 1. Nodal 2. The held-shut regime has no runtime tell. A zero-strength fault 3. Welded-probe slip shows a small rank dependence. In the 4. The redistribution target keeps its full far-field share. The 5. The unwrap coordinate-retag fix exists only on the fault path. 6. Diff size. 90 files, +16.6k lines, of which roughly a third is 7. The Surface idiom is plumbed but not practised, and crossings are a None of these are correctness defects in what the PR claims; 1 and 2 are Underworld development team with AI support from Claude Code |
Contributor
There was a problem hiding this comment.
Pull request overview
This PR introduces a user-facing split-node fault capability (2-D and 3-D, including parallel workflows) by extending the rotated strong-BC machinery to support coincident DOF-pair contact constraints and interface constitutive laws, alongside supporting meshing, visualisation, tests, and documentation.
Changes:
- Extend rotated strong-BC solving to support split-fault contact pairs and interface-law Jacobian contributions, plus correct field copy-back for inhomogeneous essential BCs.
- Add 3-D meshing support for an internal planar patch (
BoxInternalPatch) and expand visualisation helpers for hierarchies/fault facets and native P1 plotting. - Add extensive serial/parallel tests and developer/user documentation (including figures and teaching example scripts).
Reviewed changes
Copilot reviewed 67 out of 90 changed files in this pull request and generated 1 comment.
Show a summary per file
| File | Description |
|---|---|
| tests/test_1018_rotated_freeslip.py | Adds regression test for rotated-path field copy-back with inhomogeneous Dirichlet walls. |
| tests/test_0847_mesh_hierarchy_plot.py | Adds rendering-behaviour tests for plot_mesh_hierarchy (actors, faults, glyphs, legend). |
| tests/test_0847_fault_api.py | Adds end-to-end tests for the user fault API (add_fault, add_fault_bc, normals). |
| tests/test_0846_visualisation_native_mesh.py | Tests that continuous P1 fields render on native mesh connectivity (not Delaunay). |
| tests/test_0845_relax_pinned_band.py | Adds serial tests for relax(pin_bands=...) behaviour and boundary pinning. |
| tests/test_0843_edge_split_adapt.py | Adds serial tests for engine="edge_split" refinement properties and confluence contract. |
| tests/test_0840_nvb_3d_serial_adapt.py | Refactors/strengthens 3-D adapt MG-level assertions using shared ladder check. |
| tests/test_0836_nvb_graded_adapt.py | Refactors MG-level assertions and adds tests for mg_coarsening_ratio effect. |
| tests/parallel/ptest_0848_fault_split_3d_parallel.py | Adds 3-D parallel fault split tests including redistribution and SF drift checks. |
| tests/parallel/ptest_0845_relax_pinned_band_parallel.py | Adds parallel tests for partition-independent pin sets and pinned-vertex immobility. |
| tests/parallel/ptest_0845_fault_split_parallel.py | Adds parallel tests for split-node faults (collective refusal, seam/crossing cases). |
| tests/parallel/ptest_0843_edge_split_parallel.py | Adds parallel confluence tests for engine="edge_split" and adapt-tail checks. |
| tests/_mg_ladder.py | Introduces shared MG coarsening “ladder” assertion helper for adapt children. |
| src/underworld3/visualisation/init.py | Exposes new/updated visualisation helpers and constants in the public API. |
| src/underworld3/utilities/rotated_bc.py | Extends rotation assembly for fault pair blocks; adds interface assembly/tangent; fixes copy-back gap. |
| src/underworld3/utilities/nvb.py | Adds explicit TODO note about known 3-D prolongation inaccuracy. |
| src/underworld3/utilities/init.py | Exports new utilities modules (edge_split / fault_contact / fault_split / etc.). |
| src/underworld3/meshing/cartesian.py | Adds BoxInternalPatch mesh factory for embedded internal planar patches in 3-D. |
| src/underworld3/meshing/init.py | Re-exports BoxInternalPatch in meshing public namespace. |
| src/underworld3/cython/petsc_generic_snes_solvers.pyx | Adds add_fault_bc and routes fault-contact solves through rotated solve path; updates guard probes. |
| src/underworld3/cython/petsc_discretisation.pyx | Adds Cython shim petsc_dm_insert_boundary_values wrapping DMPlexInsertBoundaryValues. |
| docs/developer/subsystems/conforming-surfaces-and-fault-zones.md | New subsystem doc explaining conforming surfaces, zones, and implications. |
| docs/developer/index.md | Adds new subsystem doc to developer documentation index. |
| docs/developer/design/figures/split-node-faults/stack-progression.typ | Adds Typst/cetz source for split-node fault method figure. |
| docs/developer/design/figures/split-node-faults/split-anatomy.typ | Adds Typst/cetz source for split anatomy figure. |
| docs/developer/design/figures/split-node-faults/split-anatomy-data.json | Adds generated geometry data for split anatomy figure. |
| docs/developer/design/figures/split-node-faults/render-components-3d.py | Adds PyVista script to render 3-D split components figure from real mesh. |
| docs/developer/design/figures/split-node-faults/pair-transform.typ | Adds Typst/cetz source for pair-transform schematic. |
| docs/developer/design/figures/split-node-faults/grid-hierarchy.typ | Adds Typst/cetz source for grid hierarchy schematic. |
| docs/developer/design/figures/split-node-faults/grid-hierarchy-data.json | Adds generated geometry data for grid hierarchy schematic. |
| docs/developer/design/figures/split-node-faults/generate-stack-progression-data.py | Adds generator for stack progression figure geometry. |
| docs/developer/design/figures/split-node-faults/generate-split-anatomy-data.py | Adds generator for split anatomy figure geometry. |
| docs/developer/design/figures/split-node-faults/generate-grid-hierarchy-data.py | Adds generator for grid hierarchy figure geometry. |
| docs/developer/design/FAULT_CONTACT_DEPLOYMENT_2026-08.md | Adds deployment/architecture design note for split-node fault contacts. |
| docs/developer/ai-notes/fault-teaching-examples-handoff-2026-08.md | Adds handoff note for teaching examples and figure regeneration workflow. |
| docs/advanced/index.md | Adds new advanced-doc entries for split-node faults and teaching examples. |
| docs/advanced/figures/fault-examples/orientations.py | Adds teaching script: orientation sweep vs slip behaviour. |
| docs/advanced/figures/fault-examples/mohr_graded.py | Adds teaching script: hydrostatic loading sampled along welded faults. |
| docs/advanced/figures/fault-examples/mohr_friction.py | Adds teaching script + animation: Coulomb friction envelope vs Mohr circle. |
| docs/advanced/figures/fault-examples/mohr_cohesion.py | Adds teaching script + animation: cohesive Mohr-Coulomb behaviour. |
| docs/advanced/figures/fault-examples/mohr_circle.py | Adds teaching script: Mohr circle measured by welded faults. |
| docs/advanced/figures/fault-examples/mohr_animate.py | Adds teaching script + animation: Mohr circle build as fault rotates. |
| docs/advanced/figures/fault-examples/ladder.py | Adds teaching script: “fault-strength ladder” comparing multiple laws. |
| docs/advanced/figures/fault-examples/.gitignore | Ignores generated frame/log outputs for teaching example renders. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+451
to
+457
| from underworld3.cython.petsc_discretisation import \ | ||
| petsc_dm_insert_boundary_values | ||
| for name, var in solver.fields.items(): | ||
| sg = U.getSubVector(solver._subdict[name][0]) | ||
| solver._subdict[name][1].globalToLocal(sg, var.vec) | ||
| U.restoreSubVector(solver._subdict[name][0], sg) | ||
| petsc_dm_insert_boundary_values(solver._subdict[name][1], var.vec) |
…onverted The two pre-release follow-ups from PR #502, delivered: - The 3-D pipeline now consumes uw.meshing.FaultSurface end to end: FaultSurface.rim_polygon() extracts the ordered rim of a planar triangulated surface (boundary-edge chain; non-planar refuses toward the discrete-entity embed); BoxInternalPatch(patch_points=fault_surface) takes the object, adopts its name and stores it on the mesh; split_fault carries it (and the 2-D stored traces) onto the split child; add_fault_bc(..., normal="surface") builds the constraint frame from the surface's OWN face normals (nearest-face lookup — exact for a planar patch, unchanged for future curved sheets). - uw.meshing.prepare_fault_network converts an imported 2-D trace set into splittable offset-junction form, loudly: X crossings cut both traces, T abutments trim the abutting end, near-miss endpoints are pulled to clearance. Pull-backs are angle-corrected (ligament / sin theta — the Euclidean-clearance check caught the along-trace version being short at oblique junctions) and pieces shorter than two ligaments are dropped and reported. - Ligament sensitivity measured (T-junction, 35 degrees, both faults frictionless, ~/+Simulations/fault_junction_ligament/): branch peak slip is 1.01x / 0.96x / 0.79x its isolated value at ligament = 1h / 2h / 4h, and at FIXED ligament 0.04 the whole branch slip profile is IDENTICAL at h = 0.02 and h = 0.01 — the offset junction is a converged physical answer for a given plug size, not a resolution artifact. Keep ligaments at 1-2 local cell sizes; a true shared-vertex junction design is an accuracy refinement, not a correctness necessity. Tests: test_0849_fault_network_prep (X/T/near-miss conversion, Euclidean clearance, end-to-end add_fault) and test_0848::test_fault_surface_route (rim extraction, name adoption, surface-normal compile, non-planar refusal). User guide updated to lead 3-D with the FaultSurface idiom. Underworld development team with AI support from Claude Code
Member
Author
|
Two of the pre-release follow-ups have now landed on the branch (and in this PR): the FaultSurface object is the 3-D input idiom end to end (rim embed, stored on the mesh, The ligament approximation is now measured, not asserted: at a 35-degree T-junction with both faults slipping, the branch's peak slip is 1.01x / 0.96x / 0.79x its isolated value for ligaments of 1 / 2 / 4 local cell sizes, and at fixed ligament 0.04 the entire branch slip profile is identical at h = 0.02 and h = 0.01. The offset junction is a converged physical answer for a given plug size; keep ligaments at 1-2 local cells. A true shared-vertex junction design is therefore an accuracy refinement, staged, not a correctness gap. Underworld development team with AI support from Claude Code |
Mesh.add_fault now stores Surface OBJECTS on the split mesh (alongside the raw traces), and normal="surface" resolves in 2-D to the smoothed normal of the Surface's own control polyline — the same spelling and the same object-first idiom as the 3-D FaultSurface route. Deliberately NOT the Surface's director (signed-distance gradient): that field is discontinuous exactly ON the trace, where the fault nodes live; the director stays the right tool for TI weak zones, the polyline normal for the contact frame. Test extended with the frame-equality check. Underworld development team with AI support from Claude Code
The prepare_fault_network capability shown end to end: raw traces that genuinely intersect (a dextral trunk, a splay abutting its midpoint at 30 degrees, a conjugate crossing outright) are auto-converted to offset junctions and split in one add_fault call; the trunk and splay rupture together, the conjugate is welded. The Delta CFF map's brightest features are the junction plugs themselves — the intact ligaments take the concentrated load shed by their slipped neighbours, which is where a through-going event would break next — shown at map scale and in ligament-scale zooms at the T and X junctions. Page section added with the measured ligament-convergence statement. Underworld development team with AI support from Claude Code
prepare_fault_network was cutting the through-going trace at T abutments, contradicting its own documentation (the trunk in the branching example was interrupted twice). Now: a T abutment trims only the abutting end; X crossings cut both traces by default; and through=[names] declares MASTER faults that stay continuous even at crossings (the other trace yields on both sides; two masters crossing is a hard error). Tests updated to the corrected contract. The cost of interrupting a fault is now measured, not guessed (branching_compare.py, on the teaching page): with the trunk through-going its peak slip is 0.2590 — slightly MORE than the isolated trunk's 0.2558, because the slipping splay feeds it; cut at the crossing it drops to 0.1963 (about a quarter forfeited) and is pinned at the junction plug, but the profiles converge away from it — an interruption's reach is the segment scale, not the system scale. The both-cut state the old defect produced was the worst of all (0.144). branching.png regenerated under the corrected default so the committed figure matches the code's behaviour. Underworld development team with AI support from Claude Code
Does a near-miss tributary reproduce a genuine branch? Two findings, one per axis. The GAP is second-order: with the trunk continuous and the splay abutting, ligaments of 1h/2h/4h give the same slip on every arm except the splay's near-junction toe, which moves with the trimmed tip. The DECOMPOSITION is first-order: routing continuity through the bend instead (west arm + splay as one kinked fault) locks at the 33-degree kink — slip through a kink is geometrically incompatible with no-opening, the curved-fault mechanism appearing as genuine physics — costing both trunk arms a third of their slip. The rule for imported networks: declare the straightest path continuous; never route continuity around a corner. The continuous-trunk decomposition is admissible for the true junction and differs from it only in the ligament-scale toe, which is all a true shared-vertex junction design would buy. Underworld development team with AI support from Claude Code
Delta CFF maps for both decompositions of the Y-branch, map scale plus branch-point zooms at one colour scale. Trunk-continuous: the stress shadow runs straight through the junction, which carries no feature at all — mechanically invisible, as the declared through-going pair should be. Bend-continuous: the locked kink is a barrier — both arms' shadows terminate against it and dump load into the corner, the rupture-arresting-bend signature. Noted in the log: decomposition B's welded reference landed 139 stress units away in pressure gauge from its slipping solve — the far-field anchor removed it (printed), the strongest case yet for the anchoring discipline. Underworld development team with AI support from Claude Code
Three chains terminating at one shared vertex are refused when the second split's tip lands on the first fault's slit — but the message said "touches the domain boundary", which is what the predicate sees (one-sided facets in the tip's star), not what the user did. The refusal now detects that the offending facet carries a Plus/Minus side label and says so: a shared point would clamp every arm's slip to zero there — the kink-lock in its purest form, STIFFER than a true sector junction — and points at prepare_fault_network, whose offset form brackets the true branch (the true-branch teaching example). Test added; deliberately NOT relaxed to allow the touching case, since its mechanics is reproduced to ligament accuracy by the continuous-bend decomposition already measured. Underworld development team with AI support from Claude Code
Member
Author
Adversarial review: PR #502 — split-node faults (fault delta over #488)Scope: the 53-commit fault delta VerdictThe core mechanism is sound and well-tested: the split is a clean rebuild with Measurements
MERGE-BLOCKERS
Findings (non-blocking, severity order)
Contract checks that PASSED scrutiny
Rebase notes (08f8603 supersession)
|
…rossing machinery The shipped 2-D seam-crossing capability did not survive measurement: the cut drags partition-seam vertices onto the fault line, so the crossing design's premise (one crossing vertex, both flanking facets rank-local) fails and the shipped ptest geometry died at np=2 with a chain-fragment error and skipped at np=4 (PR #502 review, blocker 1). Rather than repair the most intricate code in the delta, give 2-D the strategy 3-D already ships and measures healthy at np=2-8: BEFORE the split, move the fault's cell star plus one growth layer to the rank that already owns most of it (shell partitioner; everything else keeps the balanced partition), so the chain is rank-interior and the split runs with serial topology. _redistribute_fault_interior is now dimension-generic and takes a list of (label, value) pairs; add_fault redistributes ONCE for a whole network, keyed on every fault together, which is also what keeps prior pairings valid (a pairing cannot yet migrate through a redistribution — split_fault still refuses that per-split case, now in both dimensions). In 2-D the move only happens when the chain actually touches the seam (_fault_labels_touch_seam, collective verdict); 3-D stays unconditional. The crossing apparatus this makes unreachable is deleted: the crossings fan walk in _take_sides, the _sf_sum census and third-rank corner rule, _pull_seam_vertices_onto_crossings, and the dup_new keyed star-forest exchange in reconnect._rebuild_point_sf (back to the simple renumbering broadcast; _install_point_sf stays). The low-level splitter now refuses ANY shared chain vertex, checked FIRST so a straddling fault always gets the honest seam verdict instead of a fragment symptom — the misdiagnosis the review measured. ptest_0845's crossing test now asserts SUCCESS (split + slip, machine-zero leak) with no tolerated refusals — measured 4 passed at np=2 and np=4, 3 passed + 1 documented low-level skip at np=3, peak slip rank-independent to ~1e-8 against serial. The organic straddling cut keeps its collective low-level refusal and now also asserts that add_fault splits the same geometry cleanly. 3-D ptest_0848 stays 9/9 at np=2 and np=4 with the generalised helper. Underworld development team with AI support from Claude Code
fault_normal_traction(solver, "A", info) on a two-fault network returned the UNION of both faults' pair nodes (measured: 44 rows where fault A has 23) — fault B's nodes interleaved into A's traction profile as bogus arc-length samples — and update_fault_state mixed the faults the same way in its slip vector, its theta store (foreign point keys written into this fault's by_point) and its convergence monitor (PR #502 review, blocker 2). The cause: _InterfaceAssembler always walked EVERY registered law fault plus the include fault. When include is given the assembler now restricts its walk to exactly those faults — a per-fault view. The SOLVE path passes no include and still assembles every law-carrying fault, untouched. Single-fault runs are bit-identical, which is why the suite never saw the bug. The new regression test (two-fault J0 network, viscous laws on both) asserts the traction profile's row count against _fault_pair_nodes — the per-fault oracle that was always correct — the arc-length span, and that a theta update for fault A writes only fault A's own points. Validated as a negative control: it fails on the unfixed assembler with exactly the measured union count (44 vs 23) and passes after the fix. Also freshens comments that still described the retired 2-D seam-crossing machinery. Underworld development team with AI support from Claude Code
Three raises inside collective paths were rank-local: the pair-straddles-ranks RuntimeError and the fault-node-on-a-wall- boundary ValueError in build_rotation's pair-block loop, and the pairing/closure-mismatch and vanishing-analytic-normal raises inside _fault_pair_nodes (reached from both build_rotation and the interface assembler). A rank that raised alone left its peers blocked in the next collective — Q assembly, the trace-mass exchange, the Newton loop — turning an error into a deadlock at np > 2 (PR #502 review, finding 3). Funnel them through the allgather-verdict pattern fault_split already uses: _collective_raise gathers each rank's local error string and raises the SAME error on every rank if any is non-empty. build_rotation wraps its whole fault block (fault_names is registration state, identical across ranks, so every rank reaches the verdict exchange); the assembler collects _fault_pair_nodes failures per fault and verdicts after its walk, before the collective mass exchange. The straddle check should now be structurally impossible — the redistribution makes every pair rank-local — but checks that "cannot fire" are exactly where this PR's blocker 1 lived, so it stays, and its verdict is now collective like the rest. Underworld development team with AI support from Claude Code
Initialised in the solver setup, never written or read anywhere in the repository — the interface laws actually live in solver._fault_interface_laws, set lazily by utilities/fault_contact.py (PR #502 review, finding 6). Deleted, and the neighbouring comment now points at the real home of the laws. Underworld development team with AI support from Claude Code
Member
Author
|
Response commits 1d48731..c0d64cd — both blockers fixed, plus findings 3 and 6:
Serial suite 28/28 (the new regression test included). Remaining known work before merge: this branch is stacked on #488 (which owns the two failing CI tests) and carries the superseded 08f8603 form of the #500 fix — both resolved at the rebase once #488 lands. Underworld development team with AI support from Claude Code |
The junction policy Louis set: declared master faults keep geometric continuity; every other junction is left as an offset gap carrying DAMAGE-ZONE material, and the stress lobes decide how the faults link. Two pieces land here: - prepare_fault_network(..., return_junctions=True) hands out every junction's kind, point, pull-back and fault pair. - uw.meshing.damage_zone_yield(mesh, junctions, tau_damage, radius) builds the composite yield-stress expression over all plugs (SHARP Piecewise regions folded with Min against a sane finite far cap). The recipe behind it is measured (~/+Simulations/fault_junction_rheology): a von Mises plug's tau_J dial is smooth over a decade and interpolates the geometric-decomposition brackets; a collinear gap grows a straight, sharp-edged yielded band tip-to-tip (verified at 6 cells across the gap with adapt-on-top refinement through the gap); TWO Picard passes match picard=8 to 0.01% (find the lobe, respond); blending a rheological parameter against a huge sentinel through any smooth mask tail contaminates the plug — hence the sharp regions; the fully-weakened band (~100x) locally recreates the thin-weak-inclusion conditioning cost, so production plugs should cap weakening near 10x. Newton via the consistent tangent on the rotated path is issue #507 (pre-existing, wall-only reproducer, NOT a split-node defect). Underworld development team with AI support from Claude Code
…eature/fault-split-node
Measured side by side on the collinear-bridge strain rate (the maintainer's challenge): continuous-P1 nodal projection renders smoothly with no ringing — the old "never P1" rule was calibrated against the kinked-normal sawtooth fields later cured at the source, and the split mesh's P1 space is naturally discontinuous across faults (doubled nodes), so the jump is representable. P0 remains right where the field carries genuine cell-scale structure: the plastic yield-zone boundary is sharp in P0 and blurred a node spacing by P1. Rule now: choose by physics; components recovered separately, invariants in numpy; one-sided Oranges on white for strain rate with mesh edges visible (the maintainer's convention, restored). Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
…split-node # Conflicts: # src/underworld3/utilities/rotated_bc.py
…#507) Differentiating an unwrapped strain-rate invariant produces half-integer powers of (grad v : grad v) whose value or derivative is 0/0 at a state of rest, so EVERY consistent-tangent assembly with edot = 0 anywhere filled the operator with NaN — measured J(0) norm = nan for a ViscoPlastic model at ANY yield stress, surfacing as GAMG's "Computed maximum singular value as zero" (error 77) on the standard, rotated free-slip, and split-node fault paths alike. Mainline, since the consistent-tangent feature existed. Why it hid, and why detection cannot fix it (both measured): - Layer 1's automatic cold-start Picard injection fires correctly and STILL fails on boundary-driven problems: one nrichardson sweep propagates boundary data a single element layer, leaving the deep interior at exactly zero strain rate. Body-force-driven problems (the yield campaigns) fill F(0) everywhere — the only class the protection ever actually covered. - The "continuation" tangent's alpha = 0 phase does not protect: the blended kernel evaluates the Newton branch pointwise and IEEE 0*NaN = NaN. - A rigidly-translating stuck region has edot = 0 at the CONVERGED solution: the singular state is physics, not a start-up artifact, so no warm-start policy can make it unreachable. The fix implements the derivative's removable-singularity limit: _jacobian_unwrap now adds 1e-36 under every half-integer power whose argument carries unknowns (+1/2 the invariant, -1/2 its reciprocal in eta_pl, -3/2 their derivatives). Jacobian sources only — the residual never passes through this function and the default Picard tangent never calls it, so both stay bit-identical by construction. The perturbation at any resolvable strain rate is under one part in 1e24; the cold Jacobian of a never-yielding probe now equals its smooth-state Jacobian to all digits. The Layer-1 warm-up remains as the convergence aid it was designed to be; its comment now states the measured facts. Gates: 53 passed across yield-homotopy, rotated free-slip, fault contact/API, plus the new cold-start regression test (test_1067_newton_cold_start: True and "continuation", solve from cold and assert J(0) finite). Every previously-failing reproducer (standard / rotated wall / fault, True / continuation) now solves. Underworld development team with AI support from Claude Code
solve_rotated_freeslip assumed the native solve() dispatch had already attached the mesh auxiliary vector to the solver DM. Callers that enter the loop directly (fault_contact.solve_with_fault) bypass that preamble, so any form referencing an auxiliary MeshVariable read a NULL aux array in the first FE kernel and segfaulted; where it did not crash, a repeat solve would have read stale auxiliary values after the field changed between calls. The loop now runs update_lvec + setAuxiliaryVec + _update_constants itself at entry (redundant and cheap when the native dispatch already did). Found by the damage-gated junction study: cell-wise damage D as a P0 auxiliary in the plug yield stress crashed identically under Newton and Picard, on split meshes only, through solve_with_fault only. Regression test covers the direct-entry solve and the repeat-solve freshness. Underworld development team with AI support from Claude Code
One object carries the validated network recipe end to end: raw traces + a seniority hierarchy in, then prepare (junctions to offset form, junior severed by senior — prepare_fault_network gains a pairwise hierarchy= alongside the absolute through=), build (graded mesh along every raw trace, split-node faults along the prepared pieces), apply_contact (no-opening pairs on all pieces), damage_yield (the junction glue: one dial, strength and xi-regularisation together, sharp Piecewise plugs sized to the junction), solve and slips (per-piece diagnostics). Defaults encode the measured rulings from the junction study (~/+Simulations/fault_junction_rheology): 2-element ligaments (1-6 all transmit within ~7%, resolution-converged at 2), dial 0.05 (zone stress proportional to the dial to ~100x contrast; the compact plug never triggers thin-inclusion conditioning), isotropic glue deliberately (TI weak-plane fabric cannot relax the corner-turning deformation a junction must accommodate — measured). User guide: docs/advanced/fault-networks.md. Tests: hierarchy flips who is severed; end-to-end network solve transmits slip on every piece under pure Newton from cold. Underworld development team with AI support from Claude Code
… glue The prepare-first architecture carries to 3-D and dissolves the hard meshing problem: the junior patch is trimmed back from the plane-plane intersection segment BEFORE meshing, so gmsh only ever embeds DISJOINT patches — surface-surface intersection never reaches the mesher, exactly as shared vertices never reach add_fault in 2-D. - fault_network_3d: planar-polygon geometry core (plane basis, convexity, half-plane clip, line-polygon clip) + crossing_segment (the junction curve, clipped to both rims) + prepare_fault_surfaces (pairwise hierarchy, ligament band trim in the junior's own plane, minimum-area gate, loud refusals: non-planar, non-convex rims, near-miss-without-crossing). - BoxInternalPatch embeds a LIST of patches (FaultSurfaces and/or (name, polygon) pairs): one physical label each, one embed call, gmsh Distance grading from every patch, and a crossing-patch guard that points the user at the preparer. Single-patch path unchanged. - FaultNetwork accepts FaultSurfaces: same API end to end in 3-D (prepare -> multi-patch build -> sequential split_fault -> contact -> tube glue -> per-patch slips). damage_yield is now dimension-general (full strain invariant; disc plugs for 2-D point junctions, TUBE plugs for 3-D segment junctions); slips uses the dimension-general tangential decomposition. Measured end to end: two crossing patches, hierarchy trims the junior, pure Newton from cold converges (125 s at h=0.06), the drive-aligned senior dominates the slip budget. Tests: analytic crossing-segment oracle, trim clearance, refusals, embed guard, end-to-end solve. Full fault battery green (29 tests incl. the unchanged single-patch 3-D suite). Underworld development team with AI support from Claude Code
_write_coordinates now takes the coordinate dimension and the actual coordinate VALUES (not the source dm and row indices) — both fault_split call sites updated with a fancy-index through _coords. Dropped the stale Null_Boundary mention from the split docstring (the label is gone, #503). Post-merge gates: fault battery 51 green (0844-0851, 1067), rotated/ MG-bundle/placed-surface 66 green (1060-1066, 1021, 0853). Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What this is
Zero-thickness split-node faults for UW3, in 2-D and 3-D, as a user-facing
capability. A fault becomes a genuine velocity discontinuity in the mesh:
the labelled facet chain (2-D) or triangulated patch (3-D) is split by
duplicating its interior nodes, a strong no-opening constraint ties the
coincident pairs together, and interface constitutive laws set the
tangential condition. There is no thin weak band and no viscosity contrast
for the solver to fight: the conditioning benchmark measures 10 outer
Schur iterations against 147 for the equivalent one-element weak inclusion
at 1e-4 contrast on the same mesh.
The user surface
Mesh.add_faultplaces tips and kinks onto mesh vertices, cuts aconforming chain, splits it, and records the coincident DOF pairing.
Lists of faults build offset networks in one call.
add_fault_bc(conds, boundary):0= frictionless,eta_f= viscousinterface;
add_coulomb_fault_bcandadd_rate_state_fault_bcfor thefriction family. Laws are sympy expressions in canonical symbols
(slip rate, effective normal stress, state); consistent tangents come
from
sympy.diff, never hand-coded. The effective normal stress isSIGNED and reaction-fed; each law clamps its own strength (a cohesive
law keeps strength into mild tension, a cohesionless one does not).
normal=on every law variant supplies a smooth fault normal — ananalytic sympy expression, or
"trace"to build it from the fault's ownpolyline. This matters for sampled curved traces: the default
facet-averaged normal zig-zags at sampling kinks and produces slip
notches and traction sawteeth that grow under refinement; the smooth
normal removes them (measured 7-17x) and restores h-convergence.
(
fault_slip,fault_pair_jumps,fault_normal_traction) — the sidesare geometrically coincident, so coordinate queries cannot see them.
BoxInternalPatchembeds a planar interior patch,now with
patch_cellSize/grading_distanceto grade resolution fromthe fault outward.
Parallel
vertex, manufactured automatically; a fault along a seam refuses
loudly and collectively. Swept at np 2-5. A seven-fault network splits
and solves at np 1-8 with per-fault slip rank-independent to ~1e-5.
split_faultredistributes first — the patch's cell star plus onegrowth layer moves to the rank that already owns most of it (shell
partitioner, SF-propagated marking), everything else keeps the balanced
partition. Splits at np 2-8 with pair topology identical to serial;
measured imbalance 1.8x at np 8 (the thin star, not the refined band).
Several 3-D faults on one parallel mesh are refused (pairing migration
through the redistribution is future work).
coordinate DM as well (
_install_point_sf) — previously a parallelHDF5 save of any rebuilt mesh wrote shared vertices as owned on every
rank and the file would not reload. The run-parallel / render-serial
workflow (solve at any np,
write_timestep, read checkpoints seriallyfor analysis) is exercised end to end by the King study scripts.
Validation
converging to the (8/3 pi)(d tau/eta) sqrt(a^2-r^2) oracle from below;
leak at machine zero (1e-17 or exactly 0) in every benchmark.
stick/slide switching with reaction-fed normal stress (2-D and 3-D);
rate-state smoke tests with exact ageing updates.
test_0845(split topology),test_0846(contact physics),test_0847(user API, analytic/trace normals with negative control),test_0848(3-D), plus parallelptest_0845/ptest_0848.docs/advanced/split-node-faults.md),method and benchmark write-up
(
docs/developer/design/SPLIT_NODE_FAULT_METHOD_2026-08.md), afault-mechanics teaching page with nine committed figures and three
animations, and the deployment design doc.
Known limitations (all refuse loudly)
Closed-loop faults; daylighting (faults reaching the domain boundary);
junctions sharing vertices; multiple 3-D faults on one parallel mesh.
Split meshes carry no geometric-MG tail (coarse levels do not contain
the fault) — solvers take algebraic-multigrid defaults. The bilateral
no-opening constraint holds a fault shut where it would physically open;
the tell is a tensile fault normal traction, documented in the guide.
Coordination
This merge brings
reconnect.py/line_cut.py(the cut and rebuildmachinery) to development for the first time. The placed-surface branch
also modifies
reconnect.py; per the agreed protocol, whichever branchmerges second reconciles — this PR merging first means the
placed-surface session rebases onto development and folds
_install_point_sfinto its rebuild paths.An adversarial review is posted as a comment below, per our review
practice.
Pre-release follow-ups (agreed direction, not in this PR)
Before importing generic fault models: (a) lead with the
uw.meshing.Surfaceidiom in examples/docs (accepted and tested today;the examples use the
(name, points)shorthand) and wireFaultSurfaceinto the 3-D split path; (b) source
normal=from the surface objects'own geometry rather than user formulae; (c) an importer that detects
crossing/abutting traces and auto-converts them to the offset-junction
(ligament) form, loudly — true shared-vertex junctions need their own
design (a crossing vertex takes four coincident copies and a non-binary
pairing) and are staged behind a ligament-sensitivity measurement.
Underworld development team with AI support from Claude Code