Skip to content

fix(mr): give the coarse levels the inward reach a clamped prediction stencil needs - #493

Open
gouarin wants to merge 4 commits into
hpc-maths:mainfrom
gouarin:bc-mesh-inward-reach
Open

fix(mr): give the coarse levels the inward reach a clamped prediction stencil needs#493
gouarin wants to merge 4 commits into
hpc-maths:mainfrom
gouarin:bc-mesh-inward-reach

Conversation

@gouarin

@gouarin gouarin commented Jul 27, 2026

Copy link
Copy Markdown
Contributor
  • I have installed pre-commit locally and use it to validate my commits.
  • The PR title follows the conventional commits convention.
  • This new PR is documented.
  • This new PR is tested.

Stacked on #492 and contains it, so #492 should go in first. Its own commit is
fix(mr): give the coarse levels the inward reach a clamped stencil needs.

Description

A prediction stencil of radius r reads r cells each side, and the coarse levels carry exactly
that margin. That is enough only while the stencil stays centred. Near a boundary - the domain's
own or a hole's - it cannot: to read only cells the domain has, it shifts inward, and it then
reaches 2r on one side. With r of margin it names cells the mesh does not hold.

Measured, on tests/test_mra.cpp's own mesh: at level 2 the reference mesh held

(j=-1,[1,5)) (j=0,[1,5)) (j=1,[1,5)) (j=2,[-1,5)) (j=3,[-1,5)) (j=4,[-1,5))

ragged, because a level exists only where the projection chain needs it - and the detail at the
parent (0,3), clamped away from the top boundary, asked for (0,1), which is not there. The
centred stencil reads the row above instead, which is an outer ghost the boundary conditions
filled: exactly the read the boundary rewrite exists to remove.

ghost_width is not the quantity at issue - it is already 2r at r = 1. Which cells each
level holds
is.

The margin becomes r outward and 2r inward, and the asymmetry is the content:

  • outward, r is all there is to read. A centred stencil reaches r, and past the domain
    those cells are the outer ghosts the boundary conditions write. Adding more adds cells nobody
    writes and nobody reads - and a cell the mesh holds but nothing writes is worse than one it does
    not hold, because it is read as a value. The symmetric 2r was tried first: it leaves NaN cells
    at min_level - 1 and at the domain's corners, which the second test below catches;
  • across a periodic boundary the band is needed on both sides, because there the cells exist
    and the periodic exchange fills them, and a stencil clamped in one direction reads further along
    the others. This is the one case the new tests catch unaided: without it, the cell at x = 0,
    level 3 names (-1,2) and the mesh has no such cell.

Two ways out were refused because both decide the stencil from what one rank happens to hold,
which makes results depend on the rank count: falling back to the centred stencil where the mesh
lacks the cells, and skipping those positions.

Cost, at the radius everything ships with: the suite goes from 32.9 s to 34.7 s (+5%), and the
reference mesh of an adapted 2D case is unchanged at 48896 cells, of a 3D one unchanged at
249834 - the band adds nothing where the margin was already wide enough. Widening symmetrically
instead costs 6.5x on the suite (215 s), all of it in the roundtrip test at radius 3 in 6
dimensions, where the stencil goes from 7^6 to 13^6 points.

Implementation note. The band is added as a separate contribution with the same shape as the
existing arms rather than folded into them with union_. Folding it in resolves the set expression
at a different level and the intervals land in the wrong level's cell list - which showed up as
cells appearing two cells outside the domain that no expansion could account for.

Related issue

None. Part of the boundary-machinery rewrite; the prerequisite the migration plan did not have,
found by switching the consumers over and watching test_mra throw.

How has this been tested?

tests/test_prediction_inward_reach.cpp, 6 tests, asserting both halves of the requirement in 1D,
2D, 3D, with a periodic direction, and across an adaptation loop:

  • every cell a (possibly clamped) stencil names, at every position where a detail is computed, is
    one the mesh holds;
  • every cell the mesh holds is written by the time the ghost update returns - the field is filled
    with NaN first, so a cell nothing writes is visible.

The periodic case fails without this change. The evidence for the rest is that test_mra, which
throws as soon as the consumer of these stencils is switched over, passes with the same consumer
commit stacked on top of this one: the failing meshes are intermediate ones inside adapt, which a
test outside adapt cannot reach.

