Discovered while fixing #291 (PR #318).
Summary
The empty-stratum-IS pattern that segfaulted `Stokes_Constrained` at
np>1 in `_constrain_interior_multipliers_in_section` also exists in five
other call sites, all using an insufficient guard of the form:
```python
if sis is None or sis.handle == 0:
... # skip
```
or the equivalent
```python
facets = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()]
```
Neither of those guards catches the case exposed by #291: on a rank owning
zero points with the given boundary label value, `getStratumIS(bvalue)`
returns a valid non-None PETSc IS with `getSize() == 0`, non-zero
`handle`, and `bool(IS) == True`. The subsequent `sis.getIndices()`
segfaults on that IS.
The empirical evidence and `bool()` behaviour are documented on #291 and
PR #318.
Sites
Grepped from `development` HEAD, all after PR #318 lands they still carry
the insufficient pattern:
`src/underworld3/utilities/boundary_flux.py`
- line 70 — in `_boundary_field_nodes`
- line 106 — in `_boundary_field_reaction` (or similar)
- line 154 — in `_boundary_facet_elements`
`src/underworld3/utilities/rotated_bc.py`
- line 78 — in `_boundary_velocity_nodes`
These fire in:
- `add_nitsche_bc` (via boundary_flux writes)
- `add_rotated_freeslip_bc` (rotated_bc)
- `boundary_flux` — `sigma_nn` recovery, dynamic-topography hand-off
- Any Consistent-Boundary-Flux path (heat flux / Nusselt, traction reactions)
Reason the segfault hasn't been triggered elsewhere yet
Most 2D-mesh partitions actually place at least one point of each boundary
on each rank. The failure requires a rank owning zero points with the
specific label value — easy to hit on a StructuredQuadBox split along one
axis where one axis-aligned wall lives entirely on the other rank
(exactly the #291 shape). Larger meshes and finer decompositions make it
more likely.
Fix
Match PR #318's guard: replace
```python
if sis is None or sis.handle == 0:
...
```
with
```python
if not (sis and sis.getSize() > 0):
...
```
(or the inverted form for the ternary sites). `bool(sis)` short-circuits
the null-handle case; `sis.getSize() > 0` catches the valid-but-empty
case that the truthiness check misses.
Suggested next step
A single small "harden empty-stratum guards" PR that touches these five
sites and adds an np=2 test that specifically exercises a boundary owned
by only one rank (via `Left`/`Right` on a boxed-mesh partition).
Underworld development team with AI support from Claude Code
Discovered while fixing #291 (PR #318).
Summary
The empty-stratum-IS pattern that segfaulted `Stokes_Constrained` at
np>1 in `_constrain_interior_multipliers_in_section` also exists in five
other call sites, all using an insufficient guard of the form:
```python
if sis is None or sis.handle == 0:
... # skip
```
or the equivalent
```python
facets = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()]
```
Neither of those guards catches the case exposed by #291: on a rank owning
zero points with the given boundary label value, `getStratumIS(bvalue)`
returns a valid non-None PETSc IS with `getSize() == 0`, non-zero
`handle`, and `bool(IS) == True`. The subsequent `sis.getIndices()`
segfaults on that IS.
The empirical evidence and `bool()` behaviour are documented on #291 and
PR #318.
Sites
Grepped from `development` HEAD, all after PR #318 lands they still carry
the insufficient pattern:
`src/underworld3/utilities/boundary_flux.py`
`src/underworld3/utilities/rotated_bc.py`
These fire in:
Reason the segfault hasn't been triggered elsewhere yet
Most 2D-mesh partitions actually place at least one point of each boundary
on each rank. The failure requires a rank owning zero points with the
specific label value — easy to hit on a StructuredQuadBox split along one
axis where one axis-aligned wall lives entirely on the other rank
(exactly the #291 shape). Larger meshes and finer decompositions make it
more likely.
Fix
Match PR #318's guard: replace
```python
if sis is None or sis.handle == 0:
...
```
with
```python
if not (sis and sis.getSize() > 0):
...
```
(or the inverted form for the ternary sites). `bool(sis)` short-circuits
the null-handle case; `sis.getSize() > 0` catches the valid-but-empty
case that the truthiness check misses.
Suggested next step
A single small "harden empty-stratum guards" PR that touches these five
sites and adds an np=2 test that specifically exercises a boundary owned
by only one rank (via `Left`/`Right` on a boxed-mesh partition).
Underworld development team with AI support from Claude Code