Skip to content

Commit 8690d64

Browse files
committed
Merge remote-tracking branch 'origin/development' into feature/placement-local-gather
2 parents 18f0b94 + 8b0d817 commit 8690d64

39 files changed

Lines changed: 1999 additions & 1470 deletions

docs/advanced/curved-boundary-conditions.md

Lines changed: 27 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -71,29 +71,48 @@ stokes.add_nitsche_bc(0.0, "Fault", direction=fault_normal, gamma=10)
7171
at 1e4 gives 0.15%).
7272

7373

74-
### 2. Penalty Free-Slip (Simple but Fragile)
75-
76-
Use the mesh-derived normals directly with a penalty parameter:
74+
### 2. Penalty Free-Slip (Simple, and Only With the Node Normal)
7775

7876
```python
79-
Gamma = mesh.Gamma
77+
n = mesh.boundary_normal("Boundary") # measure-weighted node normal
8078
penalty = 10000
81-
stokes.add_natural_bc(penalty * Gamma.dot(v.sym) * Gamma, "Boundary")
79+
stokes.add_natural_bc(penalty * n.dot(v.sym) * n, "Boundary")
8280
```
8381

82+
**Use `mesh.boundary_normal(boundary)`, not `mesh.Gamma`.** A penalty written
83+
against the per-facet normal asks a node shared by two facets to satisfy two
84+
different constraints, which on a two-component velocity leaves nothing: push
85+
the coefficient up and the boundary freezes. Measured on an annulus with an
86+
exact solution (Kramer et al. 2021), coefficient `1e6`, cell 0.15 → 0.035:
87+
88+
| normal | leak `u·n` | velocity error | surface stress error |
89+
|---|---|---|---|
90+
| `mesh.Gamma` (facet) | 1e-5 | 0.60, flat under refinement | 0.21 → 0.26, growing |
91+
| `mesh.boundary_normal` (node) | 3e-5 | 1.0e-2 → 4.9e-4 | 2.4e-2 → 1.4e-3 |
92+
93+
The facet-normal row does not converge, and the leak cannot see it: at `1e3` that
94+
penalty leaks 3e-2 and gets the surface stress right to 2e-3, while at `1e8` it
95+
leaks 1e-7 and is 26% wrong. This is the classical over-constraint that the
96+
*consistent* normal was introduced to avoid (Engelman, Sani & Gresho 1982), and
97+
`mesh.boundary_normal` is that normal.
98+
8499
**When to use:**
85100
- Quick prototyping where high accuracy isn't critical
86101
- When Nitsche is not yet available for your solver type
87102

88103
**Limitations:**
89104
- Penalty must be tuned: too small → loose constraint, too large → ill-conditioning
90105
- On spherical shells, penalty can become unstable at moderate resolution
91-
- ~25-30% error on elliptical boundaries when using raw facet normals
106+
- Check the answer, not just the leak: a constraint that is satisfied is not
107+
evidence that the solution is right
92108

93109

94-
### 3. Projected Normals (For Curved Boundaries with Penalty)
110+
### 3. Projected Normals (Superseded by `mesh.boundary_normal`)
95111

96-
Project `mesh.Gamma` onto a continuous mesh variable, which interpolates and smooths the normals:
112+
`mesh.boundary_normal(boundary)` assembles the measure-weighted node normal
113+
directly, tracks mesh deformation, and is correct in parallel, so the recipe
114+
below is kept for reference rather than recommended. Project `mesh.Gamma` onto a
115+
continuous mesh variable, which interpolates and smooths the normals:
97116

