MPI: clamp localised SubDimension thickness to the rank's local size

[ggorman]

The bug

Thickness._arg_values localises a left/right SubDimension thickness with

tkn = glb_to_loc(root, rtkn-1, side)
tkn = tkn + 1

The +1 is correct when index_glb_to_loc(offset, side) returns a relative local index, which it does whenever the offset falls inside the rank. But when the offset lies beyond the rank — i.e. the rank sits wholly inside the boundary region — the same call saturates and returns top + 1, which is a count already. The +1 is then applied twice and the rank is handed a thickness one greater than the number of points it owns.

Symptom

The generated boundary loop runs one point past the rank's extent. On a 64-point x over 4 ranks (local extent 16) with thickness 20, rank 0 gets x_ltkn = 17:

nbl	rank0 (owns 0–15)	rank1 (16–31)	rank2 (32–47)	rank3 (48–63)
20	x_ltkn=17	x_ltkn=4	x_rtkn=4	x_rtkn=17
8	x_ltkn=8	0	0	x_rtkn=8

In the generated code the 1-D temporaries are allocated for x_size=16, so rank 0's left loop writes r0[0..16] (one past the end) and rank 3's right loop writes r1[-1..15] (one before the start).

The write corrupts heap state rather than faulting at the point of the write, so it typically surfaces later as an allocator abort (double free or corruption) at free() or at MPI finalize — and on some allocator layouts it does not surface at all, silently leaving a wrong boundary region. We hit both flavours on the same code.

The fix

Clamp the localised thickness to the rank's local size. This is a no-op whenever glb_to_loc genuinely returned an index, so only the saturating case changes.

Test

tests/test_mpi.py::TestCodeGeneration::test_subdimension_thickness_localisation (mode=4) asserts the invariant directly — a localised thickness may never exceed the rank's extent. Without the clamp it fails assert 17 <= 16 on ranks 0 and 3.

Full tests/test_mpi.py + tests/test_subdomains.py on the patched tree: 679 passed, 1 xfailed. One test, TestOperatorAdvanced::test_interpolation_at_uforward[1], is deselected because it fails identically at unpatched main (pre-existing, unrelated).

Provenance

Found via a standalone MPI reproducer for a silently-wrong absorbing boundary in a downstream seismic package. With this clamp the out-of-bounds writes are gone; a separate downstream fix is needed for the profile arithmetic itself, which is not a Devito issue.

[ggorman]

For context: this had already been anticipated on Slack, and I think this PR closes out that question rather than raising a new one.

In #development on 2026-05-27, @FabioLuporini asked "why do we allow thicknesses to be overridden by the user?", adding "I think it potentially exposes us to nasty issues…", and then specifically:

Fabio: can thicknesses be negative?
Ed: They cannot be negative.

and on where overrides are handled:

Ed: The MPI domain decomposition handles the user overrides as they are all processed in Thickness._arg_values

That is exactly the function this PR patches. The thread concluded the handling was sound, and it very nearly is — the localisation is correct in the ordinary case. What it misses is the saturating branch of index_glb_to_loc, which returns a count rather than an index when a rank sits wholly inside the boundary region, so the +1 gets applied twice.

Two things worth drawing out from that exchange:

"They cannot be negative" is right, and it corrects our own earlier diagnosis. An internal note of ours had recorded this failure as "the taper loop start goes negative". It never does. The thickness comes out one too large (17 against a local extent of 16), and the loop runs one element past the end. Anyone searching for a negative bound would not have found this.

The concern was raised six days before we hit it in production. Our first double free or corruption on a 16-rank gradient was 2026-06-02. So the instinct in that thread was correct; what was missing was a reproducer tying it to an observable failure. The test added here is that reproducer, reduced to an invariant that needs no operator at all: a localised thickness must never exceed the rank's extent.

Happy to add a comment in _arg_values pointing at the saturating-versus-index distinction if reviewers think it would help the next reader — it is the non-obvious part.

