diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md index db7ff2cc1..9fec8795a 100644 --- a/.claude/skills/adapt-on-top-faults/SKILL.md +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -359,6 +359,7 @@ Measured guidance, all at matched cell count: | FMG under rotated free-slip on a plain (non-adapt) refined mesh | the rotated KSP resolves hierarchies via `custom_mg.build_transfers`: an adapt child's mesh-owned tail is picked up AUTOMATICALLY (#467 fixed the old silent GAMG fallback), but the native `dm_hierarchy` is still not read — on a plain refined mesh, `set_custom_fmg(..., field_id=0)`. | | NVB at np>1 raises NotImplementedError | native `_nvb_transform` extension not built (needs the custom-PETSc/amr env). Both `nvb` and `edge_split` are otherwise fully parallel, 2-D and 3-D. | | high stress appears in the matrix beside the fault | elements STRADDLING the weak-zone margin: one end sees high strain rate, the other high viscosity. The FE forms `mean(η)·mean(ε̇)`; the honest cell average is `mean(η ε̇)`, and the difference is `−2 Cov(η, ε̇)` across the cell. Zero for any cell wholly in or wholly out. See the band-width section below. | +| scattered 1-level refinement across the WHOLE domain; far-field quality drops | the metric's far clip ``h_far`` sits below the base mesh's cell DIAMETERS (gmsh ``cellSize`` is a target edge length; diameters run 1.2–2.5x it — measured 0.108–0.223 for cellSize 0.18+ref 1). ``edge_split`` marks on diameter, so 100% of the domain refines once and the unrequested bisection DE-CONDITIONS the grid (far-field median q 0.372 -> 0.295, measured). Set ``h_far >= 1.05 * cell_diameters(base.dm).max()``. | | refinement band narrower than the fault's INFLUENCE | measured: η still 0.07 at d=0.06 while the mesh has already coarsened 4×, so the artefact peaks on the transition flank, not on the fault. **85 % of it sits at d>0.01.** Size the flat core from the *influence* width, not the fault. | | `uw.function.evaluate` fails "Total components 8 != 6" | cached-interpolation mismatch on a mesh already carrying several solver variables. Sample the field numerically from `surface.unsigned_distance` instead. | | bare SIGSEGV, no traceback, after building a Mesh from a raw DM | `uw.discretisation.Mesh(dm, ...)` TAKES THE DM OVER. Read geometry from `child.dm`, never the handle you passed in. | diff --git a/.claude/skills/uw-visualisation/SKILL.md b/.claude/skills/uw-visualisation/SKILL.md index 87cba71fb..1a43f6ea5 100644 --- a/.claude/skills/uw-visualisation/SKILL.md +++ b/.claude/skills/uw-visualisation/SKILL.md @@ -56,6 +56,15 @@ pl.screenshot(out); pl.close() high-order DOFs → blocky. **Attach `T.data[:,0]` directly** to the DOF-cloud mesh from `meshVariable_to_pv_mesh_object` (it is correct for annulus/box/disc — do NOT avoid it). `clim` MUST be passed (default `clim=""` trips `np.any`). +5. Resampling ANY field (even P1) onto a regular pixel grid via + `uw.function.evaluate` **dapples at element boundaries** — grid points that + straddle a facet get located into a neighbouring cell with slightly-off + reference coords (Louis: "artefacts across the elements", S-fault rig). + Render derived fields NODALLY on the mesh's own triangulation instead: + evaluate at `mesh_to_pv_mesh(mesh).points` (exact at vertices for P1, + whichever cell the locator picks), attach as point_data, let VTK + interpolate WITHIN elements. On a SPLIT mesh never Delaunay the DOF cloud + (it re-triangulates across the slit) — use the mesh's own cells. ## Seeing the MESH (adaptation / moving mesh) diff --git a/docs/advanced/stress-visualisation.md b/docs/advanced/stress-visualisation.md new file mode 100644 index 000000000..0c186620f --- /dev/null +++ b/docs/advanced/stress-visualisation.md @@ -0,0 +1,197 @@ +# Visualising the Stress Tensor + +Scalar fields get colormaps and velocity gets arrows; the stress tensor +needs its own glyph. `underworld3.visualisation` provides +principal-stress glyphs — sampled at seed points, the way velocity +arrows sample the velocity — and stress trajectories, the curvilinear +net traced by the principal directions. Both work from the same +recovered stress fields you already checkpoint. + +## Principal-stress glyphs + +At each seed point we diagonalise the stress and draw one bar per +principal axis: bar length proportional to the principal-value +magnitude, blue for compressive ($\lambda < 0$), red for tensile. In +2-D this is the classical stress cross; in 3-D each seed carries three +orthogonal bars, best drawn on one or two section planes rather than a +filled volume. + +```python +import underworld3 as uw +import underworld3.visualisation as vis +import sympy + +# Recovered stress components (P1 projections of the deviatoric +# stress) plus pressure give the full stress. Project once after the +# solve; do not pass raw solver derivative expressions to a plot. +# A scalar variable's .sym is a 1x1 Matrix — index it before +# assembling the tensor. +stress = sympy.Matrix([[Txx.sym[0] - P.sym[0], Txy.sym[0]], + [Txy.sym[0], Tyy.sym[0] - P.sym[0]]]) + +pl = vis.plot_stress_glyphs(mesh, stress, num_seeds=24, + save_png=True, dir_fname="stress_glyphs.png") +``` + +Seeds default to a regular grid over the mesh bounding box, filtered +to points inside the mesh (an annulus seeds no glyphs in its hole). +Pass `seeds` explicitly to sample section planes in 3-D or to avoid +regions — in fault models, keep seeds out of the weak zones, where the +recovered stress mixes materials across the interface. + +```{figure} figures/stress_glyphs_thrust.png +:alt: Two-panel figure for a thrust-ramp model. Top panel, principal stress crosses on a 26 by 13 seed grid over a 2 by 1 section with two grey parabolic ramps rising from a basal decollement, drawn over a faint peach strain-rate wash in which the quiet wedges riding each ramp glow white and the blind ramp tips smudge darker. Far-field crosses are blue-horizontal (compression); in the wedges the crosses turn red and near-vertical (tension), largest just above the blind tips. Bottom panel, stress trajectories over the same wash: dark sigma-1 lines run horizontally and arch smoothly over each ramp tip, pale sigma-3 lines rise near-vertically between them, crossing at right angles everywhere. +:name: fig-stress-glyphs-thrust + +Principal-stress crosses (top) and stress trajectories (bottom) for a +blind-thrust model, over a faint strain-rate wash (log scale, +background-percentile limits, low opacity — the second invariant of +the strain rate from the same recovered stress). The colour +convention matches the RdBu_r field convention used across the +documentation: blue compressive, red tensile. The wash ties the +stress geometry to the deformation it drives: the near-rigid thrust +wedges glow white while the blind tips concentrate both quantities. +``` + +### The pressure gauge matters for colours, not directions + +For incompressible models the full stress is +$\sigma = \tau - p\,I$ and the pressure datum is a gauge choice. +Shifting that datum shifts every principal value equally, so it can +flip bars between red and blue — but it cannot rotate the principal +directions or reorder the principal values. State the gauge in your +caption (the examples here demean the pressure), or add the +lithostatic reference before plotting if absolute compression +matters. + +The same argument answers a common question about map-view regime +colouring in the style of the World Stress Map (red normal, green +strike-slip, blue thrust): the regime classification is +gauge-invariant, but in a 2-D incompressible plane-strain model it is +also *degenerate* — the out-of-plane deviatoric stress is zero while +the in-plane deviatoric principals are $\pm s$, so the out-of-plane +stress is always the intermediate principal stress and every map-view +point classifies as strike-slip. Regime colouring only carries +information in 3-D models. + +## Stress trajectories (2-D) + +Trajectories integrate the principal *direction* field into curves — +the classical stress-trajectory diagrams of structural geology. A +principal direction is defined only modulo 180°, so ordinary +streamline tools cannot draw this field: the integrator in +`direction_trajectories` sign-aligns each evaluated eigenvector with +the previous heading, and places lines evenly (Jobard–Lehmann +occupancy) so the figure stays legible. The two families cross at +right angles wherever both are drawn — a built-in correctness check. + +The direction callable is yours to build, which keeps the integrator +independent of how the stress is stored. From nodal arrays: + +```python +import numpy as np + +def make_direction(interpolate_stress, family="compressive"): + # interpolate_stress(p) -> (sxx, syy, sxy) at point p, or None + def direction_at(p): + values = interpolate_stress(p) + if values is None: + return None + sxx, syy, sxy = values + mean_dev = 0.5 * (sxx - syy) + radius = np.hypot(mean_dev, sxy) + if radius < 1.0e-12: # isotropic point: direction undefined + return None + angle = 0.5 * np.arctan2(sxy, mean_dev) # most-tensile axis + if family == "compressive": + angle += 0.5 * np.pi + return np.array([np.cos(angle), np.sin(angle)]) + return direction_at + +lines = vis.direction_trajectories( + direction_at, candidate_seeds, inside, + step=0.008, separation=0.04, +) +trajectory_lines = vis.trajectories_to_pv_lines(lines) +``` + +Draw the compressive family dark and the tensile family pale, as in +the figure above. In 3-D the analogue of a trajectory is a surface; +we do not attempt those — draw glyphs on section planes instead. + +## Two more regimes, same recipe + +The same two panels for the extensional and strike-slip companions of +the thrust model — nested listric normal faults above a detachment, +and en-echelon segments in dextral simple shear. Nothing changes in +the code except the checkpoint being loaded. + +```{figure} figures/stress_glyphs_listric.png +:alt: Two-panel figure for three nested listric normal faults soling into a basal detachment, over a faint peach strain-rate wash with white quiet triangles in the footwall beneath each sole. Top panel, principal stress crosses: red-horizontal tension dominates the upper plate and deepens between the fault traces; blue compression concentrates near the detachment, with strongly rotated mixed crosses hugging each curved fault. Below the detachment the crosses return to uniform red-horizontal. Bottom panel, trajectories: dark sigma-1 lines hang near-vertically in the extending upper plate, bending to meet each listric sole, while pale sigma-3 lines run horizontally; beneath the detachment the net is an undisturbed rectangular grid. +:name: fig-stress-glyphs-listric + +Nested listric normal faults (extension): horizontal tension aloft, +compression concentrating under the soles, and an undisturbed +trajectory grid beneath the detachment — the faults decouple the two +plates. +``` + +```{figure} figures/stress_glyphs_enechelon.png +:alt: Two-panel map-view figure for three en-echelon fault segments in dextral simple shear, over a faint strain-rate wash with bright white lobes at the relay steps between segment tips. Top panel, principal stress crosses at 45 degrees far from the faults, shrinking and rotating through the relay zones, with a blue compressive bridge linking overlapping tips. Bottom panel, trajectories: the dark sigma-1 family sweeps diagonally across the box and kinks sharply as it hands across each relay step; the pale sigma-3 family crosses it orthogonally. +:name: fig-stress-glyphs-enechelon + +En-echelon segments in dextral simple shear (map view): the conjugate +cross pattern of the far field collapses through the relay steps, and +the trajectory families kink as stress hands across from segment to +segment. +``` + +## 3-D glyphs + +`principal_stress_glyphs` accepts `(n, 3, 3)` tensors and returns +three bars per seed. Seed one or two planes through the feature of +interest: + +```python +u = np.linspace(0.05, 0.95, 13) +gx, gz = np.meshgrid(u, u) +plane = np.column_stack([gx.ravel(), + np.full(gx.size, 0.5), # y = centre plane + gz.ravel()]) +pl = vis.plot_stress_glyphs(mesh, stress, seeds=plane) +``` + +```{figure} figures/stress_glyphs_sinker3d.png +:alt: A unit cube drawn in outline with a grey sphere of radius 0.16 near the centre, slightly above mid-height. Three-bar principal stress glyphs are drawn on a vertical section plane and a horizontal section plane through the sphere, via the one-call plot_stress_glyphs with the cube edges and sphere added to the returned plotter. Above the sphere the bars are red and near-vertical (tension as material is pulled down behind the sinker); below and beside it they are blue (compression), fanning outward on the horizontal plane beneath the sphere. Bar length decays with distance from the sphere. +:name: fig-stress-glyphs-sinker + +Three-bar principal-stress glyphs on two section planes through a +Stokes sinker: a tensile (red) column above the sinking sphere, a +compressive (blue) fan below and around it. Rendered with +`plot_stress_glyphs(..., show=False)` so the cube outline and sphere +could be added to the returned plotter before the screenshot. +``` + +## Building blocks + +The plot function is a convenience wrapper; every step is available +separately for custom figures: + +| Function | Purpose | +|---|---| +| `tensor_fn_to_pv_points(pv_mesh, uw_fn)` | Evaluate a `dim`×`dim` sympy tensor at points | +| `principal_stress_glyphs(coords, stress, scale)` | Bar segments + `"tensile"` cell array | +| `direction_trajectories(direction_at, seeds, inside, step, separation)` | Evenly spaced mod-180° trajectories | +| `trajectories_to_pv_lines(lines)` | Bundle polylines for `add_mesh` | +| `plot_stress_glyphs(mesh, stress, ...)` | One-call cross plot | + +The figures on this page come from checkpointed models (the three +fault-interaction examples and a Stokes sinker): the solve writes the +mesh, the velocity, the pressure, and the recovered stress components +with `mesh.write_timestep`, and the glyph plots load them back with +`read_timestep` — no re-solving to restyle a figure. The faint +background in the 2-D panels is the strain-rate second invariant from +the same recovered stress ($\dot\varepsilon = \tau/2$ at matrix +viscosity 1), drawn first at low opacity with log scaling and +background-percentile colour limits, so the glyphs carry the figure +and the wash only whispers where deformation concentrates. diff --git a/docs/developer/design/fault-patch-multigrid-2026-08.md b/docs/developer/design/fault-patch-multigrid-2026-08.md new file mode 100644 index 000000000..c9d13d8ce --- /dev/null +++ b/docs/developer/design/fault-patch-multigrid-2026-08.md @@ -0,0 +1,401 @@ +# Fault-patch multigrid: the measured skeleton (2026-08, #629) + +The solver architecture for fault-system Stokes problems, settled by +measurement on a 2-D analysis rig and ready to carry into 3-D. One +sentence: **ordinary balanced multigrid solves the smooth problem — +refined or not, layered or not — and a set of small strong subdomain +solves, keyed on the fault physics and proportional to the fault trace, +restores the multigrid contract where coarse spaces cannot see.** + +This note records the architecture, the evidence, the configuration +rules that were each learned from a failure, and the parallel design +rules for the 3-D deployment. The campaign log with every number is +`~/+Simulations/place_route_health/README.md` (items 17–22) and the +issue thread is #629; the code landed on `feature/fault-outcrop-3d-cap` +(commits `ec213489` → `d4b86ea0`). + +## The problem + +A fault zone — split-node contact, a painted weak band, or a painted +transversely-isotropic band — concentrates the difficulty of the Stokes +operator into a thin cross-cutting region. Two measured pathologies: + +1. **Coarse corrections cannot see the band modes.** Level ablation at + viscosity contrast 1e4: full / no-mid / two-level tails converge in + 20 / 19 / 24 iterations — tail depth is irrelevant. Even a level that + carries the band geometry contributes nothing: an unsplit coarse + space cannot represent slip, and coefficient-blind transfers smear + the weak inclusion. +2. **Convergence is then proportional to fine-level smoothing work.** + Doubling the smoother (gmres/8 for /4) exactly halves the iteration + count. So the count grows with contrast — 7 → 10 → 36 over + 1e0 → 1e4 — where a healthy multigrid holds constant V-cycles. + +GAMG fails the same way but harder (56 → 178 over the same range, its +iteration cap by 1e6): the sharp thin inclusion is the malignant kind of +coefficient structure for any coarsening that has to guess. + +## The architecture + +* **Multigrid over everything.** The custom-P Galerkin ladder (or GAMG + where no hierarchy exists) with whole-level smoothing. This handles + every smooth mode, including in refined regions and across *extended* + (co-dimension-1, layered) viscosity jumps — measured: a sharp 1e4 + strong lid crossed by the fault leaves the iteration count unchanged, + with the lid boundary entirely outside any patch. +* **Strong patches where the physics is sharp.** The fine-level smoother + gains one ASM component whose subdomains are the fault zone(s): the + slit trace for a split fault, the painted cells for a weak or TI band. + Each block is solved directly (small shifted LU) once per smoother + application. The patch is *never* keyed on the refinement structure — + in a global subduction mesh the refined region is most of the node + budget, while the fault zone scales as the trace, + m ≈ N^(2/3) in 3-D. + +With this in place the velocity block converges in **4–9 iterations +independent of viscosity contrast from 1 to 1e6**, of the number of +along-strike blocks the trace is divided into (k = 1…16 measured), of a +strong lid of either jump or smooth profile, of fault-network topology +(a V-junction runs identically as one merged block or per-segment +blocks overlapping at the apex — the junction is a non-event), and of +the fault representation: split contact (6–7), unsplit painted weak +band (9 → 7 over 1e4 → 1e6, GAMG 40+), and painted TI band (4, GAMG +43). One cross-subsystem observation: with the velocity preconditioner +made whole, the pressure sub-solve dropped from its 200-iteration cap +to 14 — part of the recorded #625 pathology is downstream of an +incomplete velocity PC, and the Schur diagnosis should be re-baselined +against a healthy velocity block. + +**Confirmed in 3-D** (2026-08-23, pure contact, the composed +place_thin_volume ribbon + authored slit + split at 19,911 cells): 7 +velocity iterations against GAMG's 56, physics identical (slip 0.1509, +leak 0 both) — and identically 7 with the zone patch, with the +structural patch alone, and with **no finest patch at all** (the cover +gate declines the band-shaped structural patch, leaving whole-level +smoothing): under pure contact the patch is redundant, consistent with +the keying ruling above. The remaining 3-D gap is wall-clock, not +iterations, and its cause is now precise: **the 2-D composed stack's +band levels nest by construction — the 3-D ones share nothing.** The +2-D transfinite ladder at 2:1 spacing (36 → 72 rungs, rails and spine +at identical offsets) makes every mid-band vertex a fine-band vertex, +so the placed pairs behave like native refinement: native densities +(24–28 nnz/row), a small structural patch. The 3-D ribbon layers are +two independent unstructured gmsh fills at unrelated sizes (0.025 / +0.015 — 5:3), so every band row is a non-identity transfer row: fat +Galerkin products (119–173 nnz/row against the native ~90 of the +finest), a band-wide "structural" patch (~82% of the level, correctly +gated off), and an expensive V-cycle — GAMG keeps warm wall time at +rig scale (17.5 vs 48 s). The engineering consequence: to give 3-D the +2-D economics, the band levels must nest — build the fine ribbon as a +refinement of the mid ribbon (one fill, refined 2:1) or as a +structured/extruded slab, rather than as two unrelated fills. + +**The extruded ladder, measured (2026-08-24), and the corrected +diagnosis.** `place_thin_volume(mesher="ladder")` now builds the 3-D +band from the fault sheet's own discretisation — the grid offset ±w/2 +along one set of analytic normals into two prism layers, split to tets +(Dompierre; quad-diagonal compatibility proven by the analytic-skin +face count), no remesh — and it is exact where exactness is checkable: +planar volume to 1e-12, curvature inversion a refusal, and a 2:1 +subsampled grid's band sharing every vertex with the fine band all the +way through placement and split. Measured on a curved fault +(sinusoidal bulge): physics identical to GAMG, 6–7 iterations against +60–65. But vertex nesting alone did NOT collapse the chain, for two +measured reasons that correct the paragraph above: **the fill shell, +not the band, is the node majority** — the unstructured annulus gmsh +builds between the band skin and the coarse background carried 51% of +the finest level's P2 nodes (nesting 3%) against the band's 23%; and +**vertex nesting is not P2 nesting** — only 10% of the band's P2 DOFs +coincide (edge-midpoint DOFs differ between scales). Thinning the +shell (`clearance=0.3`, the 0.6·w floor governing) brought level 1 to +native density and the count to 6. + +**Correction (2026-08-24, profiled): the warm-time gap is NOT a +transfer rebuild.** The rotated path already caches the transfers and +the whole KSP/PC context across solves (`_rotated_linear_cache`: +geometry tier carries Q and `custom_Pl`; the operator tier skips +assembly and PCSetUp on an unchanged matrix) — a cProfile of the +repeat solve shows its ~46 s entirely inside the composite Krylov +solve, with no build anywhere. An earlier campaign note described the +pre-fix state and was repeated here without re-measuring; retracted. +The warm cost structure is the Schur loop times the V-cycle price: +full Schur factorization runs one velocity solve per pressure +iteration (~23 velocity solves × 6–7 V-cycles at the recorded +0.25 s/cycle ≈ the measured 46 s), against GAMG's 0.015 s cycles +(≈ 23 s). The wall-clock levers are therefore level economics (the +item-15 two-level result: near-fine-sized intermediate levels are +negative value at rig scale) and the per-cycle cost of the robust +smoother stack — plus the #625 pressure-iteration count under +contrast. **Measured on the ladder stack (tail-depth A/B, physics +identical in every arm): dropping the mid level reaches GAMG parity +warm (24.0 s vs GAMG 23.0) at 7 iterations against GAMG's 65, and the +remaining chain is entirely native density (81/75/93 nnz/row) — the +Galerkin fat was the mid level alone. The mid level's role is +production proportions (a real gap between background h and band +spacing); at rig proportions the tail of choice is [L0, L1] + finest. +With parity warm and 8x iteration headroom, any contrast regime wins +outright. The augmented-Lagrangian penalty then shortens the Schur +loop on top — at γ=1: pressure 22 → 17 with the velocity count +unmoved, warm 19.1 s against GAMG's 23.0, FMG ahead outright at flat +viscosity (the #625 mechanism: AL pays when the velocity PC is +healthy; no Schur plumbing was needed — the penalty is +viscosity-scaled, so the 1/K mass preconditioner stays spectrally +correct by design). γ=10 over-stiffens the velocity block (7 → 12) +and reverses the gain. Two caveats travel with the penalty: the +recovered p is the multiplier (use p_mech = p − λμ∇·u for +pressure-dependent rheology), and derived quantities drift — peak +slip 1% at γ=1, 3% at γ=10 at rig resolution, the same bias family as +#633's vertex-sampled dynamic topography. γ=1 is the sanctioned +choice.** Remaining mesh-side levers unchanged: graded outer sheets +on the ladder; P2-aware nesting last. + +The user-facing entry (the production ruling: conceptual simplicity is +the feature — a band built from the user's own fault-surface mesh +cannot be mis-prepped) is **`place_fault_ribbon(base, sheet_grid, +width)`**: one call, one parametrisation — derived normals, the +mid-surface labelled and split (a boundary ready for `add_fault_bc`), +the nested 2:1 bridge level returned alongside, thin-shell clearance +by default, `split=False` for painted weak/TI models on identical +geometry. + +**The grid is the fault, honoured exactly (ruling, 2026-08-24 — a +fault specified from a structural model cannot be shrunk to make a +tip margin).** The band satisfies the tip rule by EXTRAPOLATION: the +sheet is continued `margin_rings` rings outward along its own end +tangents (linear — no curvature the data never asserted), the slip +surface is labelled on precisely the supplied points, and the intact +extrapolated frame is the 3-D analogue of the junction ruling's +intact gap. Measured: the earlier inset spelling was confiscating two +rings of the requested fault and biasing peak slip by 14% (0.1523 → +0.1737 honoured, matching the full-patch embed reference ~0.174); +only the slip patch's two corner triangles erode (the splitter's +no-interior-vertex refusal — a half-cell nick). + +## Configuration rules (each learned from a failure) + +0. **The finest patch always contains the structural patch.** The patch + smoother *replaces* whole-level smoothing, so it inherits every row + the coarse level cannot represent — the cut/split-inserted DOFs that + the transfer's non-identity rows identify — whether or not the + physics zone covers them. A zone away from the cut leaves those rows + smoothed nowhere and the velocity sub-solve caps on every + application (measured three independent ways; the overlap-0 slit + stagnation was the same omission). The zone-block builder unions the + uncovered structural rows in as an additional block automatically. +1. **Basic ASM, never restricted.** Discarding the halo part of the + subdomain correction stalls the outer Krylov at 80–375 iterations + where basic runs 6 (whole-level baseline 4). `PC_ASM_RESTRICT` is the + wrong variant for this use. +2. **Shifted factorization is mandatory on split-contact patches.** The + rotated Galerkin patch block carries near-zero pivots inherited from + the constraint-zeroed transfer rows (min diagonal ~1e-5); an + unshifted factorization takes `NUMERIC_ZEROPIVOT` even as exact LU + and the KSP dies at 0 iterations with `PC_FAILED`. Use + `sub_pc_factor_shift_type nonzero`. +3. **SOR subdomain solves are the safe default; LU is the strong + option.** On the near-singular split-contact patch the LU is what + buys contrast constancy. On an unsplit painted band the two are + indistinguishable at rig scale. +4. **Patch keying.** Split meshes: coincident fine pairs — the split's + plus/minus nodes at bit-identical coordinates. (Not transfer identity + rows onto coarse nodes: a cut inserts *new* vertices at edge + crossings that coincide with no coarse node, and that detection + declines silently.) Painted models: an explicit boolean cell mask — + the modeler painted the band and knows the cells; nothing needs + detecting. **Which key for which model (ruling, 2026-08-24): the + zone key (`fac_zone`) is for VOLUMETRIC fault representations only** + — a weak / TI / damage zone whose width is physics. A split-node + fault runs the structural patch alone: splitting exists to be + efficient, its patch is the trace, and it is installed + automatically. The ribbon is never the key — it is resolution, not + rheology (measured: under pure contact the zone block is redundant — + 8 structural-only vs 9 zone-keyed iterations in 2-D — while inflating + the finest block from the trace to most of the band). +5. **Blocks per fault-network connected component; junction cells in + every adjacent block.** Distant segments couple through the medium, + which is smooth and belongs to the multigrid. Along-strike division + of one fault into blocks is measured free (7 iterations at every k), + so blocks can be sized for convenience and parallel placement. +6. **The smoother variant matters under contrast.** The fast + (richardson) smoother degrades where gmres holds (12 → 23 already at + 1e2); contrast work runs the robust variant. +7. **Check the ASM probe, not just the iteration count.** A patch mode + that declines silently leaves a plausible-looking whole-level solve; + two instrumentation bugs in this campaign were caught only because + the probe line (subdomain sizes, failed-reason codes) disagreed with + the narrative. + +## Two transfer-level results the ladder depends on + +* **Structural zeros fattened the Galerkin chain.** The point-located + transfer builders emit ~1e-16 weights that are structural nonzeros, + and PtAP fills by structure, compounding per level (90 → 265 → 481 + nnz/row measured down a 4-level tail). Dropping them + (`_drop_structural_zeros`) collapses the chain to native-like density + (138/119/173/90) and halved the warm solve. This was the largest + single win of the campaign. +* **Exact nested transfers for native `refine()` pairs.** Level pairs + tagged as native refinement get the true any-degree FE embedding + (dual-basis `W = B·M⁻¹` under the parent's affine pullback), exact to + 1e-15 at P1/P2 in 2-D/3-D. PETSc's `DMCreateInterpolation` general + path is **not** the embedding (row sums to 1.375, 1e-2 error on + quadratics) and must not be used for transfers. + +## Scaling + +The patch reduces the only direct solve from N volume DOFs to the trace +m ≈ N^(2/3). Because iterations are flat in the number of along-strike +blocks, fixed-size segment blocks make the total factorization cost +**linear in the trace and embarrassingly parallel** — no superlinear +direct-solver scaling survives. Under nonlinear fault rheology the +economics improve further: the tangent changes in the patch, so the +per-Newton refresh is the small factor, and the background ladder can be +lagged. (The end-state for stiff friction laws is nonlinear Schwarz / +ASPIN on the same subdomains; nothing in the architecture has to change +to permit it.) + +## Parallel design rules for 3-D + +1. **Split in the local frame, after distribution — then pair + co-residency is automatic.** A pair is the duplication of one owned + facet's nodes and is born on that facet's rank; no partitioner + constraint is needed (Louis's correction of an earlier draft of this + rule). The two real items instead: (a) **never distribute a + pre-split mesh** — the slit is a zero-cost cut in the mesh graph, so + a partitioner would *preferentially* separate the sides while the + contact coupling lives outside the graph; gate that pipeline rather + than engineer around it. (b) **Seam consistency**: where the fault + crosses a rank boundary along strike, on-rank duplication must + update the star forest consistently (the known np>=3 line-cut issue; + the independent-pass design is the template). Along-strike block + cuts at seams are measured free (the k-ladder). +2. **PCASM is rank-local.** Blocks become per-rank lists with local + indices; faultless ranks carry `nsd = 0`, and every IS-building path + must survive empty sets (the empty-stratum `getIndices()` SEGV + family, #589). +3. **Pair detection must go collective or be replaced by zone keying.** + Rank-local coordinate matching misses pairs that straddle a boundary + and can false-match ghosts; the zone mask is rank-local by nature and + is the recommended production key. +4. **Test at np ≥ 4 with timing probes** — np = 2 is a special case + that passes collectives which hang at 4 (#512/#596 discipline), and + `mat_block_size` divisibility is BC-dependent and mixed across ranks + (#584): gate collectively. +5. **Serial-only pieces that need parallel forms before 3-D-parallel:** + the FAC patch split, the nested-pair transfer detection, and the + zero-column repair (the parallel guard currently raises where the + serial path repairs — placed levels with orphan columns would fail + at np > 1). +6. **The pressure block is the wall-clock gate at contrast** (its gasm + sub-solve caps at 200 iterations, silently unconverged) and gasm is + the most parallel-fragile piece in the stack — #625, being addressed + separately; whatever lands there should be chosen with parallel + behaviour in mind. + +## Status of the knobs + +The campaign ran on TODO(MEASURE)-marked environment variables +(`UW_FAC_*`, `UW_CUSTOM_MG_*`, `UW_MG_SMOOTH_ITS`) plus one explicit +hook (`solver._fac_zone_cells`). The hook has since settled into API +(2026-08-23): the patch key is `set_custom_fmg(..., fac_zone=...)` — a +boolean cell mask or a list of per-segment masks, validated against the +finest mesh at registration — and the `_fac_zone_cells` attribute +spelling is retired loudly (setting it raises rather than declining +silently). `mesh.cells_labelled(name, value)` builds the mask from a +placement's cell label empty-safely (the #589 null-IS segfault is fixed +at source in `utilities/dm_labels.py`, and the dead `is None` stratum +guards in `nvb`/`reconnect`/`fault_split` are routed through it). The +transfinite ladder band is in-repo as +`place_thin_volume(..., mesher="ladder")`, and the composed 2-D +benchmark is enshrined as `tests/test_1022_composed_ribbon_fmg.py` +(tier B: level densities, 2-block finest patch, iteration bounds, +slip/leak). Still environment variables, to be settled by the remaining +measurements: the `UW_FAC_*` sub-solver choices and the +`UW_CUSTOM_MG_*` / `UW_MG_SMOOTH_ITS` A/B flags. + +## What the fault ribbon is for (ruling, 2026-08-23) + +The finite-width ribbon around a fault is a **modelling object in its +own right**, not just mesh scaffolding. Its sanctioned uses: + +* **A damage zone** — when a damage *evolution equation* is part of the + physics. Damage is then a solved field whose rheology follows from + the field; it localises where the mechanics puts it. Damage does not + exist in a model that does not evolve it. +* **A nonlinear plastic / yielding material** — the ribbon's material + given a yield rheology, so that **failure patterns can emerge within + the resolved band**. This is the mechanism by which junctions *form*: + localisation finds its own geometry inside the ribbon instead of the + modeller authoring it. (The gradient-plasticity caveats apply when + this is built: the regularising length must enter the plasticity, and + mesh objectivity needs more than one resolution — the ribbon width + supplies a geometric length scale but does not by itself regularise.) +* **A permeable zone** — a localized pathway for fluid flow, carrying + the permeability structure a fault zone has and the surrounding rock + does not. + +And a fault may itself be **represented as a weak or TI weak zone** — +a legitimate model in its own right, carried by the ribbon. The +discipline that separates it from "fake damage" is the status of the +width: **in a weak/TI-zone fault model, w is a physical parameter** — +it *is* the fault-zone width, must be chosen appropriately, and the +mesh must resolve it (~2 elements across whatever w is chosen). **In a +split-node model, w is a mesh-bridging convenience** — a resolution +transition with no physical reading, and nothing rheological may be +keyed to it. The incoherent thing — the error this campaign made three +times — is mixing the readings: split nodes *plus* band-wide weakness +with w chosen as a mesh number is neither model. Pick the +representation; w's meaning follows from the choice. In a model with +none of the physics above and a split fault, the ribbon carries +background rheology and is purely resolution. + +These roles compose with the solver design below unchanged: whichever +physics the ribbon carries, the band is the authored patch key, the +levels are the FMG bridge, and the machinery is measured +contrast-robust for whatever coefficient structure the physics +produces. + +## The ribbon is part of the FMG (ruling, 2026-08-22) + +The placed ribbon — the finite band of fault-scale resolution around +each fault, with its modest adaptation — is not an alternative to the +patch methodology; it is **part of the multigrid design**, playing +three roles at once: + +* **The bridge between the standard mesh and the patch.** The ribbon's + band is the natural intermediate structure: coarse ladder below, + band-resolved level(s) in the middle, the cut/split finest on top. + Its cell burden is minimal by construction (the band is a few local-h + wide), and its levels are exactly the composite-grid structure the + FAC machinery was built for. +* **The damage-zone identifier — authored, not detected.** The + placement machinery labels the band's cells, and the zone label + survives the cut/split as cell children. The strong patch keys on + that label directly: an exact, mesh-conforming, modeler-authored + region — never the ragged staircase of a distance mask, never the + one-cell zipper comb of ``cells_supporting`` around the split. +* **One mesh discipline for every fault representation.** The same + ribbon carries split-contact faults (cut down the transfinite + centreline — three nodes across, so nothing snaps, #595), painted + TI, or painted weak rheology; the representation choice stops + dictating the mesh strategy. + +Junctions follow the revised ruling: split segments stop a short +distance apart and the **intact gap is the linkage** — no painted core +by default (a damage core is opt-in physics). The ribbon simply +encloses the whole stepping system, gaps included, at fault resolution. + +The open engineering question is **parallel balance**: the ribbon +concentrates cells, so the partitioner must weight them. With split +surgery done in the local frame (rule 1), pair residency is automatic +and this reduces to ordinary cell-count weighting, uncoupled from any +fault-topology constraint. + +## The 2-D rig + +`~/+Simulations/place_route_health/line_mg.py`: 2-D box, interior line +cut, optional split + frictionless contact, painted weak / TI / lid +viscosity structures, V-junction fixture, ~2.5k cells, under a second +per warm solve. Every configuration reproduces with one flag. New +solver-design questions should be answered there first; 3-D is for +confirmation, not exploration. diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index 8da83e480..4ed69e815 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -772,3 +772,9 @@ def _sendData(): # Note: SymPy converter registration approach doesn't work reliably in strict mode # The better approach is to ensure UWexpression arithmetic operations return SymPy objects # This is handled by the __rmul__, __radd__ etc. methods in the mathematical mixin + +# The environment-armed hang watchdog (UW_HANG_WATCHDOG), armed HERE — +# after every module above has finished importing — because +# faulthandler's repeating C dump against a still-importing interpreter +# was measured to livelock or SIGSEGV (see mpi._watch_from_environment). +mpi._arm_environment_watchdog() diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c42bbb474..3b3da609f 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6085,19 +6085,21 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 - # The outer Krylov must be FLEXIBLE, because both sub-blocks below are - # themselves Krylov solves run to a tolerance. Inside the Schur - # factorisation the velocity sub-solve computes the search direction, so - # the operator the outer method applies differs from one outer iteration - # to the next, and plain GMRES's residual recurrence assumes it does not. - # PETSc's default is `gmres`, and that default is invisible on easy - # problems -- isoviscous SolKz converges in two outer iterations -- and - # expensive on hard ones: on the Spiegelman notch at refinement 3, 983 - # velocity iterations per step and DIVERGED_LINEAR_SOLVE, against 58 for - # fgmres and 23 with a loosened inner tolerance. Raising the velocity - # iteration cap changes nothing (byte-identical residuals), so the failure - # is inconsistency rather than too few iterations. See #576. - self.petsc_options["ksp_type"] = "fgmres" + # The OUTER Krylov must be FLEXIBLE: both sub-blocks below are + # themselves Krylov solves run to a tolerance, so the operator the + # outer method applies differs from one outer iteration to the next + # and plain GMRES's residual recurrence does not hold — the + # velocity-block FGMRES reasoning (#147), one level up. PETSc's + # default is `gmres`; invisible on easy problems (isoviscous SolKz + # converges in two outer iterations), ruinous on hard ones: 14,400 + # inner iterations / ~540 s on an 85k-cell contrast problem for both + # velocity preconditioners (#624), and on the Spiegelman notch at + # refinement 3, 983 velocity iterations per step and + # DIVERGED_LINEAR_SOLVE against 58 for fgmres — raising the velocity + # cap changes nothing (byte-identical residuals), so the failure is + # inconsistency, not iteration count (#576). Managed, so an explicit + # user ksp_type still wins. + self._push_managed_option("ksp_type", "fgmres") self.petsc_options["pc_type"] = "fieldsplit" self.petsc_options["pc_fieldsplit_type"] = "schur" @@ -7017,9 +7019,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["snes_ksp_ew"] = None self.petsc_options["snes_ksp_ew_version"] = 3 - # Flexible for the same reason as in __init__ (#576): the sub-blocks are - # inexact Krylov solves, so the outer operator varies between iterations. - self.petsc_options["ksp_type"] = "fgmres" + # Flexible for the same reason as in __init__ (#576/#624): the + # sub-blocks are inexact Krylov solves, so the outer operator varies + # between iterations. Managed, so an explicit user ksp_type wins. + self._push_managed_option("ksp_type", "fgmres") self.petsc_options["pc_type"] = "fieldsplit" self.petsc_options["pc_fieldsplit_type"] = "schur" diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 85c751fce..769dd6fa1 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2712,6 +2712,11 @@ def nuke_coords_and_rebuild( uw.pprint(f"PETScDS - (re) initialised") self._coord_array = {} + # Cleared with _coord_array because both describe the same node layout: + # the coordinates and which cell owns each of them. Topology only, so it + # survives node motion — but it is invalidated by the same re-creation + # of the DS that invalidates the coordinates. + self._cell_node_array = {} # let's go ahead and do an initial projection from linear (the default) # to linear. this really is a nothing operation, but a @@ -5466,13 +5471,14 @@ def _get_coords_for_var(self, var): self._coord_array[key] = self._get_coords_for_basis(var.degree, var.continuous) return self._coord_array[key] - def _get_coords_for_basis(self, degree, continuous): - """ - This function returns the vertex array for the - provided variable. If the array does not already exist, - it is first created and then returned. - """ + def _basis_coordinate_dm(self, degree, continuous): + """Coordinate DM carrying a degree-``degree`` Lagrange field. + Its local section defines the node layout that + :meth:`_get_coords_for_basis` reads, so anything that needs to know + WHICH node is which — as opposed to just where the nodes are — has to + come from this same DM. The caller destroys it. + """ dmold = self.dm.getCoordinateDM() dmold.createDS() dmnew = dmold.clone() @@ -5493,6 +5499,73 @@ def _get_coords_for_basis(self, degree, continuous): dmnew.setField(0, dmfe) dmnew.createDS() + dmfe.destroy() # DMSetField took its own reference + return dmnew + + def _cell_node_indices(self, degree, continuous): + """Rows of each cell's degrees of freedom, indexing exactly the array + :meth:`_get_coords_for_basis` returns for the same ``(degree, + continuous)``. + + Returns an ``(n_cells, nodes_per_cell)`` integer array; row ``k`` lists + cell ``cStart + k``'s DOF rows in closure order. The order within a row + is arbitrary but self-consistent, which is all any consumer needs: the + weights are computed from the coordinates read at these same rows, so + no reference-element node ordering is ever assumed. (Element assembly + WOULD care — the two DOFs on an edge follow the edge's own orientation, + not the cell's — so do not repurpose this for that.) + + Knowing a cell's nodes is what turns the adapt parent-CELL map into an + exact transfer at any polynomial degree: for a Lagrange element the + basis is dual to its nodal points, so a coarse cell's own DOF + coordinates determine the coarse interpolant inside it (#425). + """ + from math import comb + + key = (self.isSimplex, degree, continuous) + if key in self._cell_node_array: + return self._cell_node_array[key] + + dmnew = self._basis_coordinate_dm(degree, continuous) + section = dmnew.getLocalSection() + cStart, cEnd = self.dm.getHeightStratum(0) + cdim = self.cdim + + rows = [] + for cell in range(cStart, cEnd): + cell_rows = [] + for point in dmnew.getTransitiveClosure(cell)[0]: + ndof = section.getDof(point) // cdim + if ndof: + offset = section.getOffset(point) // cdim + cell_rows.extend(range(offset, offset + ndof)) + rows.append(cell_rows) + dmnew.destroy() + + expected = comb(degree + self.dim, self.dim) + if not self.isSimplex or any(len(r) != expected for r in rows): + # A tensor-product Q_k cell carries (k+1)^dim nodes, so the monomial + # basis of total degree <= k would not be square against them and + # the dual-basis construction does not apply. Say so here rather + # than return a ragged array the caller has to second-guess. + got = sorted({len(r) for r in rows}) + raise NotImplementedError( + f"_cell_node_indices needs a simplex mesh: expected " + f"{expected} nodes per cell for degree {degree} in {self.dim}D, " + f"got {got}") + + self._cell_node_array[key] = numpy.asarray(rows, dtype=numpy.int64) + return self._cell_node_array[key] + + def _get_coords_for_basis(self, degree, continuous): + """ + This function returns the vertex array for the + provided variable. If the array does not already exist, + it is first created and then returned. + """ + + dmold = self.dm.getCoordinateDM() + dmnew = self._basis_coordinate_dm(degree, continuous) matInterp, vecScale = dmold.createInterpolation(dmnew) coordsOld = self.dm.getCoordinates() @@ -5513,7 +5586,6 @@ def _get_coords_for_basis(self, degree, continuous): if vecScale is not None: vecScale.destroy() dmnew.destroy() - dmfe.destroy() return arrcopy @@ -7107,10 +7179,22 @@ def _wrap_coarse_level(self, dm): def _coarse_level_meshes(self): """The static coarse-mesh tail (one Mesh per base hierarchy level, coarsest..base-finest), built once and cached — they never change - because the base hierarchy is static across adapts.""" + because the base hierarchy is static across adapts. + + Each wrap is tagged with its ``(hierarchy token, level)`` slot — and so + is this mesh itself, whose ``dm`` is the finest hierarchy level. Two + levels whose slots are consecutive under the same token are a native + ``refine()`` pair, which is what lets ``custom_mg`` give that pair the + EXACT nested prolongation instead of a point-located one (#425/#629). + The token is a plain sentinel object: identity ties the family together + without a reference cycle back to this mesh.""" cached = getattr(self, "_coarse_level_meshes_cache", None) if cached is None: cached = [self._wrap_coarse_level(d) for d in self.dm_hierarchy] + token = object() + for k, w in enumerate(cached): + w._refine_slot = (token, k) + self._refine_slot = (token, len(cached) - 1) self._coarse_level_meshes_cache = cached return cached @@ -7563,13 +7647,76 @@ def cells_supporting(self, name): if label is None or label.getStratumSize(value) == 0: return zone + fS, fE = dm.getHeightStratum(1) for f in label.getStratumIS(value).getIndices(): - # `_cells_on_edge` rather than `getSupport` directly: in 2-D an edge - # IS a facet and its support is already the cells, but in 3-D the - # support holds faces and the cells are one level further up. - # Applying the 2-D walk in 3-D returns nothing at all, silently. - for c in _cells_on_edge(dm, int(f)): - zone[c - cS] = True + p = int(f) + if fS <= p < fE: + # A FACET's support is the cells, in any dimension — a 2-D + # facet is an edge, a 3-D facet is a face. (This method only + # ever saw 2-D meshes before the 3-D conforming sheet, and + # the edge walk below returns NOTHING from a 3-D face, + # silently — the same trap `_cells_on_edge` documents, one + # level up.) + for c in dm.getSupport(p): + if cS <= c < cE: + zone[c - cS] = True + else: + # A labelled EDGE in 3-D (a trace chain): cells are one + # level further up. + for c in _cells_on_edge(dm, p): + zone[c - cS] = True + return zone + + def cells_labelled(self, name, value=None): + """The cells carrying DM label ``name`` (optionally stratum ``value``). + + The cell-label partner of :meth:`cells_supporting`: where that method + derives a zone from a FACET label (a conforming surface's support), + this one reads a CELL label directly — the region a placement call + painted (:func:`~underworld3.utilities.place_surface.place_thin_volume` + labels its layer's cells), an imported region marker, or any other + authored cell stratum. It is the empty-safe way to build a fault-zone + patch key for ``set_custom_fmg(..., fac_zone=...)`` (#629): an absent + label or an empty stratum returns an all-``False`` mask rather than + touching the null IS that segfaults ``getIndices()`` (#589). + + Parameters + ---------- + name : str + The DM label name (e.g. the ``label=`` given to a placement call). + value : int or None + The stratum value; ``None`` takes the union over every value the + label carries. + + Returns + ------- + numpy.ndarray + Boolean, one entry per cell, in **plex cell order** — also the DOF + order of a ``degree=0`` :class:`MeshVariable`. Points of the label + that are not cells (faces, edges) are ignored. + + Examples + -------- + >>> zone = mesh.cells_labelled("Band") + >>> set_custom_fmg(stokes, tail, field_id=0, fac_zone=zone) + """ + from underworld3.utilities.dm_labels import label_stratum_indices + + dm = self.dm + cS, cE = dm.getHeightStratum(0) + zone = numpy.zeros(cE - cS, dtype=bool) + label = dm.getLabel(name) + if label is None: + return zone + if value is None: + vis = label.getValueIS() + values = [int(v) for v in vis.getIndices()] if vis is not None else [] + else: + values = [int(value)] + for v in values: + pts = label_stratum_indices(label, v) + pts = pts[(pts >= cS) & (pts < cE)] + zone[pts - cS] = True return zone @staticmethod @@ -7608,6 +7755,87 @@ def _repair_cut(cut_dm, lines, info, reach, verbose): f"removed, min angle now {info['min_angle']:.2f} deg") return cut_dm, info + def _adopt_cut_child(self, cut_dm, boundaries, info, mg_coarsening_ratio, + verbose): + """Wrap a DM cut at the finest level as this mesh's child. + + Shared by :meth:`add_conforming_surface` (2-D line cut) and + :meth:`add_conforming_sheet` (3-D placed sheet): the cut exists on + the finest level only, so this mesh plus everything below it is the + child's coarse multigrid tail. The coarse levels do not carry the + surface and do not need to — see the measured note in + :meth:`add_conforming_surface`. + """ + child = Mesh( + cut_dm, + simplex=self.dm.isSimplex(), + coordinate_system_type=self.CoordinateSystem.coordinate_type, + qdegree=self.qdegree, + boundaries=boundaries, + verbose=False, + ) + child.parent = self + child._relationship_kind = "refinement" + # ... but NOT a nested one. The cut moves or replaces parent vertices + # (the 2-D cut snaps them onto the surface; the 3-D carve deletes and + # refills), so a coarse DOF need not have a coincident fine DOF, and + # the injection that a bisection child's restriction relies on would + # quietly read the field at the displaced position instead. + child._refine_dofs_coincide = False + child.regions = self.regions + child._parent_mesh_version = self._mesh_version + child._surface_info = info + + # Mesh-owned custom-P geometric-MG tail. Adding a surface refines this + # mesh, so this mesh plus everything below it is a valid coarse tail and + # the solver appends the child as the finest level. The transfers are + # coordinate-based and do not need the levels to nest — just as well, + # since a cut vertex is not an edge midpoint and the exact 1/2,1/2 + # prolongation does not apply to it. + # + # A mesh that is ITSELF a child (a second surface, or an adapt child) has + # to EXTEND its own tail rather than read `dm_hierarchy`, which for a child + # holds only its own DM: reading it there would silently discard every + # level below and leave a two-level hierarchy calling itself multigrid. + # + # Tested with `is not None`, not for truthiness. A child whose own tail + # is EMPTY is still a child, and reading `dm_hierarchy` there returns + # just its own DM — the two-level collapse this comment warns about, + # reached by the one input the truth test cannot distinguish from a + # parent. + own_tail = getattr(self, "_custom_mg_coarse_meshes", None) + tail = (list(own_tail) + [self]) if own_tail is not None \ + else self._coarse_level_meshes() + + # A cut is not necessarily a refinement. It re-represents the same grid + # with the surface conformed, so `self` earns its place as a separate + # level only if the child is genuinely finer — the same question `adapt` + # asks of an engine pass, so ask it with the same routine rather than a + # second rule that could drift from it. + # + # Measured on a box fault before this: nine levels, of which the two + # added by the two cuts coarsened h by 1.11x and 1.17x on the 5th + # percentile against a threshold of 1.8 — each costing a full Galerkin + # RAP and smoother sweep for no correction. Worse, transfer 7->8, BETWEEN + # those two, is where the barycentric builder ran out of coarse DOFs with + # a fine image and fell back to the dense RBF one (#424). + # + # `_subsample_mg_levels` already does "replace the level below rather + # than append to it" for its own finest generation; handing it the pair + # (self, child) against the level beneath them puts that decision here + # too. One level back means it kept only the child. + if len(tail) >= 2: + kept, _Ps, _pc = self._subsample_mg_levels( + tail[-2].dm, [tail[-1].dm, cut_dm], [None, None], [], + ratio=mg_coarsening_ratio, verbose=verbose) + if len(kept) == 1: + tail = tail[:-1] + child._custom_mg_coarse_meshes = tail + child._custom_mg_builder = self._custom_mg_builder + + self._registered_children.add(child) + return child + def add_conforming_surface(self, surface, snap_frac=0.10, verbose=False, snap_quality=0.15, snap_dist=0.0, mg_coarsening_ratio=2.0, repair=False, @@ -7772,9 +8000,10 @@ def add_conforming_surface(self, surface, snap_frac=0.10, verbose=False, Notes ----- - Two dimensions only. A surface **ending inside** the mesh (a fault tip) is - refused rather than silently mis-meshed, as is a triangle the surface - crosses three times. + Two dimensions only — in 3-D use :meth:`add_conforming_sheet`, where + a free rim (a fault tip) is the normal case. Here a surface **ending + inside** the mesh is refused rather than silently mis-meshed, as is a + triangle the surface crosses three times. See Also -------- @@ -7814,74 +8043,83 @@ def add_conforming_surface(self, surface, snap_frac=0.10, verbose=False, f"{info['n_cut_edges']} surface facets, " f"min angle {info['min_angle']:.2f} deg") - child = Mesh( - cut_dm, - simplex=self.dm.isSimplex(), - coordinate_system_type=self.CoordinateSystem.coordinate_type, - qdegree=self.qdegree, - boundaries=boundaries, - verbose=False, - ) - child.parent = self - child._relationship_kind = "refinement" - # ... but NOT a nested one. Snapping moves parent vertices onto the - # surface, so a coarse DOF need not have a coincident fine DOF, and the - # injection that a bisection child's restriction relies on would quietly - # read the field at the displaced position instead. - child._refine_dofs_coincide = False - child.regions = self.regions - child._parent_mesh_version = self._mesh_version - child._surface_info = info + return self._adopt_cut_child(cut_dm, boundaries, info, + mg_coarsening_ratio, verbose) + + def add_conforming_sheet(self, points, triangles, name, *, + clearance=0.6, setback=0.0, size=None, + verbose=False, mg_coarsening_ratio=2.0): + r"""Add a triangulated sheet that the mesh conforms to (3-D). + + The 3-D twin of :meth:`add_conforming_surface`, and the Mesh-level + form of :func:`~underworld3.utilities.place_surface.place_sheet`: + the sheet's points become mesh vertices, every sheet triangle an + interior face labelled ``name``, and the rim is free inside the + mesh — a fault tip is the normal case here, not a refusal. + + The cut runs at the finest level ONLY. This mesh and every + multigrid level under it are untouched and become the child's + coarse tail, exactly as in 2-D: the coarse levels do not carry the + sheet and do not need to (custom-P sets ``pc_mg_galerkin=both``, + so every coarse operator is :math:`P^\mathsf{T} A P` from the fine + operator — see the measured note in :meth:`add_conforming_surface`, + whose warning about ESSENTIAL conditions on the surface applies + here unchanged; a material contrast across the sheet is the + supported use). + + Everything :func:`place_sheet` documents holds: the sheet may run + PAST the domain (it is clipped against the mesh's own boundary, + and an outcrop trace is labelled ``_trace``); ``setback`` + stops it short as a BLIND fault with the would-be intersection + returned in ``child._surface_info["surface_trace"]``; ``size`` + re-triangulates the sheet to match the mesh it cuts. - # Mesh-owned custom-P geometric-MG tail. Adding a surface refines this - # mesh, so this mesh plus everything below it is a valid coarse tail and - # the solver appends the child as the finest level. The transfers are - # coordinate-based and do not need the levels to nest — just as well, - # since a cut vertex is not an edge midpoint and the exact 1/2,1/2 - # prolongation does not apply to it. - # - # A mesh that is ITSELF a child (a second surface, or an adapt child) has - # to EXTEND its own tail rather than read `dm_hierarchy`, which for a child - # holds only its own DM: reading it there would silently discard every - # level below and leave a two-level hierarchy calling itself multigrid. - # - # Tested with `is not None`, not for truthiness. A child whose own tail - # is EMPTY is still a child, and reading `dm_hierarchy` there returns - # just its own DM — the two-level collapse this comment warns about, - # reached by the one input the truth test cannot distinguish from a - # parent. - own_tail = getattr(self, "_custom_mg_coarse_meshes", None) - tail = (list(own_tail) + [self]) if own_tail is not None \ - else self._coarse_level_meshes() + Parameters + ---------- + points, triangles : array_like + The sheet: ``(N, 3)`` vertices and ``(M, 3)`` triangle + indices. Explicit arrays rather than an object, because a + sheet is DATA — a slab model, an authored parameter-space + triangulation — whose connectivity must be embedded verbatim + (:class:`~underworld3.meshing.FaultSurface` re-derives its + triangulation, so it cannot carry an authored one). + name : str + Becomes a boundary of the returned mesh, so a solver can + resolve the facets by name and :meth:`cells_supporting` + marks the fault zone. + clearance, setback, size, verbose + Passed through to :func:`place_sheet`. + mg_coarsening_ratio : float + As in :meth:`add_conforming_surface`: the cut replaces this + mesh in the tail unless it is genuinely finer. - # A cut is not necessarily a refinement. It re-represents the same grid - # with the surface conformed, so `self` earns its place as a separate - # level only if the child is genuinely finer — the same question `adapt` - # asks of an engine pass, so ask it with the same routine rather than a - # second rule that could drift from it. - # - # Measured on a box fault before this: nine levels, of which the two - # added by the two cuts coarsened h by 1.11x and 1.17x on the 5th - # percentile against a threshold of 1.8 — each costing a full Galerkin - # RAP and smoother sweep for no correction. Worse, transfer 7->8, BETWEEN - # those two, is where the barycentric builder ran out of coarse DOFs with - # a fine image and fell back to the dense RBF one (#424). - # - # `_subsample_mg_levels` already does "replace the level below rather - # than append to it" for its own finest generation; handing it the pair - # (self, child) against the level beneath them puts that decision here - # too. One level back means it kept only the child. - if len(tail) >= 2: - kept, _Ps, _pc = self._subsample_mg_levels( - tail[-2].dm, [tail[-1].dm, cut_dm], [None, None], [], - ratio=mg_coarsening_ratio, verbose=verbose) - if len(kept) == 1: - tail = tail[:-1] - child._custom_mg_coarse_meshes = tail - child._custom_mg_builder = self._custom_mg_builder + Returns + ------- + Mesh + A new mesh; this one is not modified. Placement metadata is + on ``child._surface_info``. Call again on the result to add + another sheet — a network is built one branch at a time. - self._registered_children.add(child) - return child + See Also + -------- + add_conforming_surface : the 2-D form. + add_fault : cut AND split, for a velocity discontinuity. + cells_supporting : the fault zone of the labelled facets. + """ + from underworld3.utilities.place_surface import place_sheet + + if self.dim != 3: + raise NotImplementedError( + "add_conforming_sheet is 3-D; in 2-D use " + "add_conforming_surface.") + + boundaries = self._boundaries_with(name) + cut_dm, info = place_sheet( + self.dm, points, triangles, label=name, + label_value=boundaries[name].value, clearance=clearance, + verbose=verbose, setback=setback, size=size) + return self._adopt_cut_child(cut_dm, boundaries, info, + mg_coarsening_ratio, verbose) def add_fault(self, faults, verbose=False): """Cut AND split one or more faults; return the split mesh. @@ -8599,12 +8837,28 @@ def _relax_generation(engine_obj, carry, rcarry): # so fewer edges are independent per pass. n_pass = 8 * dim * max_levels current_dm = base_finest + _pct5_of_dm = {} for level in range(n_pass): centroids, _proxy_h, cs = cell_geometry(current_dm) + # The topology tables are read ONCE per pass and shared + # with the split — the per-cell closure walk is the + # loop's dominant cost and was paid twice (#610). The + # metric goes through marking_metric: the "nobody has + # cells" verdict is COLLECTIVE (rank-local branching on + # the raw metric is the np>1 deadlock class). + _edge_tables = None M = marking_metric(centroids) if M is not None and centroids.shape[0]: h_target = 1.0 / numpy.sqrt(M) - diameter = edge_split.cell_diameters(current_dm) + diameter, _edge_tables = edge_split.cell_diameters( + current_dm, return_tables=True) + # The marking pass has just measured this dm; the MG + # level selection re-derives the same 5th percentile + # per retained generation, so it is cached here + # rather than re-walking every level's topology. + _pct5_of_dm[id(current_dm)] = ( + float(numpy.percentile(diameter, 5)) + if diameter.size else float("inf")) sel = numpy.where(diameter > h_target)[0] if node_budget is not None and sel.size > node_budget: order = numpy.argsort(M[sel])[::-1] @@ -8615,7 +8869,7 @@ def _relax_generation(engine_obj, carry, rcarry): marked = [int(cs + j) for j in sel] _coarse_for_P = current_dm current_dm, n_split = edge_split.bisect_longest_edges( - current_dm, marked) + current_dm, marked, tables=_edge_tables) # n_split is global, so this stop is collective without a further # reduction — a rank with nothing marked still enters the split. if n_split == 0: @@ -8649,9 +8903,14 @@ def _relax_generation(engine_obj, carry, rcarry): uw.pprint(0, f"[adapt] edge_split pass {level}: repaired " f"with {n_flips} flip(s)") else: - _nested_parent_cells.append( - None if _vP is None - else _nested_parents(_coarse_for_P, current_dm, _vP)) + # Parent maps are DEFERRED to the retained MG levels: + # the subsampler discarded every per-pass map whose + # span was more than one pass, so building ~76 of + # them to keep 3-4 was almost entirely wasted work — + # and the retained multi-pass spans now get EXACT + # parents from the composed vertex transfer instead + # of falling back to the geometric builder. + _nested_parent_cells.append(None) snap_level_boundaries(current_dm) if _relax_mode == "per-generation": _mg = Mesh(current_dm.clone(), @@ -8842,7 +9101,9 @@ def _relax_generation(engine_obj, carry, rcarry): if level_dms: level_dms, _nested_Ps, _nested_parent_cells = self._subsample_mg_levels( base_finest, level_dms, _nested_Ps, _nested_parent_cells, - ratio=mg_coarsening_ratio, verbose=verbose) + ratio=mg_coarsening_ratio, verbose=verbose, + resolution_hint=(_pct5_of_dm if engine == "edge_split" + else None)) # Exact per-generation prolongations when the engine could supply them # (cell-list path). Empty for the native transform path, which falls @@ -8881,7 +9142,8 @@ def _relax_generation(engine_obj, carry, rcarry): _MG_RATIO_SLACK = 0.9 # a step of 1.92 counts as a doubling def _subsample_mg_levels(self, base_finest, level_dms, nested_Ps, - nested_parent_cells, ratio=2.0, verbose=False): + nested_parent_cells, ratio=2.0, verbose=False, + resolution_hint=None): """Keep one multigrid level per DOUBLING OF RESOLUTION, not one per pass. A refinement engine takes as many passes as it needs to reach the size @@ -8921,10 +9183,16 @@ def resolution(dm): A low percentile rather than the strict minimum, so one thin cell cannot declare a level; reduced with MIN so the finest region counts - wherever it happens to live. + wherever it happens to live. A caller that already measured a + dm during its own pass loop supplies the value through + ``resolution_hint`` instead of paying a second topology walk. """ - d = edge_split.cell_diameters(dm) - local = float(numpy.percentile(d, 5)) if d.size else float("inf") + if resolution_hint is not None and id(dm) in resolution_hint: + local = resolution_hint[id(dm)] + else: + d = edge_split.cell_diameters(dm) + local = (float(numpy.percentile(d, 5)) if d.size + else float("inf")) return uw.mpi.comm.allreduce(local, op=min) # An engine lands near the target, not on it (1.92, 1.97, 2.19 measured), @@ -8954,15 +9222,9 @@ def resolution(dm): composed, parent_cells = [], [] start = 0 + level_coarse = base_finest for i in keep: span = [P for P in nested_Ps[start:i + 1]] - # One generation -> the level IS that pass, so its parent-cell map - # still describes it. More -> no single parent per cell. Not every - # engine records the maps at all (the native transform and SBR paths - # do not), so a short list means "none for this level". - parent_cells.append(nested_parent_cells[i] - if i == start and i < len(nested_parent_cells) - else None) if any(P is None for P in span) or not span: composed.append(None) elif len(span) == 1: @@ -8974,6 +9236,26 @@ def resolution(dm): for P in span[1:]: M = _compose_prolongations(P, M) composed.append(M) + # The parent-cell map, for the RETAINED pair only. A map the + # engine recorded per pass (a single-pass span) is used as + # recorded; otherwise it is derived from the composed vertex + # transfer — nested_cell_parents is topological through the + # transfer, and a descendant's referenced coarse vertices are + # all corners of its ancestor at any depth, so multi-pass + # spans now carry exact parents instead of None. + recorded = (nested_parent_cells[i] + if i == start and i < len(nested_parent_cells) + else None) + if recorded is not None: + parent_cells.append(recorded) + elif composed[-1] is not None: + from underworld3.utilities.nvb import ( + nested_cell_parents as _parents_of) + parent_cells.append( + _parents_of(level_coarse, level_dms[i], composed[-1])) + else: + parent_cells.append(None) + level_coarse = level_dms[i] start = i + 1 if verbose: diff --git a/src/underworld3/meshing/fault_network.py b/src/underworld3/meshing/fault_network.py index b8e1a3d08..d3d888f9f 100644 --- a/src/underworld3/meshing/fault_network.py +++ b/src/underworld3/meshing/fault_network.py @@ -191,12 +191,16 @@ def _build_3d(self, h_far=None, qdegree=2, mesher="embed", The place route on a UNIFORM base is healthy and parallel-validated (ptest_0852: place -> split -> contact, - 22 s including the solve). On edge_split ADAPT CHILDREN the - composed mesh builds and converges but the solve is - pathological (measured 3850 s vs 125 s embed) — an - operator-health interaction between the graded transitions, - the placement cavity, and the split, recorded as an open work - item. Until that is resolved, place is opt-in here. + 22 s including the solve). On edge_split ADAPT CHILDREN it is + SLOW but not sick: re-measured 2026-08-20 (#621), the composed + chain built 27x the cells embed did for the same nominal sizes + and solved 41x slower — proportionate under 3-D Stokes scaling + (per-cell cost comparable, one nonlinear iteration, machine-zero + leak, agreeing slip). The over-build is this route's sizing: the + base box is built at ``cellSize=h_far`` WITH ``refinement=1``, + so the far field is ``h_far/2`` everywhere. Until #621 lands, + place is opt-in here, and expect embed to be much cheaper at + matched request. """ from underworld3.utilities.fault_split import split_fault @@ -217,12 +221,8 @@ def _build_3d(self, h_far=None, qdegree=2, mesher="embed", raise ValueError(f"mesher must be 'place' or 'embed', " f"got {mesher!r}") - from enum import Enum from .cartesian import UnstructuredSimplexBox - from underworld3.utilities.place_surface import (place_sheet, - _sheet_distance) - import underworld3 as uw - from underworld3 import discretisation + from underworld3.utilities.place_surface import _sheet_distance sheets = [(n, *self._triangulate_rim(p, h)) for n, p in self.prepared] @@ -244,23 +244,15 @@ def metric(pts_, _sheets=sheets, _h=h, _hf=h_far, _b=b, _r=r): child = base.adapt(metric, max_levels=max_levels, engine="edge_split") - members = {bd.name: bd.value for bd in child.boundaries} - for k, (n, _sp, _st) in enumerate(sheets): - members[n] = max(members.values()) + 4 - boundaries = Enum("boundaries", members) - - dm = child.dm + mesh = child for n, sp, st in sheets: # clearance 0.8 measured as the working window on # edge_split children (0.6 under-reaches the graded # transition shell and pinches; >=1.0 over-swallows). - dm, _info = place_sheet(dm, sp, st, label=n, - label_value=boundaries[n].value, - clearance=clearance) - mesh = discretisation.Mesh( - dm, simplex=True, - coordinate_system_type=child.CoordinateSystem.coordinate_type, - qdegree=qdegree, boundaries=boundaries, verbose=False) + # Mesh-level placement, so each cut child inherits the adapt + # hierarchy as its coarse multigrid tail (the 2-D contract; + # the split below still forfeits it — see add_fault). + mesh = mesh.add_conforming_sheet(sp, st, n, clearance=clearance) for n, _sp, _st in sheets: mesh = split_fault(mesh, n) self.mesh = mesh diff --git a/src/underworld3/mpi.py b/src/underworld3/mpi.py index 3e1e903d3..0805005a1 100644 --- a/src/underworld3/mpi.py +++ b/src/underworld3/mpi.py @@ -32,6 +32,22 @@ import time as _time from contextlib import contextmanager as _contextmanager +# Pre-import EVERYTHING the watchdog reporter thread can touch. The +# reporter (`_Watchdog.report` -> `_stack_dump` -> traceback/linecache) +# runs on a daemon thread, and with `UW_HANG_WATCHDOG` armed at import +# it can fire WHILE the main thread is still inside `import underworld3`. +# `traceback.format_stack` lazily imports through `linecache` (which +# reaches for `tokenize` on first use); a lazy import on a secondary +# thread while the main thread holds per-module import locks is the +# classic cross-thread import deadlock — measured: a 0.2 s watchdog +# livelocks `import underworld3` locally (the main thread pinned +# mid-import, reporter cycling), and CI's 1.0 s watchdog dies at -11 in +# the same window (test_0054). With these loaded before any reporter +# can run, the reporter never enters the import system. +import traceback as _traceback_preload # noqa: F401 +import linecache as _linecache_preload # noqa: F401 +import tokenize as _tokenize_preload # noqa: F401 + comm = _MPI.COMM_WORLD size = comm.size @@ -313,6 +329,16 @@ def _stack_dump(): interpreter lock, which mpi4py releases around blocking calls -- so a rank sitting in ``allreduce`` does still report. If some extension ever holds the lock through a block, nothing running in Python can report on it. + + ``lookup_lines=False`` is load-bearing, not cosmetic: the reporter runs + on a daemon thread, and with the environment-armed watchdog it can fire + while the MAIN thread is still importing. ``format_stack`` reads every + frame's SOURCE through ``linecache`` — file IO and loader calls against + modules mid-import, from a second thread — and that interleaving was + measured to freeze the import outright (a 0.2 s watchdog livelocked + `import underworld3` 12 times in 15; CI's 1 s watchdog died at -11 in + the same window, test_0054). File names, line numbers and function + names carry the hang report; the source text was the deadlock. """ import traceback @@ -325,10 +351,11 @@ def _stack_dump(): who = "MainThread (this is the one that is stuck)" if ident == main \ else names.get(ident, "unknown") out.append(f"\n --- thread {ident}: {who} ---") - out.extend( - " " + line.rstrip() - for line in traceback.format_stack(frames[ident]) - ) + summary = traceback.StackSummary.extract( + traceback.walk_stack(frames[ident]), lookup_lines=False) + summary.reverse() # walk_stack yields innermost first + out.extend(" " + line.rstrip() + for line in summary.format()) return "\n".join(out) @@ -387,9 +414,25 @@ def arm(self, label=None, resume=False): _faulthandler.dump_traceback_later( self.seconds, repeat=True, file=self.stream, exit=self.abort ) + self._rearm_timer() + def _rearm_timer(self): # Secondary, and only for hangs that leave the interpreter lock free. # It adds the checkpoint label, which faulthandler cannot know about. + # + # The REPORTER re-arms through this method ALONE, never through + # arm(): dump_traceback_later() internally cancels the running C + # watchdog thread and waits on its lock, and when that thread is + # mid-dump — walking frames the main thread is churning (an + # import in progress) — the wait never returns. Measured as the + # test_0054 deadlock triangle (native `sample`): the C thread + # pinned in dump_traceback, the reporter cond-waiting inside + # cancel_dump_traceback_later, the main thread starved in the + # import machinery — 12 of 15 runs frozen at a 0.2 s watchdog. + # faulthandler was armed with repeat=True; it needs no re-arm + # from the reporter. Checkpoints (watch()) still go through + # arm(), where resetting the countdown is the point and the main + # thread is in ordinary running state. if self.timer is not None: self.timer.cancel() self.timer = _threading.Timer(self.seconds, self.report) @@ -422,8 +465,8 @@ def report(self): file=self.stream, flush=True, ) - if not self.abort: - self.arm() # a no-op once cancel() has run + if not self.abort and not self.cancelled: + self._rearm_timer() # timer only — see _rearm_timer def watch(seconds=300, stream=None, abort=False): @@ -743,13 +786,21 @@ def __exit__(self, *args): def _watch_from_environment(): - """Arm the watchdog from ``UW_HANG_WATCHDOG``, at import. - - Arming here rather than from the user's script is the point: a rank that - dies, or diverges, before reaching a ``watch()`` call reports nothing, and - "before the script got going" covers mesh construction and most of the - import graph. This runs as ``underworld3.mpi`` is imported, which is about - as early as anything can. + """Arm the watchdog from ``UW_HANG_WATCHDOG``. + + Arming from the environment rather than the user's script is the point: a + rank that dies, or diverges, before reaching a ``watch()`` call reports + nothing, and arming at import time covers mesh construction and + everything after. It is invoked from the END of ``import underworld3`` + (the bottom of ``underworld3/__init__``), NOT from this module's import — + deliberately: ``faulthandler.dump_traceback_later``'s repeating C dump + walks live frames without synchronisation, and against an interpreter + that is still importing (frames churning, bytecode compiling) that walk + was measured to loop forever or die at SIGSEGV (test_0054 on CI; locally + a 0.2 s watchdog froze ``import underworld3`` on the first piped run). + The price is that a hang INSIDE the import graph itself goes unreported; + everything the tool exists for — meshing, solves, collectives — runs + after import and is covered. ``UW_HANG_WATCHDOG`` Seconds of silence that count as stuck. Unset or 0 disables. @@ -801,4 +852,11 @@ def _watch_from_environment(): #: Path this rank will dump to, or None when the environment did not ask. -environment_dump_path = _watch_from_environment() +#: Set by _arm_environment_watchdog(), called at the END of +#: ``import underworld3`` — see _watch_from_environment for why not here. +environment_dump_path = None + + +def _arm_environment_watchdog(): + global environment_dump_path + environment_dump_path = _watch_from_environment() diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index a933c0102..872194c46 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -38,6 +38,7 @@ operators from the DM hierarchy, so our explicit ``P`` is used. """ +import os from typing import NamedTuple import numpy as np @@ -156,31 +157,48 @@ def barycentric_prolongation(coarse_coords, fine_coords): def rbf_prolongation(coarse_coords, fine_coords, smooth=0.0): - """RBF prolongation: polyharmonic (r² log r) kernel + affine polynomial tail - (reproduces linear fields), Shepard row-normalised to a partition of unity. - Works for arbitrary (non-nested) point sets; software-equivalent to the - barycentric builder as an MG transfer operator.""" - import scipy.sparse as sp - from scipy.spatial.distance import cdist - - def phi(r): - # r² log r → 0 as r → 0; the clamp only keeps log(0) finite at - # coincident points — it does not perturb the kernel value. - r = np.where(r == 0.0, 1e-30, r) - return r ** 2 * np.log(r) - - nc, dim = coarse_coords.shape - Pc = np.hstack([np.ones((nc, 1)), coarse_coords]) # affine tail - Acc = phi(cdist(coarse_coords, coarse_coords)) + smooth * np.eye(nc) - M = np.block([[Acc, Pc], [Pc.T, np.zeros((dim + 1, dim + 1))]]) - B = np.hstack([phi(cdist(fine_coords, coarse_coords)), - np.ones((fine_coords.shape[0], 1)), fine_coords]) - # Solve M Xᵀ = Bᵀ rather than forming M⁻¹ explicitly (faster, more stable). M is - # symmetric, so B M⁻¹ = solve(M, Bᵀ)ᵀ. - Praw = np.linalg.solve(M, B.T).T[:, :nc] - rs = Praw.sum(axis=1, keepdims=True) - rs[np.abs(rs) < 1e-12] = 1.0 - return sp.csr_matrix(Praw / rs) + """RBF prolongation via the STANDARD local interpolator (#429). + + The sparse, linear-exact, kd-tree-based local RBF the rest of the code + uses (``kdtree.interpolation_matrix``, ``order=1``: polyharmonic + r² log r kernel + affine tail solved per target over its nearest + neighbours — see ``docs/developer/subsystems/interpolation.md``). + Replaces the original GLOBAL builder, which assembled and solved the + dense coarse-cloud kernel matrix and returned a transfer with + ``nnz/row == n_coarse`` — a rescue whose Galerkin coarse operators + were dense too, and which had no conditioning path on large clouds + (#429). Same kernel, same reproduction guarantees (constants and + linears to machine precision), sparse support. + + A row-wise kNN builder guarantees nonzeros per ROW, never per COLUMN + (#424): a coarse DOF outside every fine stencil still yields an empty + column and a singular PᵀAP. The build loop's zero-column repair + (nearest-fine-DOF injection) is the counterpart, for this builder and + the barycentric one alike. ``smooth`` is accepted for signature + compatibility and unused — locality is the conditioning here. + """ + from underworld3 import kdtree + + kdt = kdtree.KDTree(np.ascontiguousarray(coarse_coords, dtype=float)) + return kdt.interpolation_matrix( + np.ascontiguousarray(fine_coords, dtype=float), order=1) + + +def _drop_structural_zeros(P_csr, tol=1e-12): + """Remove numerically-zero transfer weights that are STRUCTURALLY nonzero. + + A fine node coincident with a coarse node gets barycentric weights + ``[1, ~1e-16, ~1e-16, ~1e-16]`` — an identity row in VALUE but a 4-entry + row in STRUCTURE, and Galerkin RAP fills by structure. On a composite + (placed/overlay) hierarchy the background is all such rows, so the junk + entries fatten every coarse operator's background block level over level + (measured 90 -> 265 -> 481 nnz/row down a 4-level tail, #629). Weights + are O(1) partition-of-unity values, so ``tol`` cuts only float noise; a + row cannot empty (its weights sum to 1).""" + P = P_csr.tocsr() + P.data[np.abs(P.data) < tol] = 0.0 + P.eliminate_zeros() + return P _BUILDERS = {"barycentric": barycentric_prolongation, "rbf": rbf_prolongation} @@ -404,7 +422,281 @@ def _reduced_transfer(coarse_coords, fine_coords, r2f_c, r2f_f, ncomp, builder): """Build one prolongation reduced(coarse) -> reduced(fine): node-level scalar P -> interleave ``ncomp`` components -> drop BC rows/cols.""" Pn = builder(coarse_coords, fine_coords) # (n_f_nodes, n_c_nodes) - return _reduced_from_node_transfer(Pn, r2f_c, r2f_f, ncomp) + return _reduced_from_node_transfer(_drop_structural_zeros(Pn), + r2f_c, r2f_f, ncomp) + + +def _is_native_refine_pair(coarse_mesh, fine_mesh): + """Are these two levels consecutive members of one native ``refine()`` + hierarchy? Decided by the ``_refine_slot`` tags ``_coarse_level_meshes`` + (and the requested-native arm of :func:`build_transfers`) stamp on the + family — token identity plus consecutive level indices. Untagged levels + (placed meshes, adapt generations, moved bases) answer False and keep + their existing transfer routes.""" + sa = getattr(coarse_mesh, "_refine_slot", None) + sb = getattr(fine_mesh, "_refine_slot", None) + return (sa is not None and sb is not None + and sa[0] is sb[0] and sb[1] == sa[1] + 1) + + +def nested_refine_pair_prolongation(coarse_mesh, fine_mesh, degree, continuous, + coarse_coords=None, fine_coords=None): + """EXACT node-level prolongation for a native ``refine()`` pair, at ANY + polynomial degree — the FE embedding of the coarse Lagrange space in the + fine one. Returns a scipy CSR ``(n_fine_nodes, n_coarse_nodes)`` matrix, + or ``None`` to decline (the caller falls back to the geometric builder). + + Why: the DOF clouds of degree >= 2 spaces do NOT nest even where the + meshes do (an L1 edge node is no L0 node), so the point-located builders + return scattered transfers whose Galerkin product runs 3-5x fatter per row + than a native coarse operator (481 vs ~90-150 nnz/row measured, #629). + PETSc's own ``DMCreateInterpolation`` general path is NOT the embedding + either (measured: row sums to 1.375, quadratic reproduction error 1e-2). + + Construction (#425, the dual-basis identity): a Lagrange basis is dual to + its nodal points, so with ``M[i,m] = mu_m(x_i)`` over a parent cell's own + ``n`` nodes and ``B[t,m] = mu_m(x_t)`` at the fine nodes inside it, the + weight block is ``W = B M^-1`` — no reference element, no per-degree + formulae, no point location. Coordinates are pulled back through the + parent's affine map first: raw monomials on a cell of size ``h`` give + ``cond(M) ~ h^-k``, while under the pullback conditioning depends only on + (degree, dim) — and the pulled-back coordinate is the barycentric vector, + so the "is this fine node inside its parent?" guard is free. + + The parent relation itself is recovered from the two DMs: every fine + vertex of a uniform refinement is an inherited coarse vertex or an exact + edge midpoint (identified by bit equality, so a snapped/relaxed hierarchy + declines rather than guesses), and the parent cell follows topologically + (:func:`~underworld3.utilities.nvb.nested_cell_parents`). Structurally + full rank: every coarse node is inherited into the fine level with weight + 1, so the zero-column failure (#424) cannot arise here. + """ + import itertools + + import scipy.sparse as sp + from underworld3.utilities.nvb import (nested_cell_parents, + nested_prolongation_from_dms) + + cdm, fdm = coarse_mesh.dm, fine_mesh.dm + dim = cdm.getDimension() + if cdm.getCoordinateDim() != dim: # embedded surface: no affine pullback + return None + vP = nested_prolongation_from_dms(cdm, fdm) + if vP is None: + return None + parents = nested_cell_parents(cdm, fdm, vP) + if parents is None: + return None + try: + cn = coarse_mesh._cell_node_indices(degree, continuous) + fn = fine_mesh._cell_node_indices(degree, continuous) + except (NotImplementedError, AttributeError): + return None + + Xc = (np.asarray(coarse_coords) if coarse_coords is not None + else np.asarray(coarse_mesh._get_coords_for_basis(degree, continuous))) + Xf = (np.asarray(fine_coords) if fine_coords is not None + else np.asarray(fine_mesh._get_coords_for_basis(degree, continuous))) + n_c, n_f = Xc.shape[0], Xf.shape[0] + + # One containing fine cell per fine node — any one: the parent interpolant + # is continuous across parent faces, so the weights agree either way. + owner = np.full(n_f, -1, dtype=np.int64) + for k in range(fn.shape[0]): + owner[fn[k]] = k + if (owner < 0).any(): + return None + + ccS, ccE = cdm.getHeightStratum(0) + cvS, cvE = cdm.getDepthStratum(0) + vxy = np.ascontiguousarray( + cdm.getCoordinatesLocal().array.reshape(-1, dim)) + cell_verts = np.empty((ccE - ccS, dim + 1), dtype=np.int64) + for c in range(ccS, ccE): + vv = [q - cvS for q in cdm.getTransitiveClosure(c)[0] + if cvS <= q < cvE] + if len(vv) != dim + 1: + return None + cell_verts[c - ccS] = vv + + # Monomials of total degree <= degree: exactly nodes-per-cell many on a + # simplex, so M is square (guaranteed by _cell_node_indices' own check). + E = np.asarray([e for e in itertools.product(range(degree + 1), repeat=dim) + if sum(e) <= degree]) + if E.shape[0] != cn.shape[1]: + return None + + def vander(L): + return np.prod(L[:, None, :] ** E[None, :, :], axis=2) + + parent_of = parents[owner] - ccS # coarse cell index per fine node + order = np.argsort(parent_of, kind="stable") + runs = np.flatnonzero(np.diff(parent_of[order])) + 1 + rows, cols, vals = [], [], [] + for grp in np.split(order, runs): + c = int(parent_of[grp[0]]) + verts = vxy[cell_verts[c]] + A = (verts[1:] - verts[0]).T + try: + Ainv = np.linalg.inv(A) + except np.linalg.LinAlgError: + return None + Lc = (Xc[cn[c]] - verts[0]) @ Ainv.T + Lf = (Xf[grp] - verts[0]) @ Ainv.T + # The pullback is the barycentric vector: outside the parent means the + # attribution (or the geometry) is off — decline, never extrapolate. + if Lf.min() < -1e-8 or (1.0 - Lf.sum(axis=1)).min() < -1e-8: + return None + try: + W = np.linalg.solve(vander(Lc).T, vander(Lf).T).T # B M^-1 + except np.linalg.LinAlgError: + return None + nodes_c = cn[c] + for t_i, t in enumerate(grp): + keep = np.flatnonzero(np.abs(W[t_i]) > 1e-12) + rows.extend([int(t)] * keep.size) + cols.extend(int(nodes_c[j]) for j in keep) + vals.extend(float(W[t_i, j]) for j in keep) + + P = sp.csr_matrix((vals, (rows, cols)), shape=(n_f, n_c)) + # Lagrange partition of unity: every row must sum to 1 exactly. + if np.abs(np.asarray(P.sum(axis=1)).ravel() - 1.0).max() > 1e-8: + return None + return P + + +def _fac_patch_split(P_csr, coords_c, coords_f, map_c, map_f, nc, + w_tol=1e-8, x_tol=1e-9, cover_max=0.75): + """FAC patch/background split of one level's reduced fine DOFs (#629). + + A locally-refined hierarchy is a COMPOSITE grid: levels differ only in + the refined patch, and the background falls through unchanged. Smoothing + the whole level is the measured pathology (a V-cycle smooths ~175% of the + fine level's nonzeros); the FAC/MLAT answer is to smooth each level only + on its patch plus an interface halo, with the background smoothed once at + the level that owns it. This function finds that patch algebraically from + the level's own transfer — builder-agnostic, so placed levels work too. + + A fine reduced DOF is BACKGROUND when its transfer row is an identity row + (one effective weight-1 entry) AND its node coincides with the referenced + coarse node (within ``x_tol``): the node falls through the pair. Anything + else is PATCH. Split-duplicated nodes — two fine nodes onto one coarse + node, a fault slit — count as patch even though each row is identity: the + pair changes topology there. The HALO is the one-coarse-cell layer of + background whose coarse node is referenced by some patch row. Membership + is decided per NODE (all components together), so the split is invariant + under the rotated path's per-node Q rotation. + + Returns ``(owned_rows, subdomain_rows)`` — reduced fine row indices of + the patch and of patch+halo — or ``None`` when patch+halo covers more + than ``cover_max`` of the level: a uniform pair refines everywhere, and + whole-level smoothing is then the right configuration. + """ + P = P_csr.tocsr() + n = P.shape[0] + if n == 0: + return None + # TODO(MEASURE): #629 campaign knob for the whole-level gate; remove + # once the cover threshold is settled. + cover_max = float(os.environ.get("UW_FAC_COVER", cover_max)) + map_f = np.asarray(map_f, dtype=np.int64) + map_c = np.asarray(map_c, dtype=np.int64) + node_f = map_f // nc + indptr, cols, data = P.indptr, P.indices, P.data + + bg_col = np.full(n, -1, dtype=np.int64) # coarse NODE of background rows + patch_row = np.zeros(n, dtype=bool) + for r in range(n): + sl = slice(indptr[r], indptr[r + 1]) + w = data[sl] + if w.size == 0: + patch_row[r] = True + continue + aw = np.abs(w) + j = int(np.argmax(aw)) + if abs(w[j] - 1.0) > w_tol or (aw.sum() - aw[j]) > w_tol: + patch_row[r] = True + continue + cnode = int(map_c[cols[sl][j]] // nc) + if np.max(np.abs(coords_f[int(node_f[r])] - coords_c[cnode])) > x_tol: + patch_row[r] = True + continue + bg_col[r] = cnode + + nn = int(node_f.max()) + 1 + node_is_patch = np.zeros(nn, dtype=bool) + node_is_patch[node_f[patch_row]] = True + # split duplicates: >1 distinct background fine node onto one coarse node + dup_nodes = np.zeros(nn, dtype=bool) + bgr = np.flatnonzero(bg_col >= 0) + if bgr.size: + pairs = np.unique(np.stack([bg_col[bgr], node_f[bgr]], axis=1), axis=0) + counts = np.bincount(pairs[:, 0]) + dup = np.flatnonzero(counts > 1) + if dup.size: + dup_nodes[pairs[np.isin(pairs[:, 0], dup), 1]] = True + node_is_patch |= dup_nodes + + # TODO(MEASURE): #629 campaign knob — "slit" keys the strong patch on the + # PHYSICS (the split-duplicated fault nodes), not the refinement + # structure: the smooth refined bulk stays on whole-level smoothing (it + # is well served by ordinary multigrid), and only the fault trace gets + # the subdomain solve. The halo comes from ASM overlap through the + # operator sparsity (set UW_FAC_OVERLAP). Remove when settled. + if os.environ.get("UW_FAC_PATCH") == "slit": + # The fault trace = COINCIDENT FINE PAIRS — the split's plus/minus + # nodes at bit-identical coordinates. Detected in the fine cloud + # itself, not through the transfer: a cut inserts NEW vertices at + # edge crossings, so slit nodes need not coincide with any coarse + # node and the transfer-dup rule (above) can be empty. + slit_nodes = np.zeros(nn, dtype=bool) + cf = np.ascontiguousarray(np.asarray(coords_f, dtype=float)) + _, first, counts = np.unique( + cf.round(decimals=9)[:nn], axis=0, + return_index=True, return_counts=True) + if (counts > 1).any(): + keys = {cf.round(decimals=9)[i].tobytes() + for i in first[counts > 1]} + for i in range(nn): + if cf.round(decimals=9)[i].tobytes() in keys: + slit_nodes[i] = True + slit_nodes |= dup_nodes + rows = np.flatnonzero(slit_nodes[node_f]) + if rows.size == 0: + return None # unsplit level: whole-level smoothing + rows = rows.astype(np.int64) + # TODO(MEASURE): #629 segmentation experiment — split the fault + # subdomain into k along-strike blocks (contiguous chunks of the + # coordinate along the trace's principal direction) to measure the + # iteration cost of block-wise solving vs one trace-wide solve; the + # discriminator for the fault-network coarse-space question. + k = int(os.environ.get("UW_FAC_SEGMENTS", "1")) + if k <= 1: + return rows, rows + pts = coords_f[node_f[rows]] + c0 = pts - pts.mean(axis=0) + # leading principal component = the strike direction + _, _, Vt = np.linalg.svd(c0, full_matrices=False) + s = c0 @ Vt[0] + order = np.argsort(s, kind="stable") + return [(chunk, chunk) for chunk in + (rows[idx] for idx in np.array_split(order, k)) + if chunk.size] + + row_is_patch = node_is_patch[node_f] + # halo: background rows whose coarse node is referenced by a patch row + entry_rows = np.repeat(np.arange(n), np.diff(indptr)) + csel = np.zeros(coords_c.shape[0], dtype=bool) + csel[map_c[cols[row_is_patch[entry_rows]]] // nc] = True + halo_rows = (~row_is_patch) & (bg_col >= 0) & csel[np.clip(bg_col, 0, None)] + node_in_sub = node_is_patch.copy() + node_in_sub[node_f[halo_rows]] = True + + sub = np.flatnonzero(node_in_sub[node_f]) + if sub.size > cover_max * n: + return None + owned = np.flatnonzero(row_is_patch) + return owned.astype(np.int64), sub.astype(np.int64) # --------------------------------------------------------------------------- # @@ -467,7 +759,8 @@ def _build_parallel_transfer(coarse_coords, fine_coords, coarse_layout, Constrained coarse DOFs (``l2g == -1``) drop out -> reduced->reduced.""" l2g_c = coarse_layout.l2g l2g_f, fstart, fend = fine_layout.l2g, fine_layout.rstart, fine_layout.rend - Pn = builder(coarse_coords, fine_coords).tocsr() # (n_f_nodes, n_c_nodes), local + Pn = _drop_structural_zeros( + builder(coarse_coords, fine_coords)) # (n_f_nodes, n_c_nodes), local nloc_f = fend - fstart nloc_c = coarse_layout.rend - coarse_layout.rstart @@ -531,7 +824,7 @@ def _build_crosspart_transfer(coarse_coords, fine_coords, coarse_layout, coords_u, cols_u = _gather_coarse_cloud(coarse_coords, coarse_layout, ncomp, comm) - Pn = builder(coords_u, fine_coords).tocsr() # (n_f_nodes_local, Nu) + Pn = _drop_structural_zeros(builder(coords_u, fine_coords)) # (n_f_local, Nu) nloc_f = fend - fstart nloc_c = coarse_layout.rend - coarse_layout.rstart @@ -595,8 +888,49 @@ def _assert_no_zero_columns_serial(P_csr, level): f"operator would be singular.") +def _repair_zero_columns_serial(P_csr, coords_c, coords_f, map_c, map_f, + nc, level): + """Give each unreached coarse DOF a nearest-fine-DOF entry (serial). + + The barycentric builder has LOCAL support, so on NON-NESTED level pairs + — two independently PLACED meshes (#626), a relaxed child (#424) — a + coarse DOF can lose every fine image and its Galerkin column goes + singular. The dense-RBF rescue fixes that globally at a performance + cliff; when only a handful of columns are empty, a surgical injection + entry (weight 1 at the nearest fine DOF of the same component) makes + the RAP nonsingular at zero cost. These are preconditioner transfers, + not the discretisation — a few injected rows cost iterations at worst, + never correctness. Returns the (possibly repaired) matrix and the + repair count; a column it cannot repair is left empty for the guard + to refuse loudly. + """ + colsum = np.asarray((P_csr != 0).sum(axis=0)).ravel() + zero = np.flatnonzero(colsum == 0) + if not len(zero): + return P_csr, 0 + from scipy.spatial import cKDTree + + n_full_f = coords_f.shape[0] * nc + inv_f = -np.ones(n_full_f, dtype=np.int64) + inv_f[np.asarray(map_f, dtype=np.int64)] = np.arange(len(map_f)) + tree = cKDTree(coords_f) + P = P_csr.tolil() + repaired = 0 + for j in zero: + full_c = int(map_c[j]) + node_c, comp = divmod(full_c, nc) + k = min(8, coords_f.shape[0]) + for i in np.atleast_1d(tree.query(coords_c[node_c], k=k)[1]): + full_f = int(i) * nc + comp + if inv_f[full_f] >= 0: + P[int(inv_f[full_f]), int(j)] = 1.0 + repaired += 1 + break + return P.tocsr(), repaired + + def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None, - ksp=None): + ksp=None, patch_rows=None): """Reconfigure ``pc`` as a fresh PCMG (FMG F-cycle) driven by the supplied reduced->reduced prolongations ``Ps``, Galerkin RAP for coarse operators. @@ -621,12 +955,41 @@ def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None, ``setMGInterpolation`` persists through ``setFromOptions``; the first ``PCSetUp`` builds the coarse operators from our P (no ``MatProductReplaceMats`` shape bug, since the PCMG is fresh and P's size is - fixed).""" + fixed). + + ``patch_rows`` is the FAC configuration (#629): per PCMG level, ``None`` + (whole-level smoothing, the classical setup) or ``(owned, subdomain)`` + reduced-row index arrays from :func:`_fac_patch_split`. A level with a + patch gets its smoother PC switched to ASM with that ONE user subdomain + — the smoother relaxes only the refined patch plus its interface halo, + the background falls through to the level that owns it, and the V-cycle's + smoothing cost follows the patch sizes' geometric series instead of + levels-times-whole-mesh. Residual and transfer work stay global. The + per-level ``mg_levels__pc_type`` keys written here are returned so a + per-solve caller (the rotated path) can drop them from the DB.""" nlev = len(Ps) + 1 prefix = pc.getOptionsPrefix() or "" opts = PETSc.Options() multigrid_options.geometric_mg_bundle(coarse=coarse, smoother=smoother).apply( opts, prefix, owned=owned) + # TODO(MEASURE): #629 contrast-campaign knob — smoother iteration count + # override (the bundle's gmres/4 vs /8 discriminator); remove when settled. + _sm_its = os.environ.get("UW_MG_SMOOTH_ITS") + if _sm_its: + opts[prefix + "mg_levels_ksp_max_it"] = _sm_its + # Per-level override BEFORE setFromOptions: the numbered key beats the + # generic mg_levels_pc_type from the bundle, and having it in the DB keeps + # any later setFromOptions from reverting the live setType below. + fac_keys = [] + _only = os.environ.get("UW_FAC_LEVELS") # TODO(MEASURE): #629 knob + _only = {int(t) for t in _only.split(",")} if _only else None + if patch_rows: + for l in range(1, nlev): + if (l < len(patch_rows) and patch_rows[l] is not None + and (_only is None or l in _only)): + key = f"mg_levels_{l}_pc_type" + opts[prefix + key] = "asm" + fac_keys.append(key) # ``ksp_type`` is in the bundle (#514: a Krylov smoother makes this PC vary # between applications, so its KSP must judge convergence flexibly), but on # the top-level path the KSP consumed its options long before this @@ -643,9 +1006,72 @@ def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None, for l in range(1, nlev): pc.setMGInterpolation(l, Ps[l - 1]) pc.setFromOptions() - - -def _install_transfers(solver, Ps, verbose=False): + # FAC subdomains go on the LIVE smoother PCs — an IS cannot ride the + # options DB. setFromOptions above has already pushed the bundle (and the + # per-level asm keys) into the level KSPs, so the objects are stable now. + if patch_rows: + # TODO(MEASURE): #629 campaign knobs — ASM variant and subdomain + # solver for the FAC smoother; remove once the configuration settles. + # BASIC, not restricted, ASM: the correction must land on the halo + # too. Measured (banded Poisson, 4-level tail): restrict stalls the + # outer KSP at 80-375 iterations where basic runs 6 against a + # whole-level baseline of 4 — discarding the subdomain solve's halo + # correction breaks the interface error systematically. + _asm_type = os.environ.get("UW_FAC_ASM_TYPE", "basic") + _sub_pc = os.environ.get("UW_FAC_SUB_PC") + _whole = os.environ.get("UW_FAC_WHOLE") # asm WITHOUT subdomain + for l in range(1, nlev): + entry = patch_rows[l] if l < len(patch_rows) else None + if entry is None or (_only is not None and l not in _only): + continue + if _whole: + pc.getMGSmoother(l).getPC().setType("asm") + continue + # An entry is one (owned, subdomain) pair, or a LIST of them — + # one ASM subdomain per block (the segmented-fault shape). + blocks = entry if isinstance(entry, list) else [entry] + sm = pc.getMGSmoother(l) + spc = sm.getPC() + spc.setType("asm") + is_subs = [PETSc.IS().createGeneral( + np.asarray(sub_rows, dtype=PETSc.IntType), + comm=PETSc.COMM_SELF) for _own, sub_rows in blocks] + is_owns = [PETSc.IS().createGeneral( + np.asarray(own_rows, dtype=PETSc.IntType), + comm=PETSc.COMM_SELF) for own_rows, _sub in blocks] + # Subdomain = patch + halo; BASIC applies corrections on the + # whole subdomain (restricted stalls — see the class note). + # Overlap 0: the transfer-graph halo is the overlap; the env + # knob adds operator-sparsity layers at PCSetUp instead. + if _asm_type == "basic": + spc.setASMType(PETSc.PC.ASMType.BASIC) + spc.setASMLocalSubdomains(len(blocks), is_subs) + else: + spc.setASMLocalSubdomains(len(blocks), is_subs, is_owns) + # Extra operator-sparsity overlap layers on top of the transfer + # halo (PCASM extends via MatIncreaseOverlap at PCSetUp). + spc.setASMOverlap(int(os.environ.get("UW_FAC_OVERLAP", "0"))) + # Subdomain solver: SOR, the patch-restricted twin of the + # whole-level smoother — no factorization, so no pivot to hit. + # PCASM's default sub-solve (ILU-0) takes a NUMERIC_ZEROPIVOT on + # the rotated Galerkin patch block (measured: min |diag| 8e-5 near + # the constraint; PC_FAILED -11 before the first iteration). + opts[prefix + f"mg_levels_{l}_sub_pc_type"] = _sub_pc or "sor" + fac_keys.append(f"mg_levels_{l}_sub_pc_type") + if (_sub_pc or "sor") in ("lu", "ilu", "cholesky"): + # The rotated Galerkin patch block carries near-zero pivots + # (constraint-zeroed transfer rows leave weakly-attached + # coarse DOFs, min diag ~1e-5); an unshifted factorization + # takes NUMERIC_ZEROPIVOT even as exact LU. + key = f"mg_levels_{l}_sub_pc_factor_shift_type" + opts[prefix + key] = "nonzero" + fac_keys.append(key) + for _is in is_subs + is_owns: # the PC holds its own references + _is.destroy() + return fac_keys + + +def _install_transfers(solver, Ps, verbose=False, patch_rows=None): """Configure the managed PCMG block to use the supplied prolongations. Two paths, keyed by ``solver._pc_option_prefix``: @@ -671,7 +1097,8 @@ def _install_transfers(solver, Ps, verbose=False): ksp.setDMActive(PETSc.KSP.DMActive.OPERATOR, False) _configure_pcmg(ksp.getPC(), Ps, smoother=solver._mg_smoother_variant, - owned=solver._managed_pc_options, ksp=ksp) + owned=solver._managed_pc_options, ksp=ksp, + patch_rows=patch_rows) if verbose: from underworld3 import mpi mpi.pprint(f"[{solver.name}] custom FMG installed: {nlev} levels, " @@ -682,10 +1109,12 @@ def _install_transfers(solver, Ps, verbose=False): raise NotImplementedError( f"custom_mg install: unsupported PC prefix '{pfx}'.") - _install_velocity_block_transfers(solver, Ps, verbose=verbose) + _install_velocity_block_transfers(solver, Ps, verbose=verbose, + patch_rows=patch_rows) -def _install_velocity_block_transfers(solver, Ps, verbose=False): +def _install_velocity_block_transfers(solver, Ps, verbose=False, + patch_rows=None): """Stokes velocity-block install (mechanism A: reset + fresh PCMG). Preconditions: ``solver._build`` + ``setFromOptions`` + ``_attach_stokes_nullspace`` @@ -740,7 +1169,7 @@ def _install_velocity_block_transfers(solver, Ps, verbose=False): vel_pc.reset() vel_pc.setOperators(A_vv, P_vv) _configure_pcmg(vel_pc, Ps, smoother=solver._mg_smoother_variant, - owned=solver._managed_pc_options) + owned=solver._managed_pc_options, patch_rows=patch_rows) vel_pc.setUp() # 4. re-attach the coupled Stokes nullspace (operator state was touched) @@ -777,10 +1206,18 @@ class CustomMGHierarchy: (default) uses the fast path and, if it produces a zero-column transfer (the signature of a cross-partition point-location miss), rebuilds that level cross-partition. Serial builds ignore this. + fac_zone : array_like of bool, or list of such, or None + The fault-zone patch key for the FINEST level's strong (ASM) patch + smoother (#629): a boolean mask over the finest mesh's cells — the + painted / placed band the modeler authored — or a LIST of masks for + one ASM block per fault segment (masks may overlap; a junction cell + in two masks sits inside both blocks). The structural patch (the + cut/split-inserted rows the transfer cannot represent) is unioned in + automatically. ``None`` keys the patch on the structural rows alone. """ def __init__(self, level_meshes, builder="barycentric", field_id=None, - cross_partition="auto"): + cross_partition="auto", fac_zone=None): if builder not in _BUILDERS: raise ValueError("builder must be 'barycentric' or 'rbf'") if len(level_meshes) < 2: @@ -792,8 +1229,33 @@ def __init__(self, level_meshes, builder="barycentric", field_id=None, self.builder_name = builder self.field_id = field_id self.cross_partition = cross_partition + self.fac_zone = self._validated_fac_zone(fac_zone) self.transfers = None + def _validated_fac_zone(self, fac_zone): + """Normalise ``fac_zone`` to None or a list of boolean cell masks. + + Validated against the FINEST level's cell count at construction — + a mask of the wrong length would otherwise surface as a silently + empty (or wrong) patch, which is exactly the declining-quietly + failure mode the ASM probe discipline exists to catch.""" + if fac_zone is None: + return None + masks = (list(fac_zone) if isinstance(fac_zone, (list, tuple)) + else [fac_zone]) + finest = self.level_meshes[-1] + cS, cE = finest.dm.getHeightStratum(0) + out = [] + for k, mk in enumerate(masks): + arr = np.asarray(mk, dtype=bool) + if arr.ndim != 1 or arr.shape[0] != cE - cS: + raise ValueError( + f"fac_zone mask {k}: expected a 1-D boolean mask over the " + f"finest mesh's {cE - cS} cells (plex cell order), got " + f"shape {arr.shape}.") + out.append(arr) + return out + def _recorded_node_transfer(self, level, nlev, degree, n_coarse, n_fine): """The EXACT nested prolongation recorded by ``mesh.adapt`` for this transfer, or ``None`` to fall back to point location. @@ -906,6 +1368,11 @@ def build(self, solver): self._assert_finest_matches_operator(solver, maps[-1], parallel) Ps = [] + # FAC patch smoothing (#629): one entry per PCMG level; level ``l``'s + # patch is read off its own transfer ``Ps[l-1]``. Serial-only for now, + # like the rest of the custom-P specifics; parallel leaves every entry + # None, which _configure_pcmg reads as whole-level smoothing. + self.level_patch_rows = [None] * nlev comm = solver.dm.comm for l in range(1, nlev): if parallel: @@ -924,14 +1391,101 @@ def build(self, solver): _assert_no_zero_columns_parallel(P, comm) Ps.append(P) else: - Pn = self._recorded_node_transfer( - l, nlev, degree, coords[l - 1].shape[0], coords[l].shape[0]) + # A native refine() pair gets the EXACT nested embedding at the + # field's own degree (#629 item 1): the point-located builders + # scatter across the non-nested >P1 node clouds and fatten the + # Galerkin product 3-5x. Declines (None) fall through to the + # recorded/geometric routes unchanged. + Pn = None + # TODO(MEASURE): the env flag is an A/B affordance for the #629 + # benchmark campaign (nested-native vs geometric on the base + # ladder); remove once the comparison is settled. + if (not os.environ.get("UW_CUSTOM_MG_DISABLE_NESTED_NATIVE") + and _is_native_refine_pair(self.level_meshes[l - 1], + self.level_meshes[l])): + Pn = nested_refine_pair_prolongation( + self.level_meshes[l - 1], self.level_meshes[l], + degree, continuous, + coarse_coords=coords[l - 1], fine_coords=coords[l]) + if Pn is None: + Pn = self._recorded_node_transfer( + l, nlev, degree, coords[l - 1].shape[0], + coords[l].shape[0]) if Pn is not None: Pr = _reduced_from_node_transfer(Pn, maps[l - 1], maps[l], nc) else: Pr = _reduced_transfer(coords[l - 1], coords[l], maps[l - 1], maps[l], nc, self.builder) + Pr, n_rep = _repair_zero_columns_serial( + Pr, coords[l - 1], coords[l], maps[l - 1], maps[l], + nc, l) + if n_rep: + import warnings + warnings.warn( + f"custom_mg: transfer {l - 1}->{l} had {n_rep} coarse " + f"DOF(s) with no fine image (non-nested levels); " + f"repaired by nearest-fine-DOF injection. Costs " + f"iterations at worst, never correctness.") _assert_no_zero_columns_serial(Pr, l) + # TODO(MEASURE): A/B affordance for the #629 campaign, like + # the nested-native flag above; remove when settled. + if not os.environ.get("UW_CUSTOM_MG_DISABLE_FAC"): + # An EXPLICIT zone beats detection: a painted weak/TI + # band has no split topology to detect, but the modeler + # knows exactly which cells were painted. The fac_zone + # masks key the finest level's strong patch to those + # cells' DOFs (#629 — the patch solve is rheology- + # agnostic; only detection was split-specific). + if getattr(solver, "_fac_zone_cells", None) is not None: + raise RuntimeError( + "solver._fac_zone_cells is retired; pass the " + "mask(s) as set_custom_fmg(..., fac_zone=...). " + "(A silently ignored zone would leave the patch " + "structural-only — the declining-quietly failure " + "the #629 campaign was bitten by.)") + masks = self.fac_zone + if masks is not None and l == nlev - 1: + # one mask -> one block; a LIST of masks -> one ASM + # block per mask. Masks may overlap (a junction cell + # in two masks sits inside both blocks — the + # junction-coupling design). + cn = self.level_meshes[l]._cell_node_indices( + degree, continuous) + node_of_row = np.asarray(maps[l]) // nc + blocks = [] + covered = np.zeros(len(node_of_row), dtype=bool) + for mk in masks: + node_in = np.zeros(coords[l].shape[0], dtype=bool) + node_in[np.unique( + cn[np.asarray(mk, dtype=bool)])] = True + rows = np.flatnonzero( + node_in[node_of_row]).astype(np.int64) + if rows.size: + blocks.append((rows, rows)) + covered[rows] = True + # The patch smoother REPLACES whole-level smoothing, + # so it inherits every row the coarse level cannot + # represent — the STRUCTURAL (non-identity) transfer + # patch — whether or not the physics zone covers it. + # A zone away from the cut/split leaves those rows + # smoothed nowhere and the velocity solve stalls at + # its cap (measured: any off-cut zone, #629). Add + # the uncovered structural rows as one more block. + split = _fac_patch_split( + Pr, coords[l - 1], coords[l], maps[l - 1], + maps[l], nc, cover_max=1.01) + if split is not None: + _own, sub = split + extra = sub[~covered[sub]].astype(np.int64) + if extra.size: + blocks.append((extra, extra)) + self.level_patch_rows[l] = ( + None if not blocks + else blocks[0] if len(blocks) == 1 else blocks) + else: + self.level_patch_rows[l] = _fac_patch_split( + Pr, coords[l - 1], coords[l], maps[l - 1], + maps[l], nc) Ps.append(_to_petsc_aij(Pr)) self.transfers = Ps return Ps @@ -976,7 +1530,12 @@ def _assert_finest_matches_operator(solver, finest_map, parallel): def install(self, solver, verbose=False): if self.transfers is None: raise RuntimeError("call build() before install()") - _install_transfers(solver, self.transfers, verbose=verbose) + _install_transfers(solver, self.transfers, verbose=verbose, + patch_rows=getattr(self, "level_patch_rows", None)) + # Record what is live on this solver's PC, so a repeat solve can + # skip the re-install (see _pcmg_still_installed). + solver._custom_mg_live = {"h": id(self), + "nlev": len(self.transfers) + 1} class _DMLevelView: @@ -1025,12 +1584,25 @@ def _get_coords_for_basis(self, degree, continuous): from underworld3.discretisation import Mesh return Mesh._get_coords_for_basis(self, degree, continuous) + def _basis_coordinate_dm(self, degree, continuous): + from underworld3.discretisation import Mesh + return Mesh._basis_coordinate_dm(self, degree, continuous) + + def _cell_node_indices(self, degree, continuous): + # Same unbound-delegation pattern as _get_coords_for_basis; the cache + # dict a full Mesh initialises in _setup_ds is created on demand here. + from underworld3.discretisation import Mesh + if not hasattr(self, "_cell_node_array"): + self._cell_node_array = {} + return Mesh._cell_node_indices(self, degree, continuous) + # --------------------------------------------------------------------------- # # Entry points # --------------------------------------------------------------------------- # def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", - field_id=None, cross_partition="auto", verbose=False): + field_id=None, cross_partition="auto", fac_zone=None, + verbose=False): """Generalized custom-P FMG with BC-per-level reduction (the correct path). Registers a :class:`CustomMGHierarchy` on the solver so that the next @@ -1043,12 +1615,21 @@ def set_custom_fmg(solver, coarse_meshes, *, builder="barycentric", ``cross_partition`` selects the parallel (np>1) transfer strategy (see :class:`CustomMGHierarchy`); the default ``"auto"`` handles both nested and - non-nested coarse tails.""" + non-nested coarse tails. + + ``fac_zone`` keys the finest level's strong patch smoother on the fault + zone (#629): a boolean mask over the solver mesh's cells — e.g. + ``mesh.cells_labelled("Band")`` for a placed ribbon, or + ``mesh.cells_supporting("Fault")`` for a split surface's support — or a + list of masks for per-segment ASM blocks. The structural (cut/split) + rows are unioned in automatically; the patch is solved strongly once + per smoother application (basic ASM, shifted factorization).""" solver._custom_mg = { "mode": "hierarchy", "hierarchy": CustomMGHierarchy(list(coarse_meshes) + [solver.mesh], builder=builder, field_id=field_id, - cross_partition=cross_partition), + cross_partition=cross_partition, + fac_zone=fac_zone), "verbose": verbose, } solver.is_setup = False @@ -1151,6 +1732,16 @@ def build_transfers(solver, field_id=None): "remesh between build and solve)") return None, None level_tail = [_DMLevelView(dm, solver.mesh) for dm in hierarchy_dms[:-1]] + # Stamp the refine-family slots so consecutive pairs (including the + # finest pair, whose fine level is the solver's own mesh) take the + # exact nested transfer. Reuse the mesh's token if _coarse_level_meshes + # already stamped one — the level indices coincide by construction. + _slot = getattr(solver.mesh, "_refine_slot", None) + _token = _slot[0] if _slot is not None else object() + for _k, _v in enumerate(level_tail): + _v._refine_slot = (_token, _k) + if _slot is None: + solver.mesh._refine_slot = (_token, len(hierarchy_dms) - 1) builder = "barycentric" else: return None, None # nothing to inject @@ -1283,6 +1874,38 @@ def auto_inject_custom_mg(solver, field_id=None): "auto_cached": True} +def _pcmg_still_installed(solver, h): + """Is THIS hierarchy still live on the solver's managed PC block? + + The marker written by :meth:`CustomMGHierarchy.install` says an install + happened; it cannot say the PC still carries it — a rebuilt SNES, an + explicit ``preconditioner=`` change, or anything that reset the PC leaves + the marker stale. So the marker is only the cheap first test, and the + verdict comes from the LIVE object: the managed block must exist, be a + PCMG, and have the hierarchy's level count. Any doubt (unreachable PC, + un-set-up fieldsplit) answers False — the cost of a wrong False is one + redundant install, the cost of a wrong True is a solve on a stale PC. + """ + mark = getattr(solver, "_custom_mg_live", None) + if not mark or mark.get("h") != id(h): + return False + try: + pfx = solver._pc_option_prefix or "" + if pfx == "": + pc = solver.snes.getKSP().getPC() + elif pfx == "fieldsplit_velocity_": + outer = solver.snes.getKSP().getPC() + if outer.getType() != "fieldsplit": + return False + pc = outer.getFieldSplitSubKSP()[0].getPC() + else: + return False + return (pc.getType() == "mg" + and pc.getMGLevels() == mark.get("nlev")) + except PETSc.Error: + return False + + def inject_custom_mg(solver): """Build + install the custom-P FMG. Called from ``solve()`` (after ``_build``, before the SNES solve) when ``solver._custom_mg`` is set. Dispatches: @@ -1293,6 +1916,16 @@ def inject_custom_mg(solver): if isinstance(cfg, dict) and cfg.get("mode") == "hierarchy": h = cfg["hierarchy"] + if _pcmg_still_installed(solver, h): + # The hierarchy is already live on this solver's PC. PETSc + # re-Galerkins the coarse operators when the fine operator's + # values change at the next PCSetUp, so a repeat solve needs + # NO re-install — and the install is expensive: measured + # (#622) at 55 s of a 61 s warm Stokes solve on an 85k-cell + # cut child, 49.5 s of it a DUPLICATE Jacobian assembly done + # only to make the fieldsplit reachable, 5.8 s rebuilding + # transfers that depend only on the meshes. + return h.build(solver) # parallel-capable (nested co-partitioned) h.install(solver, verbose=cfg.get("verbose", False)) return diff --git a/src/underworld3/utilities/dm_labels.py b/src/underworld3/utilities/dm_labels.py new file mode 100644 index 000000000..6663bf9bf --- /dev/null +++ b/src/underworld3/utilities/dm_labels.py @@ -0,0 +1,44 @@ +"""Empty-safe DMLabel stratum reads (#589). + +``DMLabel.getStratumIS`` hands back a null IS when the stratum is empty on +this rank — the NORMAL case in parallel (a rank owning no part of a fault), +and routine in serial too (a label whose value was never assigned here). +Calling ``getIndices()`` on that null IS is a hard SEGV, and probing +``getStratumIS(v)`` for a value outside the label's live set can abort +outright on some labels (the "Centre" pseudo-label). ``getStratumSize`` is +safe in both situations, so it is the one gate every stratum read goes +through. + +This module is a LEAF — numpy only, no underworld3 imports — so both the +discretisation layer and the utilities can use it without cycles. +""" + +import numpy as np + +__all__ = ["label_stratum_indices"] + + +def label_stratum_indices(label, value): + """The point indices of ``label``'s stratum ``value``, always an array. + + Parameters + ---------- + label : PETSc.DMLabel or None + The label, as returned by ``dm.getLabel(name)`` (which is ``None`` + when the DM has no label of that name — also handled here). + value : int + The stratum value. + + Returns + ------- + numpy.ndarray + The stratum's points (this rank's), dtype int64; EMPTY when the + label is missing or the stratum has no points here. Never raises + for an absent stratum and never touches a null IS. + """ + if label is None or label.getStratumSize(int(value)) == 0: + return np.empty(0, dtype=np.int64) + iss = label.getStratumIS(int(value)) + if iss is None: + return np.empty(0, dtype=np.int64) + return np.asarray(iss.getIndices(), dtype=np.int64) diff --git a/src/underworld3/utilities/edge_split.py b/src/underworld3/utilities/edge_split.py index 877ddb91d..2d1073d21 100644 --- a/src/underworld3/utilities/edge_split.py +++ b/src/underworld3/utilities/edge_split.py @@ -90,15 +90,28 @@ def _cell_edges(dm): return out -def cell_diameters(dm): +def cell_diameters(dm, return_tables=False): """Longest edge length of every cell, in plex cell order. This is the quantity the interpolation error of a linear element depends on, and the one this engine marks against. + + ``return_tables=True`` also returns ``(cell_edges, edge_lengths)`` so + a caller marking AND splitting in the same pass reads the topology + once — the per-cell closure walk is the expensive part (#610), and + the adapt loop was paying it twice per pass. """ L = _edge_lengths(dm) eS, _eE = dm.getDepthStratum(1) - return np.array([L[edges - eS].max() for edges in _cell_edges(dm)]) + edges_of = _cell_edges(dm) + # A simplex has a fixed edge count, so the ragged list stacks and the + # per-cell max vectorises. + E = (np.stack(edges_of) if len(edges_of) + else np.empty((0, 1), dtype=np.int64)) + diam = L[E - eS].max(axis=1) if len(edges_of) else np.zeros(0) + if return_tables: + return diam, (edges_of, L) + return diam def _sf_logical_or(dm, flag): @@ -229,7 +242,7 @@ def _cells_on_edge(dm, edge): return sorted(seen) -def bisect_longest_edges(dm, cells): +def bisect_longest_edges(dm, cells, tables=None): """Split the longest edge of as many of ``cells`` as one pass allows. Parameters @@ -251,13 +264,21 @@ def bisect_longest_edges(dm, cells): ----- Independence caps one pass, so a cell marked here may still exceed the metric afterwards. Re-mark from the returned mesh and call again. + + ``tables`` is the ``(cell_edges, edge_lengths)`` pair from + ``cell_diameters(dm, return_tables=True)`` computed on the SAME dm — + the marking caller already walked the topology, and walking it a + second time here doubled the adapt loop's dominant cost (#610). """ _register_transform() cS, _cE = dm.getHeightStratum(0) eS, _eE = dm.getDepthStratum(1) - L = _edge_lengths(dm) - edges_of = _cell_edges(dm) + if tables is not None: + edges_of, L = tables + else: + L = _edge_lengths(dm) + edges_of = _cell_edges(dm) wanted = {int(edges_of[int(c) - cS][np.argmax(L[edges_of[int(c) - cS] - eS])]) for c in cells} diff --git a/src/underworld3/utilities/fault_split.py b/src/underworld3/utilities/fault_split.py index 3734c4df1..841363546 100644 --- a/src/underworld3/utilities/fault_split.py +++ b/src/underworld3/utilities/fault_split.py @@ -62,6 +62,7 @@ from petsc4py import PETSc import underworld3 as uw +from underworld3.utilities.dm_labels import label_stratum_indices from underworld3.utilities.reconnect import ( _TOPOLOGY_LABELS, _cell_vertices_and_seam, _coords, _copy_labels, _rebuild_point_sf, _shared_points, _write_coordinates) @@ -285,10 +286,10 @@ def _clone_labels(new, dm, clone_map): if values is None: continue for val in values.getIndices(): - points = source.getStratumIS(int(val)) - if points is None: - continue - for p in points.getIndices(): + # Empty-safe (#589): an absent stratum hands back a NULL IS + # wrapper, never None — the old `is None` guard was dead and + # getIndices() on the wrapper segfaults. + for p in label_stratum_indices(source, val): for q in twins.get(int(p), ()): target.setValue(q, int(val)) diff --git a/src/underworld3/utilities/nvb.py b/src/underworld3/utilities/nvb.py index c891d4534..b035e322f 100644 --- a/src/underworld3/utilities/nvb.py +++ b/src/underworld3/utilities/nvb.py @@ -39,6 +39,8 @@ import numpy as np from petsc4py import PETSc +from underworld3.utilities.dm_labels import label_stratum_indices + __all__ = ["NVBMesh", "TaggedBisectionMesh", "write_tagged_state_label"] #: Per-cell refinement-state label consumed by the native 3D driver @@ -521,10 +523,9 @@ def from_dm(cls, dm, boundaries=(), regions=()): for name, value in boundaries: if not dm.hasLabel(name): continue - iset = dm.getStratumIS(name, value) - if iset is None: - continue - for p in iset.getIndices(): + # Empty-safe (#589): an absent stratum hands back a NULL IS + # wrapper, never None — getIndices() on it segfaults. + for p in label_stratum_indices(dm.getLabel(name), value): if eS <= p < eE: # a boundary edge a, b = (int(v - vS) for v in dm.getCone(p)) self.edge_label[_fs(a, b)] = int(value) @@ -534,10 +535,9 @@ def from_dm(cls, dm, boundaries=(), regions=()): for name, value in regions: if not dm.hasLabel(name): continue - iset = dm.getStratumIS(name, value) - if iset is None: - continue - for p in iset.getIndices(): + # Empty-safe (#589): an absent stratum hands back a NULL + # IS wrapper, never None — getIndices() on it segfaults. + for p in label_stratum_indices(dm.getLabel(name), value): if p in cell_of_pt: self.region[cell_of_pt[p]] = int(value) return self @@ -885,10 +885,8 @@ def from_dm(cls, dm, boundaries=(), regions=()): for name, value in boundaries: if not dm.hasLabel(name): continue - iset = dm.getStratumIS(name, value) - if iset is None: - continue - for p in iset.getIndices(): + # Empty-safe (#589), as in the 2-D from_dm above. + for p in label_stratum_indices(dm.getLabel(name), value): if fS <= p < fE: # a labelled facet clos = dm.getTransitiveClosure(p)[0] fverts = tuple(sorted( @@ -900,10 +898,9 @@ def from_dm(cls, dm, boundaries=(), regions=()): for name, value in regions: if not dm.hasLabel(name): continue - iset = dm.getStratumIS(name, value) - if iset is None: - continue - for p in iset.getIndices(): + # Empty-safe (#589): an absent stratum hands back a NULL + # IS wrapper, never None — getIndices() on it segfaults. + for p in label_stratum_indices(dm.getLabel(name), value): if p in cell_of_pt: self.region[cell_of_pt[p]] = int(value) return self diff --git a/src/underworld3/utilities/place_surface.py b/src/underworld3/utilities/place_surface.py index e285dddb7..c656e7869 100644 --- a/src/underworld3/utilities/place_surface.py +++ b/src/underworld3/utilities/place_surface.py @@ -412,20 +412,28 @@ def _inside_polygon(P, q): return inside -def _gmsh_fill_2d(Xall, ring, chain, holes=()): +def _gmsh_fill_2d(Xall, ring, chain, holes=(), size_of=None): """Triangulate the cavity with gmsh: the ring verbatim, the chain embedded. The ring — the cavity boundary, anticlockwise — goes in as a discrete curve carrying its existing segmentation; the surface's chain as a second - discrete curve embedded in the plane surface. A chain end that IS a ring - vertex (a crossing, or an end on the wall) is expressed by having the - chain's elements reference the ring's own node tag — no duplicate node, no - snapping — and a free end is a tip, gmsh's ordinary free-end embed. + discrete curve embedded in the plane surface. ``chain`` is one polyline + of ``Xall`` indices or a sequence of polylines (a trace crossing a + collar piece in several runs), each its own discrete curve. A chain end + that IS a ring vertex (a crossing, or an end on the wall) is expressed + by having the chain's elements reference the ring's own node tag — no + duplicate node, no snapping — and a free end is a tip, gmsh's ordinary + free-end embed. ``holes`` are further closed loops of ``Xall`` indices excluded from the fill — the 2-D thin volume's skin, meshed elsewhere and sewn on — each a discrete curve of its own, verbatim like the ring. + ``size_of`` is an optional callable ``(x, y) -> h`` giving the target + mesh size at a point of the plane — a VARIABLE-resolution fill (a + sheet fine near its shallow tip, coarse at depth). Without it the + fill takes its size from the constrained curves' own segmentation. + Everything is gated, because a fill that looks plausible and is not conforming is worse than a refusal: zero moved nodes, every input segment an edge of the triangulation, triangles present and returned anticlockwise. @@ -437,9 +445,15 @@ def _gmsh_fill_2d(Xall, ring, chain, holes=()): import gmsh ring = [int(v) for v in ring] - chain = [int(v) for v in (chain if chain is not None else [])] + if chain is None or not len(chain): + chains = [] + elif isinstance(chain[0], (list, tuple, np.ndarray)): + chains = [[int(v) for v in c] for c in chain] + else: + chains = [[int(v) for v in chain]] + flat_chain = [v for c in chains for v in c] holes = [[int(v) for v in loop] for loop in holes] - if len(set(chain)) != len(chain): + if len(set(flat_chain)) != len(flat_chain): raise RuntimeError("the surface's chain repeats a vertex") tag_of = {v: i + 1 for i, v in enumerate(ring)} nxt = len(ring) + 1 @@ -452,7 +466,7 @@ def _gmsh_fill_2d(Xall, ring, chain, holes=()): tag_of[v] = nxt hole_nodes.append(v) nxt += 1 - interior = [v for v in chain if v not in tag_of] + interior = [v for v in flat_chain if v not in tag_of] for v in interior: tag_of[v] = nxt nxt += 1 @@ -480,38 +494,45 @@ def discrete_loop(loop_verts): ring_tag = discrete_loop(ring) hole_tags = [discrete_loop(loop) for loop in holes] - line_tag = None - if chain: + chain_tags = [] + added = set(ring) | set(hole_nodes) + for c in chains: line_tag = gmsh.model.addDiscreteEntity(1) - if interior: + own = [v for v in c if v not in added] + added.update(own) + if own: gmsh.model.mesh.addNodes( - 1, line_tag, [tag_of[v] for v in interior], - np.column_stack([Xall[interior], - np.zeros(len(interior))]) + 1, line_tag, [tag_of[v] for v in own], + np.column_stack([Xall[own], np.zeros(len(own))]) .reshape(-1).tolist()) cseg = np.array([[tag_of[a], tag_of[b]] - for a, b in zip(chain[:-1], chain[1:])], + for a, b in zip(c[:-1], c[1:])], dtype=np.int64) gmsh.model.mesh.addElementsByType(line_tag, 1, [], cseg.reshape(-1).tolist()) + chain_tags.append(line_tag) loops = [gmsh.model.geo.addCurveLoop([ring_tag])] loops += [gmsh.model.geo.addCurveLoop([t]) for t in hole_tags] surf = gmsh.model.geo.addPlaneSurface(loops) gmsh.model.geo.synchronize() - if line_tag is not None: - gmsh.model.mesh.embed(1, [line_tag], 2, surf) + if chain_tags: + gmsh.model.mesh.embed(1, chain_tags, 2, surf) # Sizes bracketing what is already there: fine enough to accept the # chain's own spacing, coarse enough not to refine the cavity beyond - # the surviving mesh around it. - constrained = (ring + ring[:1], chain, + # the surviving mesh around it. A size callback, when given, takes + # over the interior grading (the bracket stays as the guard rail). + constrained = (ring + ring[:1], *chains, *[loop + loop[:1] for loop in holes]) lengths = np.concatenate( [np.linalg.norm(np.diff(Xall[c], axis=0), axis=1) for c in constrained if len(c) > 1]) gmsh.option.setNumber("Mesh.MeshSizeMin", 0.5 * float(lengths.min())) gmsh.option.setNumber("Mesh.MeshSizeMax", 2.0 * float(lengths.max())) + if size_of is not None: + gmsh.model.mesh.setSizeCallback( + lambda dim, tag, x, y, z, lc: float(size_of(x, y))) gmsh.model.mesh.generate(2) out_tags, xyz, _ = gmsh.model.mesh.getNodes() @@ -546,7 +567,8 @@ def discrete_loop(loop_verts): for e in ((int(a), int(b)), (int(b), int(c)), (int(c), int(a))): edges.add((min(e), max(e))) wanted = list(zip(ring, ring[1:] + ring[:1])) - wanted += list(zip(chain, chain[1:])) + for c in chains: + wanted += list(zip(c, c[1:])) for loop in holes: wanted += list(zip(loop, loop[1:] + loop[:1])) missing = [(a, b) for a, b in wanted @@ -1351,6 +1373,81 @@ def _interior_face_counts_3d(dm, exclude=()): return counts +def _reach_query(X, reach): + """A reach-aware spatial hash over ``X``: query(lo, hi) -> indices. + + Points are binned on uniform grids, one per octave of their + ``reach`` (the per-point distance beyond which the caller does not + care), so a box query returns every point whose reach-padded + position could intersect ``[lo, hi]`` — and, on a graded mesh, only + those: fine-reach points bin finely and are returned only near the + box, the few coarse-reach points bin coarsely. This is what keeps + sheet-against-mesh sweeps O(near-band) instead of O(domain x sheet) + (#613 — measured 3,217 s of a 3,705 s placement without it). + """ + X = np.asarray(X, dtype=float) + reach = np.broadcast_to(np.asarray(reach, dtype=float), (len(X),)) + if not len(X): + return lambda lo, hi: np.empty(0, dtype=np.int64) + origin = X.min(axis=0) + r_pos = np.maximum(reach, 1e-30) + octave = np.floor(np.log2(r_pos / r_pos.min())).astype(np.int64) + grids = [] + for b in np.unique(octave): + sel = np.flatnonzero(octave == b) + pitch = float(r_pos[sel].max()) + keys = np.floor((X[sel] - origin) / pitch).astype(np.int64) + bins = {} + for i, k in zip(sel, map(tuple, keys)): + bins.setdefault(k, []).append(int(i)) + grids.append((pitch, bins)) + + def query(lo, hi): + out = [] + for pitch, bins in grids: + k_lo = np.floor((np.asarray(lo) - pitch - origin) + / pitch).astype(np.int64) + k_hi = np.floor((np.asarray(hi) + pitch - origin) + / pitch).astype(np.int64) + for kx in range(k_lo[0], k_hi[0] + 1): + for ky in range(k_lo[1], k_hi[1] + 1): + for kz in range(k_lo[2], k_hi[2] + 1): + got = bins.get((kx, ky, kz)) + if got: + out.extend(got) + return np.asarray(out, dtype=np.int64) + + return query + + +def _sheet_distance_within(X, pts, tris, reach): + """Point-to-sheet distance, EXACT wherever it is below ``reach``. + + ``reach`` (scalar or per-point) is the threshold the caller compares + against; a point farther than its reach from every triangle gets a + sentinel above its reach, never the true distance. Every caller + thresholds where it reads this, so decisions match + :func:`_sheet_distance` exactly at O(near-band) cost (#613). + """ + X = np.asarray(X, dtype=float) + tris = np.asarray(tris, dtype=np.int64) + reach = np.broadcast_to(np.asarray(reach, dtype=float), + (len(X),)).copy() + best = 2.0 * np.maximum(reach, 1e-30) + if not len(X) or not len(tris): + return best + query = _reach_query(X, reach) + for t in tris: + A, B, C = pts[t[0]], pts[t[1]], pts[t[2]] + P3 = np.array([A, B, C]) + idx = query(P3.min(axis=0), P3.max(axis=0)) + if not len(idx): + continue + best[idx] = np.minimum( + best[idx], _sheet_distance(X[idx], pts, [t])) + return best + + def _sheet_distance(X, pts, tris): """Distance from each point of ``X`` to a triangulated sheet. @@ -1384,6 +1481,38 @@ def _sheet_distance(X, pts, tris): return best +def _nearest_facet(X, pts, tris): + """``(distance, facet index)`` from each point of ``X`` to a + triangulated sheet — :func:`_sheet_distance` keeping the argmin, for + when the answer is WHICH facet a point lies on, not how far it is.""" + best = np.full(len(X), np.inf) + which = np.zeros(len(X), dtype=np.int64) + for k, t in enumerate(tris): + A, B, C = pts[t[0]], pts[t[1]], pts[t[2]] + ab, ac = B - A, C - A + n = np.cross(ab, ac) + nn = float(n @ n) + rel = X - A + d00, d01, d11 = float(ab @ ab), float(ab @ ac), float(ac @ ac) + d20, d21 = rel @ ab, rel @ ac + denom = d00 * d11 - d01 * d01 + v = (d11 * d20 - d01 * d21) / denom + w = (d00 * d21 - d01 * d20) / denom + inside = (v >= 0.0) & (w >= 0.0) & (v + w <= 1.0) + d_plane = np.abs(rel @ n) / np.sqrt(nn) + d_edge = np.full(len(X), np.inf) + for P, Q in ((A, B), (B, C), (C, A)): + e = Q - P + u = np.clip(((X - P) @ e) / float(e @ e), 0.0, 1.0) + d_edge = np.minimum( + d_edge, np.linalg.norm(X - (P + u[:, None] * e), axis=1)) + d = np.where(inside, d_plane, d_edge) + closer = d < best + best[closer] = d[closer] + which[closer] = k + return best, which + + def _propagate_vertex(dm, chart_values, mpi_op, np_combine): """Reconcile a chart-length per-point array over the point star-forest. @@ -1529,7 +1658,7 @@ def star_of_marked(m): def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, held_cells, h_vertex, on_wall, shared_chart, - open_wall=None): + open_deletable=None, open_near=None): """Victims, dropped tets and the closed cavity shell around the sheet. The same two-part rule as 2-D — vertices within the clearance go, and any @@ -1541,19 +1670,21 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, guarantees (and this function asserts) that the whole region is interior to this rank — the gather's contract. - ``open_wall = (axis, value)`` is the OUTCROP bowl: wall vertices lying - in that plane may be victims (the cap over the cavity is remeshed), and - the shell may open there — the wall faces of dropped cells come back as - ``cap_faces`` for the caller to pre-mesh. Returns + An OUTCROP passes its frame rule as two vertex masks + (:func:`_outcrop_frame_3d`, the same contract as + :func:`_carve_around_volume_3d`): ``open_deletable`` — wall vertices the + carve may take — and ``open_near`` — wall vertices near the outcrop, + where the cavity may open into cap faces; the wall faces of dropped + cells come back as ``cap_faces`` for the caller to pre-mesh. Returns ``(victims, drop_ids, shell, cap_faces)``. """ from underworld3.utilities.edge_split import cell_diameters h_cell = cell_diameters(dm) - d_sheet = _sheet_distance(X, sheet_pts, sheet_tris) - on_open = np.zeros(len(X), dtype=bool) - if open_wall is not None: - on_open = X[:, open_wall[0]] == open_wall[1] + d_sheet = _sheet_distance_within(X, sheet_pts, sheet_tris, + clearance * h_vertex) + on_open = (open_deletable if open_deletable is not None + else np.zeros(len(X), dtype=bool)) held_vertex = np.zeros(len(X), dtype=bool) if held_cells: @@ -1564,18 +1695,25 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, drop = victim[cells].any(axis=1) cen = X[cells].mean(axis=1) + cen_query = _reach_query(cen, h_cell) for t in sheet_tris: A, B, C = sheet_pts[t[0]], sheet_pts[t[1]], sheet_pts[t[2]] n = np.cross(B - A, C - A) n = n / np.linalg.norm(n) diam = max(np.linalg.norm(B - A), np.linalg.norm(C - A), np.linalg.norm(C - B)) - near = np.linalg.norm(cen - (A + B + C) / 3.0, axis=1) < diam + h_cell - if not near.any(): + centre = (A + B + C) / 3.0 + cand = cen_query(centre - diam, centre + diam) + if not len(cand): + continue + near_c = (np.linalg.norm(cen[cand] - centre, axis=1) + < diam + h_cell[cand]) + sub0 = cand[near_c] + if not len(sub0): continue - s = (X[cells[near]] - A) @ n + s = (X[cells[sub0]] - A) @ n straddle = (s.max(axis=1) > 1e-12) & (s.min(axis=1) < -1e-12) - sub = np.flatnonzero(near)[straddle] + sub = sub0[straddle] if not len(sub): continue rel = cen[sub] - A @@ -1594,7 +1732,8 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, # and a facet intersect" on a sheet re-placed into a cleared region). # The centroid rule from the volume carve closes the gap: any cell whose # centre is within reach of the BOUNDED sheet is dropped. - drop |= _sheet_distance(cen, sheet_pts, sheet_tris) < 0.6 * h_cell + drop |= (_sheet_distance_within(cen, sheet_pts, sheet_tris, + 0.6 * h_cell) < 0.6 * h_cell) if held_cells: drop[list(held_cells)] = False @@ -1612,18 +1751,22 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, shell, cap_faces, drop = _closed_shell_3d( dm, X, cells, drop, victim, held_cells, shared_chart, "sheet", - open_wall=open_wall) + open_vertex=open_near) # The straddle rule — or the shell growth — can drop the whole star of # a vertex that is NOT a victim; such a vertex is on no shell face and # would come through the rebuild as an ISOLATED point (global Euler 2, # not 1 — the measured defect class). Every surviving vertex must have - # a surviving cell. + # a surviving cell — except a PROTECTED wall vertex near the outcrop, + # which is a COLLAR node: the cap is re-triangulated through it and the + # fill's tets reference it (the volume carve's measured rule). referenced = np.zeros(len(X), dtype=bool) if (~drop).any(): referenced[cells[~drop].ravel()] = True orphan = ~referenced & ~victim - if orphan[on_wall & ~on_open].any(): + if open_near is not None: + orphan &= ~(orphan & on_wall & ~on_open & open_near) + if (orphan & on_wall & ~on_open).any(): raise RuntimeError( "the sheet's cavity would strand a domain-wall vertex; the sheet " "must be interior, with clearance to spare") @@ -1631,87 +1774,513 @@ def _carve_cavity_3d(dm, X, cells, sheet_pts, sheet_tris, clearance, return np.flatnonzero(victim), np.flatnonzero(drop), shell, cap_faces -def _clip_sheet_to_box(pts, tris, lo, hi): - """Clip a triangulated sheet to the axis-aligned box ``[lo, hi]``. +def _orient_boundary_complex(verts, tris): + """Orient a closed boundary complex so every facet normal points OUT of + the domain. - The specify-long contract (ruling, 2026-08-11): fault surfaces are - defined generously PAST the domain and prep trims them. Each triangle is - Sutherland–Hodgman-clipped against the six half-spaces; cut points snap - exactly onto the plane, and a cut point is computed from its edge's - endpoints in a canonical order so the two triangles sharing the edge - produce the bitwise-identical point (the dedup relies on it). Returns - ``(pts, tris, on_plane)`` with ``on_plane[k] = 2*axis + side`` for a - point lying exactly on a wall plane, else ``-1``. + The gathered complex carries no orientation, so it is propagated by + shared-edge parity within each connected component, and each component's + global sign comes from its enclosed signed volume: the component + enclosing the greatest volume is the outer boundary (normals away from + the domain, signed volume positive), every other component bounds an + interior cavity — a spherical shell's inner surface — whose + domain-outward normals point INTO the cavity, so its signed volume is + made negative. Returns the reoriented ``(nf, 3)`` triangle array. + """ + tris = [list(map(int, t)) for t in np.asarray(tris)] + edge_owner = {} + for k, (a, b, c) in enumerate(tris): + for e in ((a, b), (b, c), (c, a)): + edge_owner.setdefault(frozenset(e), []).append(k) + seen = np.zeros(len(tris), dtype=bool) + comp = np.full(len(tris), -1, dtype=np.int64) + n_comp = 0 + for start in range(len(tris)): + if seen[start]: + continue + stack = [start] + seen[start] = True + comp[start] = n_comp + while stack: + k = stack.pop() + a, b, c = tris[k] + for e in ((a, b), (b, c), (c, a)): + for j in edge_owner[frozenset(e)]: + if j == k or seen[j]: + continue + ja, jb, jc = tris[j] + # Consistent orientation: the shared edge must appear in + # OPPOSITE directions in its two facets. + if e in {(ja, jb), (jb, jc), (jc, ja)}: + tris[j] = [ja, jc, jb] + seen[j] = True + comp[j] = n_comp + stack.append(j) + n_comp += 1 + T = np.asarray(tris, dtype=np.int64) + P = verts[T] + vol6 = np.einsum("ij,ij->i", np.cross(P[:, 0], P[:, 1]), P[:, 2]) + comp_vol = np.array([vol6[comp == c].sum() for c in range(n_comp)]) + outer = int(np.argmax(np.abs(comp_vol))) + for c in range(n_comp): + if (comp_vol[c] > 0.0) != (c == outer): + rows = comp == c + T[rows] = T[rows][:, [0, 2, 1]] + return T + + +def _boundary_signed_distance(X, verts, tris): + """Signed distance to an ORIENTED closed complex: NEGATIVE inside. + + The sign comes from the nearest facet's outward normal; a tie at a + crease is broken by whichever facet is truly nearest, which is + unambiguous for query points off the surface by more than rounding. + """ + d, at = _nearest_facet(X, verts, tris) + T = verts[np.asarray(tris)[at]] + n = np.cross(T[:, 1] - T[:, 0], T[:, 2] - T[:, 0]) + n = n / np.linalg.norm(n, axis=1)[:, None] + side = np.einsum("ij,ij->i", X - T[:, 0], n) + return np.where(side >= 0.0, d, -d) + + +def _concave_crease_depth(verts, tris): + """Per-facet DEPTH of the deepest locally concave crease it touches. + + A crease between two facets of an ORIENTED complex is concave when + either facet's far vertex lies on the OUTWARD side of the other's + plane, and its depth is that offset: zero for a convex or coplanar + boundary (a box, a uniformly-faceted sphere); the SAGITTA MISMATCH — + metres at Earth scale — where a refined-and-snapped curved boundary + changes resolution; the facet size itself on an inner boundary. The + sequential plane clip over-cuts by at most this depth, so the caller + weighs it against the setback rather than refusing on existence. The + noise floor sits at rounding. + """ + tris = np.asarray(tris, dtype=np.int64) + T = verts[tris] + n = np.cross(T[:, 1] - T[:, 0], T[:, 2] - T[:, 0]) + n = n / np.linalg.norm(n, axis=1)[:, None] + depth = np.zeros(len(tris)) + edge_first = {} + for t, tri in enumerate(tris): + opp = {frozenset((int(tri[0]), int(tri[1]))): int(tri[2]), + frozenset((int(tri[1]), int(tri[2]))): int(tri[0]), + frozenset((int(tri[2]), int(tri[0]))): int(tri[1])} + for e, far in opp.items(): + if e not in edge_first: + edge_first[e] = (t, far) + continue + s, far_s = edge_first[e] + ea, eb = (int(v) for v in e) + scale = float(np.linalg.norm(verts[eb] - verts[ea])) + d = max(float((verts[far] - T[s, 0]) @ n[s]), + float((verts[far_s] - T[t, 0]) @ n[t])) + if d > 1e-12 * scale: + depth[t] = max(depth[t], d) + depth[s] = max(depth[s], d) + return depth + + +def _split_safe_triangulation(pts, tris): + """Give every all-rim face an interior vertex, by interior-edge splits. + + A face with ALL THREE vertices on the sheet's rim cannot be + split-node duplicated (its two copies would carry the same vertex + triple), and both the clip's corner polygons and a gmsh planar mesh + of a polygon produce a few however fine the triangulation. Splitting + the face's longest INTERIOR edge at its midpoint bisects both + sharing faces — the midpoint is a new interior vertex, rim edges + (the trace included) are untouched, and no child is worse-shaped + than its parent (a centroid inside a sliver corner face was tried + and drove the fill to 1e-23 cell volumes). + """ + from collections import Counter + + P = [np.asarray(p, dtype=float) for p in pts] + T = [[int(v) for v in t] for t in tris] + while True: + edge_use = Counter() + for a, b, c in T: + for e in ((a, b), (b, c), (c, a)): + edge_use[(e[0], e[1]) if e[0] < e[1] else (e[1], e[0])] += 1 + rim_v = {v for e, k in edge_use.items() if k == 1 for v in e} + bad = next((k for k, t in enumerate(T) + if all(v in rim_v for v in t)), None) + if bad is None: + break + a, b, c = T[bad] + splittable = [(u, w) for u, w in ((a, b), (b, c), (c, a)) + if edge_use[(u, w) if u < w else (w, u)] == 2] + if not splittable: + raise RuntimeError( + "a sheet face has no interior vertex and no interior " + "edge; the sheet is a single sliver and cannot be placed " + "as a splittable fault") + u, w = max(splittable, + key=lambda e: float(np.linalg.norm(P[e[0]] - P[e[1]]))) + m = len(P) + P.append(0.5 * (P[u] + P[w])) + fresh = [] + for row in T: + if u in row and w in row: + fresh.append([m if v == w else v for v in row]) + row[:] = [m if v == u else v for v in row] + T.extend(fresh) + return np.asarray(P, dtype=float), np.asarray(T, dtype=np.int64) + + +def _resample_planar_sheet(pts, tris, size): + """Re-triangulate a PLANAR sheet at the given size, rim geometry kept. + + The authored triangulation is DATA — its spacing is whatever the + source provided — but the embedded fault's resolution must match the + mesh it cuts, or the verbatim embed forces sliver cells around every + mismatched sheet triangle (the ruling, 2026-08-19). The sheet's rim + is compressed to its corners (exact, collinear runs removed — this + is where a too-fine authored rim coarsens), each straight run is + resampled at the target size, and the interior is re-meshed by gmsh + in the sheet's own plane. Cut corners and quality are gmsh's + business afterwards, not the clip's. Non-planar sheets refuse — a + curved surface needs a parametric remesh, which is not built. + + ``size`` is a length, or a callable ``(x, y, z) -> h`` for a + VARIABLE resolution — the middle ground: a fault fine only where it + approaches the surface (where the carve needs its clearance) and + coarsening with depth to match the mesh it cuts, so the cut never + slices a giant cell yet the deep fault is not over-resolved. """ - lo = np.asarray(lo, dtype=float) - hi = np.asarray(hi, dtype=float) pts = np.asarray(pts, dtype=float) - key_of = {} - out_pts = [] - corner_ids = {} - - def intern(p): - k = tuple(float(x) for x in p) - if k not in key_of: - key_of[k] = len(out_pts) - out_pts.append(np.asarray(p, dtype=float)) - return key_of[k] - - def cut(a_id, b_id, A, B, axis, value): - # Canonical order so both triangles sharing the edge agree bitwise. - if (a_id, tuple(A)) > (b_id, tuple(B)): - a_id, b_id, A, B = b_id, a_id, B, A - t = (value - A[axis]) / (B[axis] - A[axis]) - p = A + t * (B - A) - p[axis] = value # exactly on the plane - return p + tris = np.asarray(tris, dtype=np.int64) + centre = pts.mean(axis=0) + rel = pts - centre + _u_svd, _s, vt = np.linalg.svd(rel, full_matrices=False) + normal = vt[2] + span = float(np.linalg.norm(rel, axis=1).max()) + off = rel @ normal + if float(np.abs(off).max()) > 1e-9 * max(span, 1.0): + raise NotImplementedError( + "resampling a non-planar sheet is not built (deviation " + f"{float(np.abs(off).max()):.2e} from the fitted plane); " + "author the sheet at the target size instead.") + e0, e1 = vt[0], vt[1] + + from collections import Counter + edge_use = Counter() + for a, b, c in tris: + for e in ((int(a), int(b)), (int(b), int(c)), (int(c), int(a))): + edge_use[tuple(sorted(e))] += 1 + loops = _skin_loops([e for e, k in edge_use.items() if k == 1], + what="the sheet's rim") + + size_at = size if callable(size) else (lambda x, y, z: float(size)) + + resampled = [] + for loop in loops: + corners = _compress_collinear_loop(pts[np.asarray(loop)]) + ring = [] + n_c = len(corners) + for i in range(n_c): + A, B = corners[i], corners[(i + 1) % n_c] + # Local target along the run: sampled at both ends and the + # middle, subdivided to the finest of them — a straight run + # spanning shallow to deep takes the shallow (fine) budget, + # which errs toward resolution, never toward a giant segment. + probes = [A, 0.5 * (A + B), B] + h_run = min(float(size_at(*p)) for p in probes) + n_seg = max(1, int(np.ceil(np.linalg.norm(B - A) / h_run))) + ring += [A + (k / n_seg) * (B - A) for k in range(n_seg)] + resampled.append(np.asarray(ring)) + + def area2(P3): + q = np.column_stack([(P3 - centre) @ e0, (P3 - centre) @ e1]) + return 0.5 * float(q[:, 0] @ np.roll(q[:, 1], -1) + - q[:, 1] @ np.roll(q[:, 0], -1)) + + order = sorted(range(len(resampled)), + key=lambda k: -abs(area2(resampled[k]))) + P3 = np.vstack([resampled[k] for k in order]) + P2 = np.column_stack([(P3 - centre) @ e0, (P3 - centre) @ e1]) + def oriented(ids, anticlockwise): + q = P2[np.asarray(ids)] + a = 0.5 * float(q[:, 0] @ np.roll(q[:, 1], -1) + - q[:, 1] @ np.roll(q[:, 0], -1)) + return ids if (a > 0.0) == anticlockwise else ids[::-1] + + ring_ids, holes, start = [], [], 0 + for j, k in enumerate(order): + ids = list(range(start, start + len(resampled[k]))) + start += len(resampled[k]) + if j == 0: + ring_ids = oriented(ids, True) + else: + holes.append(oriented(ids, False)) + size_2d = (None if not callable(size) + else (lambda x, y: size_at(*(centre + x * e0 + y * e1)))) + new_tris, extra2 = _gmsh_fill_2d(P2, ring_ids, None, holes=holes, + size_of=size_2d) + lifted = np.vstack([P3, centre + extra2 @ np.vstack([e0, e1])]) \ + if len(extra2) else P3 + return _split_safe_triangulation( + lifted, np.asarray(new_tris, dtype=np.int64)) + + +def _clip_sheet_to_boundary(pts, tris, dom_verts, dom_tris, tol=1e-12, + setback=0.0): + """Clip a triangulated sheet against the mesh's OWN boundary complex. + + ``setback > 0`` clips against the boundary OFFSET INWARD by that + distance (each region plane shifted along its inward normal): the + sheet stops deliberately short of the surface — a BLIND fault, whose + rim is strictly interior and therefore splittable — rather than + outcropping. On a curved boundary the offset is the shifted faceted + planes' intersection, within a sagitta of the true offset surface, + which is the right contract for "stop an element or two below". + + The specify-long contract (ruling, 2026-08-11) on a general boundary: + fault surfaces are defined generously PAST the domain and prep trims + them. Sheet vertices classify by signed distance to the oriented + complex; a mixed triangle is cut sequentially by the planes of the + boundary facets it can cross, keeping the inside — exact where the + crossed boundary is locally CONVEX (a box wall, a sphere's outer + surface), which covers the outcrop cases. A locally concave crossing + over-cuts by at most the crease's DEPTH + (:func:`_concave_crease_depth`), so it is allowed only when the + setback dwarfs that depth (a resolution transition of a snapped + curved boundary, metres under a kilometres-deep blind tip) and + refused otherwise (an inner boundary; ANY concavity at setback zero, + where cut nodes must land on the complex exactly — the true polyline + cut is not built). The cutting planes are the COPLANAR REGIONS' + (:func:`_coplanar_regions`) — a box wall cuts as one plane however it + is faceted — and cut nodes are computed from canonically ordered edge + endpoints and interned by the ``(edge, region)`` identity, so the two + triangles sharing a cut edge produce the SAME node and the clipped + sheet stays conforming; each cut node lands exactly on its region's + plane. Uncut geometry is preserved verbatim. Every kept node is gated + inside the domain on exit — a violation means the crossing was not + locally convex at the sheet's scale and is a refusal, not a tolerance. + Returns ``(pts, tris)``. + """ + pts_in = np.asarray(pts, dtype=float) + tris_in = np.asarray(tris, dtype=np.int64) + setback = float(setback) + oriented = _orient_boundary_complex(dom_verts, dom_tris) + sd = _boundary_signed_distance(pts_in, dom_verts, oriented) + inside = sd < tol - setback + if inside.all(): + return pts_in, tris_in + + DT = dom_verts[oriented] + f_lo = DT.min(axis=1) + f_hi = DT.max(axis=1) + concave_depth = _concave_crease_depth(dom_verts, oriented) + # The cutting planes are the coplanar regions': two coplanar facets of + # one wall must cut with the IDENTICAL plane, or the two triangles + # sharing a cut edge key their cut by different facets and the sheet + # tears along a duplicated node (measured on the box top wall). + region, planes_of = _coplanar_regions(dom_verts, oriented) + + out_pts = [p for p in pts_in] + cut_id = {} + + def cut_on(key, A, B, r): + # Interned by identity, endpoints canonical: the two triangles + # sharing a cut compute the bitwise-identical point, ONE node, and + # the clipped sheet stays conforming. The plane is the region's, + # shifted inward by the setback. + if key not in cut_id: + anchor, nrm = planes_of[int(r)] + oa = float((A - anchor) @ nrm) + setback + ob = float((B - anchor) @ nrm) + setback + p = A + (oa / (oa - ob)) * (B - A) + p -= (float((p - anchor) @ nrm) + setback) * nrm + cut_id[key] = len(out_pts) + out_pts.append(p) + return cut_id[key] + + def edge_key(a, b): + return (a, b) if a < b else (b, a) out_tris = [] - for tri in tris: - # The polygon starts as the triangle, with vertex identity carried - # so shared cut edges intern to the same point. - poly = [(int(v), pts[int(v)].copy()) for v in tri] - for axis in range(3): - for side, value, keep in ((0, lo[axis], 1.0), (1, hi[axis], -1.0)): - if not poly: - break - nxt = [] - for i in range(len(poly)): - (ai, A), (bi, B) = poly[i], poly[(i + 1) % len(poly)] - a_in = keep * (A[axis] - value) >= 0.0 - b_in = keep * (B[axis] - value) >= 0.0 - if a_in: - nxt.append((ai, A)) - if a_in != b_in: - p = cut(ai, bi, A, B, axis, value) - nxt.append((("cut", min(ai, bi), max(ai, bi), - axis, side), p)) - poly = nxt + for tri in tris_in: + flags = inside[tri] + if flags.all(): + out_tris.append([int(v) for v in tri]) + continue + if not flags.any(): + continue + # The facet planes this triangle can cross: bbox overlap, widened + # by the setback — a shifted plane cuts a triangle that never + # comes within bbox reach of the facet itself. + P3 = pts_in[tri] + t_lo = P3.min(axis=0) - tol - setback + t_hi = P3.max(axis=0) + tol + setback + near = np.flatnonzero(((f_lo <= t_hi) & (f_hi >= t_lo)).all(axis=1)) + # The sequential plane clip over-cuts a locally concave crossing + # by at most the crease DEPTH. Weighed against the setback: metres + # of sagitta mismatch on a refined-and-snapped sphere are harmless + # under a kilometres-deep blind tip, while an inner boundary's + # facet-scale concavity refuses at any realistic setback — and at + # setback zero (an outcrop, where cut nodes must land ON the + # complex exactly) any concavity at all refuses, as before. + delta = float(concave_depth[near].max()) if len(near) else 0.0 + if delta > 0.2 * setback + 1e-15: + raise NotImplementedError( + f"the sheet crosses the domain boundary near a locally " + f"concave crease of depth {delta:.2e} (an inner boundary, " + f"or a resolution transition of a snapped curved " + f"boundary); the sequential plane clip can over-cut by " + f"that much, which setback={setback:.2e} does not cover. " + "Raise the setback above ~5x the depth, refine the " + "boundary uniformly under the crossing, or wait for the " + "polyline cut.") + # Each polygon vertex carries the ORIGINAL sheet edges it lies on: + # a cut on a triangle side is re-derived from that side's original + # endpoints and interned by (edge, region), so the neighbouring + # triangle — whose side may be truncated differently — produces the + # bitwise-identical node (a sub-segment interpolation differs by + # rounding, measured as ~1e-17 duplicate pairs on the sphere). A + # cut on a CHORD (a previous cut's trace across the interior, i.e. + # a crease crossing) is triangle-local and keys by its endpoints. + a0, b0, c0 = (int(v) for v in tri) + sides = {a0: {edge_key(a0, b0), edge_key(a0, c0)}, + b0: {edge_key(a0, b0), edge_key(b0, c0)}, + c0: {edge_key(a0, c0), edge_key(b0, c0)}} + poly = [pts_in[v].copy() for v in (a0, b0, c0)] + rows = [a0, b0, c0] + srcs = [sides[a0], sides[b0], sides[c0]] + for r in sorted({int(region[f]) for f in near}): + anchor, nrm = planes_of[r] + offs = [float((q - anchor) @ nrm) + setback for q in poly] + if all(o > -tol for o in offs): # wholly outside (or on) + poly, rows, srcs = [], [], [] + break + if all(o < tol for o in offs): # wholly inside: no cut + continue + new_poly, new_rows, new_srcs = [], [], [] + m = len(poly) + for i in range(m): + j = (i + 1) % m + oi, oj = offs[i], offs[j] + if oi < tol: + new_poly.append(poly[i]) + new_rows.append(rows[i]) + new_srcs.append(srcs[i]) + if (oi < -tol and oj > tol) or (oi > tol and oj < -tol): + common = srcs[i] & srcs[j] + if common: + e, = common + w = cut_on((e[0], e[1], r), pts_in[e[0]], + pts_in[e[1]], r) + src = {e} + else: + lo_r, hi_r = edge_key(rows[i], rows[j]) + w = cut_on(('x', lo_r, hi_r, r), + np.asarray(out_pts[lo_r]), + np.asarray(out_pts[hi_r]), r) + src = set() + new_poly.append(np.asarray(out_pts[w])) + new_rows.append(w) + new_srcs.append(src) + poly, rows, srcs = new_poly, new_rows, new_srcs if len(poly) < 3: continue - ids = [intern(p) for _tag, p in poly] - for k in range(1, len(ids) - 1): - if len({ids[0], ids[k], ids[k + 1]}) == 3: - out_tris.append((ids[0], ids[k], ids[k + 1])) + for k in range(1, len(rows) - 1): + if len({rows[0], rows[k], rows[k + 1]}) == 3: + out_tris.append([rows[0], rows[k], rows[k + 1]]) if not out_tris: raise ValueError("the sheet lies entirely outside the domain") - out_pts = np.array(out_pts) - on_plane = np.full(len(out_pts), -1, dtype=np.int64) - for axis in range(3): - on_plane[out_pts[:, axis] == lo[axis]] = 2 * axis - on_plane[out_pts[:, axis] == hi[axis]] = 2 * axis + 1 - return out_pts, np.array(out_tris, dtype=np.int64), on_plane - - -def _outcrop_chain(pts, tris, on_plane): - """The sheet's boundary polyline on ONE wall plane, ordered. - - Boundary edges of the clipped sheet whose BOTH ends lie on the same - plane form the outcrop. One open chain on one wall is the supported - case; anything else (two walls, a closed loop, several chains) is - refused with the reason — box-edge outcrops are a later phase. - Returns ``(chain_point_ids, wall_code)`` or ``(None, None)``. + used = sorted({v for t in out_tris for v in t}) + remap = {v: i for i, v in enumerate(used)} + new_pts = np.array([out_pts[v] for v in used]) + new_tris = np.array([[remap[v] for v in t] for t in out_tris], + dtype=np.int64) + + new_pts, new_tris = _split_safe_triangulation(new_pts, new_tris) + + sd_out = _boundary_signed_distance(new_pts, dom_verts, oriented) + if (sd_out > 1e-9 - setback).any(): + raise RuntimeError( + "the clip kept a sheet node outside the domain (or inside the " + "setback strip); the boundary is not locally convex at the " + "crossing's scale, or the crossing spans facets the clip did " + "not see. A defect, not a tolerance.") + return new_pts, new_tris + + +def _sheet_boundary_intersection_chain(pts, tris, dom_verts, dom_tris): + """The sheet's intersection polyline with the boundary — the + surface-trace LOCATOR, by contouring the boundary signed distance + over the sheet's own triangulation (marching triangles). + + Robust by construction on ANY boundary — concave, resolution-graded, + snapped — because every crossing point is interpolated on a sheet + EDGE from that edge's two vertex distances: the two triangles + sharing the edge produce the identical point, so the polyline chains + exactly with no tolerance welding (a direct triangle-triangle + intersection was tried first and fragmented into hundreds of + components at crease endpoints). Accuracy is the linear interpolant's + — O(spacing^2 / R), metres at Earth scale — which is the locator + contract: the deliberate blind-fault workflow needs to know WHERE + the sheet would daylight so it can stop short of it and put the + damage region there, not to mesh against it. Returns the ordered + ``(n, 3)`` polyline, or ``None`` when the sheet never reaches the + boundary; refuses when the intersection is not one open chain. + """ + pts = np.asarray(pts, dtype=float) + oriented = _orient_boundary_complex(dom_verts, dom_tris) + sd = _boundary_signed_distance(pts, dom_verts, oriented) + outside = sd >= 0.0 # the tie sits with outside: 0 or 2 + if outside.all() or not outside.any(): + return None # never reaches, or never inside + + cut = {} + + def crossing(a, b): + key = (a, b) if a < b else (b, a) + if key not in cut: + oa, ob = float(sd[key[0]]), float(sd[key[1]]) + w = oa / (oa - ob) + cut[key] = pts[key[0]] + w * (pts[key[1]] - pts[key[0]]) + return key + + adj = {} + for t in np.asarray(tris, dtype=np.int64): + a, b, c = (int(v) for v in t) + crossed = [crossing(u, w) for u, w in ((a, b), (b, c), (c, a)) + if outside[u] != outside[w]] + if not crossed: + continue + k0, k1 = crossed # a consistent tie-break gives 0 or 2 + adj.setdefault(k0, []).append(k1) + adj.setdefault(k1, []).append(k0) + if not adj: + return None + ends = [k for k, ns in adj.items() if len(ns) == 1] + if len(ends) != 2 or any(len(ns) > 2 for ns in adj.values()): + raise NotImplementedError( + "the sheet's intersection with the boundary is not one open " + "chain; multiple or closed traces are not built.") + chain, prev, cur = [ends[0]], None, ends[0] + while cur != ends[1]: + ns = adj[cur] + nxt = ns[0] if ns[0] != prev else ns[1] + chain.append(nxt) + prev, cur = cur, nxt + return np.asarray([cut[k] for k in chain]) + + +def _outcrop_chain(pts, tris, dom_verts, dom_tris, tol=1e-9): + """The clipped sheet's boundary polyline ON the domain boundary, ordered. + + Boundary edges of the clipped sheet whose ends AND midpoint lie on the + domain's boundary complex form the outcrop trace — the membership rule + of :func:`_split_skin_trace` one dimension down, metric but unambiguous: + the clip put cut nodes ON the crossed facets, while an interior rim node + is a sheet spacing away. One open chain is the supported case; multiple + chains or a closed loop are refused with the reason. Returns the ordered + chain of point rows, or ``None``. """ from collections import Counter @@ -1719,18 +2288,19 @@ def _outcrop_chain(pts, tris, on_plane): for a, b, c in tris: for e in ((int(a), int(b)), (int(b), int(c)), (int(c), int(a))): edge_count[tuple(sorted(e))] += 1 - walls = set(int(w) for w in on_plane[on_plane >= 0]) - if not walls: - return None, None - if len(walls) > 1: - raise NotImplementedError( - "the sheet meets more than one domain wall; box-edge outcrops " - "are not built. Clip the sheet to a single wall.") - wall = walls.pop() - chain_edges = [e for e, k in edge_count.items() if k == 1 - and on_plane[e[0]] == wall and on_plane[e[1]] == wall] + boundary_edges = [e for e, k in edge_count.items() if k == 1] + if not boundary_edges: + return None + on_v = _sheet_distance(pts, dom_verts, dom_tris) < tol + both_on = [e for e in boundary_edges if on_v[e[0]] and on_v[e[1]]] + if not both_on: + return None + mids = 0.5 * (pts[[e[0] for e in both_on]] + + pts[[e[1] for e in both_on]]) + mid_on = _sheet_distance(mids, dom_verts, dom_tris) < tol + chain_edges = [e for e, ok in zip(both_on, mid_on) if ok] if not chain_edges: - return None, None + return None adj = {} for a, b in chain_edges: adj.setdefault(a, []).append(b) @@ -1746,11 +2316,11 @@ def _outcrop_chain(pts, tris, on_plane): nxt = ns[0] if ns[0] != prev else ns[1] chain.append(nxt) prev, cur = cur, nxt - return chain, wall + return chain def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, - noun, open_wall=None): + noun, open_vertex=None): """The cavity shell over ``drop``, GROWN at pinch edges until manifold. Shared by every 3-D carve. The union of victim stars around any object @@ -1761,16 +2331,18 @@ def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, original mesh did not). Growing the drop at every non-manifold edge merges the wedges; dropping more cells only enlarges the fill. - ``open_wall = (axis, value)`` lets the cavity OPEN onto one flat wall — - the outcrop bowl. A dropped cell's wall face lying in that plane becomes - a CAP face (returned separately; the caller pre-meshes the cap); any - other wall contact still refuses. The manifold check runs on shell - ∪ cap, which together must close. + ``open_vertex`` (a vertex mask) lets the cavity OPEN onto the boundary + near an outcrop — the bowl. A dropped cell's wall face whose vertices + all carry the mask becomes a CAP face (returned separately; the caller + pre-meshes the cap); any other wall contact still refuses. The mask is + the caller's frame rule: the wall plane for a box-framed sheet, the + band-touched coplanar regions for a general zone. The manifold check + runs on shell ∪ cap, which together must close. Refusals stay per-object worded via ``noun``; a wall or seam contact and a growth that would need a held cell are refusals, not growth. Returns ``(shell, cap_faces, drop)`` with ``drop`` possibly grown; ``cap_faces`` - is empty when ``open_wall`` is None. + is empty when ``open_vertex`` is None. """ from collections import Counter @@ -1805,8 +2377,8 @@ def _closed_shell_3d(dm, X, cells, drop, victim, held_cells, shared_chart, "after the gather; the region marking under-reached. " "A defect, not a configuration error.") verts = face_verts[f] - if open_wall is not None and all( - X[v][open_wall[0]] == open_wall[1] for v in verts): + if open_vertex is not None and all( + open_vertex[v] for v in verts): cap_faces.append((f, verts)) continue raise RuntimeError( @@ -2248,7 +2820,7 @@ def _validity_and_orientation_gates(new, comm): def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1, - clearance=0.6, verbose=False): + clearance=0.6, verbose=False, *, setback=0.0, size=None): """Embed a triangulated sheet in a 3-D mesh by placing its points. The 3-D form of :func:`place_along_lines`: the sheet's points become mesh @@ -2278,9 +2850,15 @@ def place_sheet(dm, points, triangles, label=CUT_LABEL, label_value=1, A 3-D simplex mesh, serial or distributed. **Not modified.** points, triangles : array_like The sheet: ``(N, 3)`` vertices and ``(M, 3)`` triangle indices — the - form :class:`~underworld3.meshing.FaultSurface` carries. Interior to - the domain, non-self-intersecting, at least a cell from any embedded - surface. + form :class:`~underworld3.meshing.FaultSurface` carries. + Non-self-intersecting, at least a cell from any embedded surface. + The sheet may run PAST the domain: it is clipped against the mesh's + own boundary, and a trace left ON the boundary becomes the sheet's + OUTCROP — the trace chain's edges carry ``