Full suite: 452/452 pass (6 of them new), pre-commit clean.

Code of Conduct

  • I agree to follow this project's Code of Conduct

@codacy-production

Copy link
Copy Markdown

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 131 complexity · 22 duplication

Metric Results
Complexity 131
Duplication 22

View in Codacy

NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer
TIP This summary will be updated as you push new changes.

@gouarin

gouarin commented Jul 27, 2026

Copy link
Copy Markdown
Contributor Author

Application-level cost, measured after the fact on finite-volume-advection-2d
(--min-level 2 --max-level 10 --mr-eps 1e-4 --Tf 0.2, interleaved runs):

timer before after
total runtime 4.85 s 4.98 s (+2.2%)
update_sub_mesh 0.371 s 0.475 s (+28%)
ghost update 0.947 s 0.954 s (noise)
detail computation 0.392 s 0.395 s (unchanged)

The whole cost is in mesh construction and nothing else moves - which is consistent with the
cell counts being unchanged: the band adds no cell the solver has to touch, it adds a second
traversal per level.

Restricting the band to the 2r shell around the boundary, which is where the only cells it
can usefully add live, does not recover it (0.482 s -> 0.475 s). That identifies the cause
rather than fixing it: the extra pass walks the whole projected footprint whatever survives the
intersection, exactly as #476's measurements found for predict_bc (32 target positions,
almost the same cost as project_bc's 872). The shell is kept anyway, because it is what the
requirement actually asks for and it adds fewer cells.

Recovering it means folding the band into the existing traversal instead of adding a second
one. The obvious way - one expression with union_ - is what produced the wrong-level
intervals described in the commit message, so it needs samurai's level resolution for mixed
expressions to be understood first. That is its own piece of work, and it is not on the
critical path: the boundary machinery is 17.5 % of runtime in 2D and 18 % in 3D, of which
project_bc + predict_bc alone are 61-79 %, and those are what the next phase deletes.

@gouarin

gouarin commented Jul 27, 2026

Copy link
Copy Markdown
Contributor Author

MPI: invariant verified, and the remaining failure located

Rank independence holds. The CI's own comparison, run locally: finite-volume-advection-2d --Tf 0.05 --nfiles 20 at 1, 2, 3 and 4 ranks gives identical files (all 19 outputs, via
python/compare.py), on this branch and on #494. That is acceptance item 4, and it is the first
time it has been exercised rather than argued: the shift is classified against the global
replicated domain, so no communication is needed and none was added.

linux-mpi-mamba fails only through petsc, and only from 3 ranks. Reproduced locally with
MPI + petsc:

ranks finite-volume-heat --init-sol crenel -pc_type lu --min-level 3 --max-level 8
1 passes
2 passes
3, 4 fails

Failing at 3 but not 2 is the signature of two ranks sharing ghosts without being registered as
neighbours of each other - impossible with a single pair.