98117
```python
99118
import sympy

docs/advanced/performance.md

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,71 @@ Tools for identifying bottlenecks:
2222
- Batch operations for better performance
2323
- Parallel scaling
2424

25+
## Choosing MPI ranks for a MatMult-bound Stokes solve (#635)
26+
27+
A large nonlinear Stokes solve is usually **memory-bandwidth bound, not
28+
compute bound**. In a controlled Gadi campaign on a 2-D power-law Barr–Houseman
29+
fault benchmark (P2/P1, cell size 0.01, `n=3`, viscosity contrast 1e5, tolerance
30+
1e-8), roughly **97% of solve time was inside PETSc `MatMult`**, and every
31+
configuration ran the identical numerical workload — four SNES solves, 24 KSP
32+
solves, the same mesh SHA-256.
33+
34+
That is the useful diagnostic: when the solver path is fixed and the runtime still
35+
moves, you are looking at throughput, not convergence.
36+
37+
| queue / launch | ranks | Stokes solve (s) | peak memory | `MatMult` (Mflop/s) |
38+
|---|---:|---:|---:|---:|
39+
| `normal`, default placement | 8 | 4932 | 8.99 GB | 6707 |
40+
| `normal`, explicit core binding | 8 | 5751 | 9.13 GB | 5738 |
41+
| `normal`, binding + no-ML attempt | 8 | 5708 | 8.97 GB | 5782 |
42+
| `normal`, binding + no-ML attempt | 4 | 8964 | 7.31 GB | 3842 |
43+
| `normal`, binding + no-ML attempt | 12 | **3288** | 10.1 GB | **9995** |
44+
| `normalsr` (Sapphire Rapids) | 8 | 3992 | 9.25 GB | 8308 |
45+
46+
What this supports:
47+
48+
- **Runtime tracks `MatMult` throughput**, not nonlinear behaviour. Twelve ranks gave
49+
the best wall time in this sample (56:02); four ranks cut peak memory to 7.3 GB but
50+
took 2:30:48.
51+
- **Fewer ranks buy memory, not speed.** If a job is near a memory ceiling, dropping
52+
ranks is a legitimate trade — but expect the wall time to move roughly inversely.
53+
- **Newer nodes help.** The `normalsr` 8-rank run beat the `normal` 8-rank run, though
54+
it was still slower than 12 ranks on `normal`.
55+
56+
```{warning}
57+
**Explicit core binding did not help here**, and this campaign cannot say whether it
58+
ever does. `--map-by core --bind-to core` was *slower* than default placement at 8
59+
ranks (5751 s vs 4932 s), but the configurations ran concurrently on different nodes,
60+
and a repeat of the same default-placement case varied by 8.6% on its own (4540 s vs
61+
4932 s). Node-to-node variation is the same size as the effect. A controlled same-node
62+
comparison, or a socket-distributed mapping, is needed before drawing any affinity
63+
conclusion — do not read this table as a recommendation against binding.
64+
```
65+
66+
**Harmless Open MPI noise.** Messages of the form
67+
68+
```text
69+
[LOG_CAT_ML] component basesmuma is not available but requested in hierarchy
70+
[LOG_CAT_ML] ml_discover_hierarchy exited with error
71+
```
72+
73+
are **nonfatal**. Setting `OMPI_MCA_coll=^ml` neither suppressed them nor changed the
74+
timings materially.
75+
76+
**Serialise your BLAS.** Every run above set one thread per rank, which is what you
77+
want when MPI already owns the cores:
78+
79+
```bash
80+
export OPENBLAS_NUM_THREADS=1
81+
export OMP_NUM_THREADS=1
82+
export MKL_NUM_THREADS=1
83+
export NUMEXPR_NUM_THREADS=1
84+
```
85+
86+
**Still open.** A recommended Gadi rank placement for bandwidth-bound solves, whether
87+
UW3's examples should ship explicit mapping/binding options, and a small repeatable
88+
PETSc timing benchmark to accompany the setup notes — see #635.
89+
2590
## Related Documentation
2691

2792
- [Developer: Performance Guidelines](../developer/guidelines/performance-optimization.md)

docs/advanced/troubleshooting.md

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,3 +50,55 @@ everywhere.
5050

5151
See the `Stokes.penalty` property docstring for the full description, and the
5252
`CONSTRAINED_FREESLIP_MULTIPLIER` design note for the derivation.
53+
54+
## A long output path can segfault parallel HDF5 output (#645)
55+
56+
**Symptom.** A native segmentation fault inside `mesh.write()` or
57+
`mesh.write_timestep()` on an HPC filesystem, with no Python traceback. The mesh,
58+
labels, fields, solve, MPI size and output calls are all valid, and the same script
59+
succeeds when only the output directory is renamed.
60+
61+
**Cause.** The failure follows the *length of the full generated filename*, not the
62+
validity of the path. UW3 appends its own suffixes to whatever you pass — a mesh write
63+
becomes `output.mesh.00000.h5`, and a P2 velocity variable becomes
64+
`output.mesh.U.00000.h5` — so a descriptive case directory can push the complete path
65+
past what the native PETSc/HDF5/MPI-I/O stack tolerates.
66+
67+
```{warning}
68+
This happens **well below** the advertised limits. In the reported case every path
69+
component was under `NAME_MAX` (255), the filesystem allowed `PATH_MAX=4096`, and PETSc
70+
was configured with `PETSC_MAX_PATH_LEN=4096`. A 286-character filename still crashed.
71+
```
72+
73+
**Reported thresholds.** On Gadi (PETSc 3.25.4, HDF5 1.12.2p, Open MPI 4.1.7, MPI-enabled
74+
h5py 3.16.0, GPFS) a staged test failed at a **259-character** filename, and 252
75+
characters is the shortest length observed to fail on that stack. The exact first unsafe
76+
length was not established, so treat these as observations on one stack rather than a
77+
portable threshold. A macOS stack (PETSc 3.25.0, HDF5 1.14.6, Open MPI 5.0.10) did not
78+
reproduce the original failure, though the full threshold matrix was not repeated there.
79+
80+
**What to do.**
81+
82+
1. Keep the output root and the generated case identifier **compact**, especially on HPC.
83+
2. Put fixed configuration in HDF5 metadata or a metrics file rather than encoding every
84+
parameter in the directory name. A directory named for the two or three parameters
85+
that actually vary across a campaign is enough to tell runs apart.
86+
3. If output segfaults natively, **print the full generated filename first** — including
87+
the UW3 suffixes — before investigating the mesh, the labels or the solver.
88+
89+
Shortening the longest filename from 286 to 189 characters made an otherwise unchanged
90+
eight-rank benchmark pass end to end, and a subsequent 192-rank run at 1/64 resolution
91+
completed every mesh, Stokes, HDF5, XDMF and postprocessing stage.
92+
93+
**One source of path growth has since been removed.** The snapshot backend used to
94+
expand `mesh.name` — which is the full source path for a mesh loaded from file — into
95+
every generated filename, so a 123-character mesh name became a 346-character output
96+
filename on its own. It now writes deterministic `mesh_0000` stems and calls native
97+
PETSc/HDF5 from the artifact directory using short relative names, keeping the original
98+
name in the wrapper metadata. Snapshot writes are therefore no longer a path-length
99+
amplifier.
100+
101+
That fix does not make the advice above unnecessary. It covers the snapshot backend
102+
only: `Mesh.write_timestep()` and anything else that builds its own output filenames
103+
still passes whatever you give it straight through to the native stack, and the
104+
underlying limit — wherever it actually sits — has not moved.

docs/developer/CHANGELOG.md

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,110 @@ This log tracks significant development work at a conceptual level, suitable for
66

77
## 2026 Q3 (July – September)
88

9+
### The Multiplier Was Not the Whole Traction (August 2026)
10+
11+
**`Stokes_Constrained.topography()` now returns the traction the boundary is
12+
actually held with**, and a new `traction()` exposes it directly. The momentum
13+
row carries `λ + r(n·u − g)`, so the bare multiplier is short by the
14+
augmented-Lagrangian share — `r` times the discrete constraint residual. With the
15+
viscosity-weighted default `r = 1e4·μ(x)` that share is a few per cent of the
16+
surface traction on a uniform-viscosity annulus and most of it across a `1e6`
17+
viscosity step, where `λ` alone reads a tenth of the exact SolCx topography and
18+
is anti-correlated with it. `multiplier()` still returns `λ` and now says what it
19+
is not.
20+
21+
The defect survived because the validation scored a **correlation** (0.9999)
22+
between the multiplier and the recovered normal stress. A correlation is
23+
scale-free and cannot see a systematic amplitude deficit, which is precisely what
24+
a missing share of the load is. The new guard,
25+
`tests/test_1063_constrained_traction.py`, scores a relative `l2` against the
26+
exact SolCx surface topography and carries the bare multiplier as its negative
27+
control.
28+
29+
The corrected quantity is the consistent boundary flux: at convergence
30+
`M_Γ(λ + r(n·u − g))` balances the volume residual restricted to the boundary,
31+
which is the CBF nodal load (Zhong, Gurnis & Hulbert 1993). So the multiplier
32+
route and the rotated constraint's `boundary_normal_traction` are the same
33+
computation, and they agree to 3–5% — inside each route's own error against the
34+
exact answer.
35+
36+
Documentation: `docs/advanced/curved-boundary-conditions.md` now writes the
37+
penalty free-slip recipe against `mesh.boundary_normal` rather than `mesh.Gamma`.
38+
A penalty against the per-facet normal over-constrains the shared nodes and does
39+
not converge — measured on an annulus at coefficient `1e6`, the velocity error
40+
stays at 0.60 and the surface-stress error grows from 0.21 to 0.26 as the mesh is
41+
refined, while the leak reads 1e-5 throughout. (underworld3#607, #608, #614)
42+
43+
### A Singular Recovery Mass, Mistaken for a Penalty Defect (August 2026)
44+
45+
**The grad-div penalty default stays off**, but the reason it was held off turned out to
46+
be a defect somewhere else entirely (#633) — so the objection that had blocked it is gone,
47+
and a different one took its place.
48+
49+
With the penalty at 10, the spherical dynamic topography recovered from the rotated
50+
free-slip reaction dropped 28% at *vertices* while the facet-integrated value stayed
51+
correct. The natural reading — that grad-div augmentation corrupts the de-smearing from
52+
reaction loads to pointwise stress — was wrong.
53+
54+
The de-smearing mass for a 3-D **P2 triangular** trace has vertex rows that sum to
55+
**exactly zero**. Those rows annihilate a constant, so solving `M σ = R` amplifies any
56+
perturbation of the nodal load at vertices by O(1) — and, being an instability rather
57+
than a discretisation error, does so independently of mesh resolution. The recovery was
58+
already 7.6% low with no penalty at all; the penalty only made it large enough to fail a
59+
test whose 12% tolerance had been hiding it.
60+
61+
The discrimination needed a case that was curved but not 3-D. A 2-D annulus reproduces
62+
the signature exactly and then parts company under refinement: its error falls ~O(h²)
63+
while the shell's stays flat at ~0.28 over a 3.2× node-count range. The 2-D P2 **line**
64+
mass has positive vertex row sums, which is why 2-D never showed the defect and why
65+
dimension, curvature and the rotated constraint were all red herrings.
66+
67+
- The zero-mean P2 vertex basis was **already known and documented in #414**, which
68+
recorded the same drift-away-under-refinement we re-measured here. Its mechanism is
69+
the sharper one and is adopted: because the vertex basis has zero surface mean, the
70+
vertex reaction carries essentially only the O(h) facet-normal/geometry error, and the
71+
consistent solve *faithfully reconstructs that error* — it is not amplifying noise.
72+
What #633 adds is the separation from the grad-div penalty (which was blamed for it)
73+
and the fix below, which is #414's own unactioned recommendation (2).
74+
- So `mass="auto"` stops asking. On a 3-D P2 trace it now takes the **consistent** solve,
75+
keeps its superconvergent midpoints, and **reconstructs the vertices from them**: the
76+
three midpoints of a facet determine a unique linear function, so a vertex reads its
77+
two adjacent midpoints and subtracts the opposite one, averaged over incident facets.
78+
Worst-node error against the analytic coefficient, over cellSize 0.25 → 0.11:
79+
80+
| | 0.25 | 0.20 | 0.16 | 0.13 | 0.11 |
81+
|---|---|---|---|---|---|
82+
| surface, P1-projected | 0.041 | 0.026 | 0.018 | 0.013 | 0.008 |
83+
| surface, reconstructed | 0.016 | 0.012 | 0.012 | 0.003 | 0.004 |
84+
| CMB, P1-projected | 0.116 | 0.067 | 0.047 | 0.030 | 0.025 |
85+
| CMB, reconstructed | 0.094 | 0.058 | 0.043 | 0.024 | 0.015 |
86+
87+
Better at every resolution on both boundaries, by 1.8x to 4.9x, and converging. The
88+
simpler P1-projected recovery stays available as `mass="p1"` — it is sound, it just
89+
discards the good data along with the bad.
90+
- `FreeSurface` already used the P1-projected recovery in 3-D, so production dynamic
91+
topography was never affected. The exposure was `mass="auto"`.
92+
- The spherical topography test is refined (cellSize 0.25 → 0.13) and its tolerances
93+
tightened from 0.10/0.12 to 0.01/0.05, set from measured discretisation error with
94+
~2× headroom, and now assert every node class rather than the aggregate — the failure
95+
was confined to one class and an aggregate assertion passed straight through it.
96+
- Filed #637: 3-D recovery accepts only P1/P2 triangular traces, so dynamic topography
97+
has exactly one supported discretisation there and cannot be cross-validated. That
98+
blocked the P3/hex arm of this investigation.
99+
- `Stokes.DEFAULT_PENALTY` was flipped to 10 on the #625 evidence and then **reverted**.
100+
Three tier-A/B tests fail at 10 and pass at 0, and the same three run **8.3x slower**
101+
(9.27 s to 76.92 s, warm cache both ways). The #625 win needs an FMG hierarchy; without
102+
`refinement>=1` the velocity block falls back to GAMG, which is where grad-div
103+
augmentation drives the solve into its iteration cap. The default path is the one
104+
without a hierarchy, so the default serves it; set `penalty=10` explicitly where FMG
105+
is available.
106+
- Two of those three failures are not penalty defects. The Nitsche free-slip leak
107+
(1.234e-4 against a 1e-4 bound) is augmentation perturbing a *weakly* imposed
108+
constraint — a strong rotated constraint is untouched. The swarm one exposed #641:
109+
`evaluate` returns −0.4976 for `sqrt((E**2).trace()/2)` at in-domain points near the
110+
lid-corner singularity, at `penalty=0` as well; the penalty merely moved an accumulated
111+
total across zero.
112+
9113
### The Free Surface Reaches the Spherical Shell (July 2026)
10114

11115
**`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same

0 commit comments

Comments
 (0)