[FabioLuporini]

"and returns a count already"

then it means the conceptual flaw is in glb_to_loc, which btw simply calls into https://github.com/devitocodes/devito/blob/main/devito/data/decomposition.py#L137

that function should always return indices and not counts

so ask the AI to determine where the inconsistency is

[ggorman]

You were right, and the inconsistency is a bit worse than either of us put it. Fix moved to the root; the call-site clamp is gone.

Where it is

index_glb_to_loc(offset, side) had two return conventions:

if side == LEFT:
    rel_ofs = glb_min + abs_ofs - base
    if abs_ofs >= base and abs_ofs <= top:
        return rel_ofs          # a local index, relative to `base`
    elif abs_ofs > top:
        return top + 1          # not an index

top + 1 and the RIGHT branch's glb_max - base + 1 are one past the end of the valid local range, on every subdomain including the first — so they are not indices in any frame. The LEFT form compounds it by mixing frames: top is absolute while rel_ofs is relative to base, so the overshoot grows with the subdomain's offset rather than staying at one.

On a 64-point dimension over 4 ranks (local extent 16), offset 40:

rank	base	top	in-range LEFT	saturating LEFT	last valid index
0	0	15	—	16	15
1	16	31	—	32	15
2	32	47	8	—	15

Both saturating branches now return top - base, the subdomain's last local index. Callers convert to a count themselves, so Thickness._arg_values needs no clamp.

A second inconsistency this exposed

Thickness._arg_values converts two different ways. left/right passes rtkn-1 and adds one; middle passed rtkn and used the result directly as a count. That only worked while saturation returned a count — with the root fix in place, a wholly covered rank's middle thickness came back one short. Both branches now convert identically. This is in the diff, and test_subdimension_middle_thickness_localisation pins it.

Why nothing caught it

test_glb_to_loc_w_side only exercised in-range offsets — the saturating branch had no coverage at all. It now has the saturating cases, plus test_glb_to_loc_w_side_is_always_an_index, which asserts the invariant over every subdomain, both sides and all offsets. Each new test was confirmed to fail without its corresponding fix.

tests/test_mpi.py is 229 passed; tests/test_data.py 169 passed, 2 skipped.

One thing I did not touch

The rel=False behaviour of this two-argument form looks independently inconsistent — e.g. a subdomain [16, 31] covered from the RIGHT returns 15 rather than its absolute last index 31, and the in-range predicates don't line up with the absolute frame either. No caller uses rel=False here, so I've documented the two-argument form as defined for the default relative case and left the rest alone rather than widen the change. Happy to fix it separately if you'd like it defined.

(Correcting myself on the earlier comment: I described the old return as "a count", and said the discrepancy was invisible on rank 0. Neither is quite right — it is an out-of-range bound on every rank, and rank 0 is off by one too. The clamp happened to paper over both.)

[codecov]

Codecov Report

❌ Patch coverage is 54.34783% with 21 lines in your changes missing coverage. Please review.
✅ Project coverage is 83.53%. Comparing base (e401545) to head (1e8cb2f).

Files with missing lines	Patch %	Lines
tests/test_mpi.py	16.00%	21 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #2993      +/-   ##
==========================================
- Coverage   83.54%   83.53%   -0.02%     
==========================================
  Files         257      257              
  Lines       53922    53965      +43     
  Branches     4613     4622       +9     
==========================================
+ Hits        45047    45077      +30     
- Misses       8076     8093      +17     
+ Partials      799      795       -4     

Flag	Coverage Δ
pytest-gpu-aomp-amdgpuX	68.56% <66.66%> (+<0.01%)	⬆️
pytest-gpu-gcc-	78.16% <54.34%> (-0.01%)	⬇️
pytest-gpu-icx-	78.09% <54.34%> (-0.02%)	⬇️
pytest-gpu-nvc-nvidiaX	69.18% <66.66%> (-0.02%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.