From e8312ad96224f6b17a8f5eca46e046cc6dfafc29 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 5 Jul 2026 22:27:04 +1000 Subject: [PATCH] fix(boundary_flux, rotated_bc): guard against empty-but-valid stratum IS (#319) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Same class of latent parallel segfault as #291 (PR #318), at four more sites in `boundary_flux.py` and `rotated_bc.py`: * `boundary_flux.py:70` — in `_boundary_field_nodes` * `boundary_flux.py:106` — in `_boundary_reaction_by_facet` * `boundary_flux.py:154` — in `_boundary_facet_elements` * `rotated_bc.py:78` — in `_boundary_velocity_nodes` Each uses the guard if sis is None or sis.handle == 0: ... which catches None and null-handle IS values, but MISSES the case that crashed `_constrain_interior_multipliers_in_section` in #291: on a rank owning zero points with a given boundary label value, `dm.getLabel("UW_Boundaries").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. These sites fire in: * `add_nitsche_bc` (writes go through boundary_flux), * `add_rotated_freeslip_bc` (rotated_bc), * Consistent-Boundary-Flux traction / heat-flux / Nusselt recovery, * σ_nn / dynamic-topography hand-off. Empirically none has bitten yet because most 2D partitions place at least one point of each boundary on each rank. Fine decompositions or axis-aligned splits on axis-aligned walls (exactly the #291 shape) can still hit it. Fix: replace each guard with if not (sis and sis.getSize() > 0): ... `bool(sis)` short-circuits the null-handle case (project convention); `.getSize() > 0` catches the valid-but-empty case that the None / handle checks miss. test_1018_rotated_freeslip.py: 14/14 pass. test_1019_boundary_flux.py: 1/1 pass. No behaviour change on the happy path. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/boundary_flux.py | 6 +++--- src/underworld3/utilities/rotated_bc.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index df9b81f38..a06ab808d 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -67,7 +67,7 @@ def _boundary_field_nodes(solver, boundary, field_id=0): v0, v1 = dm.getDepthStratum(0) fS, fE = dm.getHeightStratum(1) sis = _boundary_stratum_is(dm, solver.mesh, boundary) - if sis is None or sis.handle == 0: + if not (sis and sis.getSize() > 0): return [], lsec, csec, cvec, v0, v1 facets = [int(z) for z in sis.getIndices()] seen = set(); out = [] @@ -103,7 +103,7 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): if normal is None: # accumulate area-weighted facet normals to the closure nodes sis = _boundary_stratum_is(dm, solver.mesh, boundary) - facets = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] + facets = [] if not (sis and sis.getSize() > 0) else [int(z) for z in sis.getIndices()] fS, fE = dm.getHeightStratum(1) acc = {} for f in facets: @@ -151,7 +151,7 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True): def vcoord(q): return cvec[csec.getOffset(q) // dim] nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)} sis = _boundary_stratum_is(dm, solver.mesh, boundary) - strat = [] if (sis is None or sis.handle == 0) else [int(z) for z in sis.getIndices()] + strat = [] if not (sis and sis.getSize() > 0) else [int(z) for z in sis.getIndices()] local_elems = [] for e in [q for q in strat if e0 <= q < e1]: a, b = (int(c) for c in dm.getCone(e)) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index e6afd60f1..36b6d4229 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -75,7 +75,7 @@ def _boundary_velocity_nodes(solver, boundary, normal=None): # In parallel a rank may own NO part of this boundary → a null IS; guard and return # no local nodes (calling getIndices() on a null IS would segfault). sis = _boundary_stratum_is(dm, solver.mesh, boundary) - if sis is None or sis.handle == 0: + if not (sis and sis.getSize() > 0): return [] facets = [int(z) for z in sis.getIndices()] fS, fE = dm.getHeightStratum(1) # facets (edges in 2D, faces in 3D)