What it is not. The reach a rank advertises for neighbour discovery was understated (the
prediction margin counted as r, now 2r), and correcting it - committed separately here so it
does not read as a fix - does not fix the failure. Nor is it the projection row's
off-diagonal reservation (o_nnz = proj_stencil_size / 2, whose comment admits assuming "half
the stencil can be on other processes"): raising it to the full stencil changes nothing for the
failing rows, so that change is not in this PR either, though the assumption is false for a
coarse ghost whose 2^dim children are all remote and is worth its own look.

What it looks like. The failing entries are (455,13367) (456,13368) (457,13369) …: the
column is consistently the row plus one rank's worth of rows. That is a ghost's local-to-global
column mapping pointing at another rank's copy of the same cell - i.e. the case ci.yml
already documents as unsupported next to a disabled 3-rank test: "multiple ranks own the same
ghost cell, which is not supported by the current implementation of the ghost cells"
. The
ownership arbitration itself succeeds locally (no Mismatch for cell) and fails in CI, which
runs with ParMETIS and therefore a different partition - two symptoms of one area.

So the band pushes the mesh into a configuration petsc's ghost ownership does not support. The
fix belongs there and is independent of the boundary rewrite; this PR cannot land before it.

…oth sides to agree

The exchange that gives a ghost its petsc global index is positional: the n-th value received
is the n-th value the neighbour sent. Both sides walked their own version of the sequence:

- the sender pushed a value only for the cells **it** owned, the receiver read one only for the
  cells it believed **the neighbour** owned. Those are two ranks' opinions of ownership, and
  where they differ - which is exactly what `compute_cell_ownership`'s correction passes exist
  to repair, and cannot always - the sequences drift and every later ghost is handed another
  cell's global index;
- the level range came from the sender's own mesh (`min_level..max_level`) on one side and from
  `0..max_level` on the other, and the intersection was built with the operands in the opposite
  order on each side.

The corruption is silent. A shifted index is a valid-looking one in another rank's range, so
nothing complains until petsc refuses a column it was not preallocated for - one rank's worth of
rows away from the row, which is the signature I chased.

Now: one value per shared unknown whether owned or not (UNSET for what the sender does not own),
a level range derived from both meshes, and an operand order derived from the rank numbers. The
sequences cannot drift, and an actual ownership disagreement is reported instead of being turned
into a wrong index.

`has_duplicates`, the check that catches this class of defect, was O(n^2), which is why it could
only live inside an `assert` - and asserts are compiled out of every build CI runs. Sorting a
copy makes it O(n log n), so it is now affordable enough that keeping it out of release builds is
a choice rather than a necessity.

Measured, `finite-volume-heat --init-sol crenel -pc_type lu --min-level 3 --max-level 8`, with
asserts **on** (`-DCMAKE_CXX_FLAGS_RELEASE=-O1`):

| | 1 rank | 2 | 3 | 4 |
|---|---|---|---|---|
| main | ok | ok | **duplicate indices** | ok |
| main + this commit | ok | ok | **duplicate indices** | ok |
| with the boundary rewrite's wider ghost band, before this commit | ok | ok | fails | **fails** |
| with that band, after | ok | ok | fails | **ok** |

So this repairs the case the wider band exposed, and it leaves a **pre-existing** defect standing:
at 3 ranks the mapping already has duplicate global indices on main. That one is invisible in CI
because the guard is an assert, and it deserves its own investigation - the exchange is no longer
where it comes from.
gouarin added 3 commits July 27, 2026 17:25
numeric/prediction_coefficients.hpp answers what the coefficients are once the shift is
known; nothing yet answers what the shift is at a given cell. This adds that query. No
consumer uses it, so nothing changes behaviour.

for_each_prediction_shift_run() walks one interval and hands back the maximal runs over
which the shift is constant, one shift per direction. Three properties it is built around:

- it classifies against the global, replicated domain(level), never against the cells one
  rank happens to hold, which is what makes the answer partition independent with no
  communication. Whether those cells are present locally is a separate halo question;
- it is one query per interval rather than per cell, and the bulk of an interval comes
  back as a single run with every shift zero, so a consumer keeps its hoisted kernel there
  and cannot move an interior value;
- it is exact on a holed domain. An interval passing over the edge of a hole is split,
  because classifying the whole interval by its worst cell would shift cells that need no
  shift.

A periodic direction has no boundary: stepping off the end of one reaches the cells the
periodic exchange fills from, so nothing is clamped there. The caller passes the wrap per
direction - the same quantity update_ghost_periodic shifts by - and 0 where the direction
is not periodic. A hole still clamps in a periodic direction, because only stepping off
the end of the domain wraps.

Tested without a mesh, the answer being a property of the domain alone: the shift at every
distance from a boundary in 1D, 2D and 3D at radius 1 and 2; the decomposition of an
interval crossing a hole, cell for cell; periodicity per direction; and a cross-check
against a slow per-cell implementation over a deliberately awkward domain - two holes, one
of them one cell wide so that no stencil fits across it, one biting into the edge, and a
block hanging off the side - with a guard asserting the comparison actually meets every
shift a radius-1 stencil can take, a cell outside the domain, and a cell that does not fit.
…l needs

A prediction stencil of radius `r` reads `r` cells each side, and the coarse levels carry
exactly that margin. That is enough only while the stencil stays centred. Near a boundary - the
domain's own or a hole's - it cannot: to read only cells the domain has, it shifts inward, and
it then reaches `2r` on one side. With `r` of margin it names cells the mesh does not hold.

Measured, on `tests/test_mra.cpp`'s own mesh: at level 2 the reference mesh held

```
(j=-1,[1,5)) (j=0,[1,5)) (j=1,[1,5)) (j=2,[-1,5)) (j=3,[-1,5)) (j=4,[-1,5))
```

ragged, because a level exists only where the projection chain needs it - and the detail at the
parent `(0,3)`, clamped away from the top boundary, asked for `(0,1)`, which is not there. The
centred stencil reads the row above instead, which is an outer ghost the boundary conditions
filled: exactly the read the boundary rewrite exists to remove.

`ghost_width` is not the quantity at issue - it is already `2r` at `r = 1`. **Which cells each
level holds** is. Both mesh types are fixed, `MRMesh` and the AMR `Mesh`; the AMR one provisions
its prediction ghosts by hand and needed the same band.

The margin becomes **`r` outward and `2r` inward**, and the asymmetry is the content:

- **outward, `r` is all there is to read.** A centred stencil reaches `r`, and past the domain
  those cells are the outer ghosts the boundary conditions write. Adding more adds cells nobody
  writes and nobody reads - and a cell the mesh holds but nothing writes is worse than one it
  does not hold, because it is read as a value. The symmetric `2r` was tried first: it leaves NaN
  cells at `min_level - 1` and at the domain's corners, which the second test below catches;
- **across a periodic boundary the band is needed on both sides**, because there the cells exist
  and the periodic exchange fills them, and a stencil clamped in one direction reads further
  along the others. This is the one case the new tests catch unaided: without it, the cell at
  `x = 0`, level 3 names `(-1,2)` and the mesh has no such cell.

Two ways out were refused because both decide the stencil from what one rank happens to hold,
which makes results depend on the rank count: falling back to the centred stencil where the mesh
lacks the cells, and skipping those positions.

Restricting the band to the `2r` shell around the boundary, where its only useful cells live,
was tried and abandoned. It is an edge: a stencil at distance 0 is shifted by `r` and reads
distance `2r` **inclusive**, so the shell must be `2r + 1` wide; built one cell short, it aborts
eight demos on a missing cell. It bought 1.5 % of `update_sub_mesh`, which is not worth it.

**Cost**, at the radius everything ships with: the suite goes from 32.9 s to 34.0 s (+3%), and
the reference mesh of an adapted 2D case is *unchanged* at 48896 cells, of a 3D one unchanged at
249834 - the band adds nothing where the margin was already wide enough. Widening symmetrically
instead costs **6.5x** on the suite (215 s), all of it in the roundtrip test at radius 3 in 6
dimensions, where the stencil goes from `7^6` to `13^6` points. On
`finite-volume-advection-2d` (max_level 10): total runtime +2.2 %, all of it in
`update_sub_mesh` (+28 %), with `ghost update` and `detail computation` unmoved and the cell
counts unchanged - the cost is one extra traversal per level, not extra cells.

**Implementation note.** The band is added as a *separate contribution* with the same shape as
the existing arms rather than folded into them with `union_`. Folding it in resolves the set
expression at a different level and the intervals land in the wrong level's cell list - which
showed up as cells appearing two cells outside the domain that no expansion could account for.

tests/test_prediction_inward_reach.cpp asserts both halves of the requirement in 1D, 2D, 3D,
with a periodic direction, and across an adaptation loop: every cell a stencil names is one the
mesh holds, and every cell the mesh holds is written by the time the ghost update returns. The
periodic case fails without this change. The evidence for the rest is that test_mra, which throws
as soon as the consumer of these stencils is switched over, passes with the same consumer commit
stacked on top of this one: the failing meshes are intermediate ones inside adapt, which a test
outside adapt cannot reach.

452/452 pass, pre-commit clean.
`ghost_physical_reach()` is what a rank advertises so that others discover it as a neighbour,
and it counted the prediction margin as `r`. The coarse levels now carry `2r` there, so the
advertised reach understated the truth and a rank whose ghosts overlap another's could go
undiscovered - the situation `compute_cell_ownership` calls an ownership mismatch.

**This does not fix the petsc-MPI failure of this branch**, and it is committed separately so
that it does not read as if it did. Measured: the 3-rank `finite-volume-heat` case still fails
with the reach corrected, and its symptom is a matrix column pointing one rank's worth of rows
away from its own (column = row + the local row count, consistently), which is the
local-to-global column mapping of a ghost, not neighbour discovery. The reach is corrected here
because it is wrong on its own terms - the function's own comment says overestimating it only
costs a few extra neighbours in the exchange lists.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant