feat(mr): predict and take details with the boundary-clamped prediction stencil by gouarin · Pull Request #494 · hpc-maths/samurai

gouarin · 2026-07-27T13:33:37Z

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 #493 and contains them, so they should go in first. Its own commits are
the last two: feat(mr): predict and take details with the boundary-clamped stencil and
test: regenerate the demo references the clamped stencil moves.

Description

The prediction stencil now shifts inward near a boundary instead of reading the outer ghosts the
boundary conditions filled. This is the numerical extension of the boundary rewrite, and it is
where the wavelet stops being fed values it cannot reproduce.

compute_detail_op (1d/2d/3d) and prediction_op (both arms) split their interval into the runs
of constant shift that #492's query returns, and pick the coefficient family per run and per
direction. The hoists survive: the coefficient table and the stencil's storage offsets are
rebuilt per run rather than per interval, and the inner loop is untouched, so a run with every
shift zero runs exactly the code it ran before. Outside the domain, and where no shift fits, the
stencil stays centred - that is today's behaviour, and those ghosts are what the boundary
conditions write.

The periodic wrap comes from the mesh, not from the caller. Threading it through the call chain
was tried and abandoned: prediction_fn is a public customization point
(make_MRAdapt(prediction_fn, ...)), so widening its arity is a user-visible API break that
belongs to a later phase. Mesh_base therefore caches its domain's bounding box - a linear scan
that update_ghost_periodic and MPI neighbour discovery each recompute today - and
prediction_period() reads it per interval for a few integer operations.

Two constants become 2r, which is what they always meant: the max_stencil_radius floor (its
comment read "2 is because prediction_stencil_radius=1, if >1 we don't know what to do" - the
answer is 2r) and the two < 2 throws in reconstruction() and transfer(). At r = 1,
which is what everything ships with, nothing changes; at r >= 2 the mesh now provisions the
ghosts a clamped stencil needs, without which the roundtrip test at radius 2 in 3D returns a
wrong value rather than failing loudly.

Cost at r = 1, on finite-volume-advection-2d (max_level 10, interleaved runs): total
runtime +3.6 %, detail computation +8 %, the prediction inside ghost update +9 %,
update_sub_mesh unchanged - one position query per interval and one coefficient table per run,
paid exactly in the two kernels that changed. The query re-does its row lookups for every
interval of a row, and the intervals of a row are visited together, so hoisting them per
transverse index is where that 8-9 % would come back from.

The two s2_3D roundtrip tolerances go from 1e-13 / 5e-13 to 3e-12, with the reason recorded
next to the table: a shifted stencil has a larger l1 mass than a centred one (10/3 against 4/3
in 1D), so the roundoff floor is position dependent, worst where every direction is clamped at
once - up to 2.5^dim, 15.6 in 3D. The identity the test asserts is untouched, both operators
reproducing the polynomial exactly; only its floor moved, and the observed errors (2.8e-13 and
5.7e-13) sit far inside that bound.

Related issue

None. Part of the boundary-machinery rewrite; follows #492 and #493.

How has this been tested?

tests/test_prediction_boundary_reproduction.cpp carries both halves of the acceptance bar:

the point of the change: a polynomial of degree 2r now has no detail anywhere, the
boundary included, and one degree above still has one. It fails on main;
the interior is bit-identical, asserted with EXPECT_EQ on doubles against the centred
formula recomputed in the test with the kernel's own summation order - a different order
gives a different double, which is the whole point. On a non-polynomial field, so the details
compared are large.

The demo references are regenerated, and the second commit says what justifies it rather than
just doing it. On a uniform raster both meshes can be compared on:

case	cells before / after	difference of the two solutions	error against the exact solution
advection_2d, `Tf = 0.01`	35320 / 35380	Linf 4.4e-16	identical to 5 digits
advection_2d, `Tf = 0.35`, disc crossing the outflow boundary	14722 / 15034	Linf 5.7e-05, at the boundary, on values ~2e-4	identical to 5 digits
level-set MRA, `Tf = 0.1`	5104 / 5188	Linf 2.8e-02 on `phi`, `rho` identical	-

The level-set number is the honest one: 2.8e-02 on a field of order 0.35 is not roundoff. It is
the scale of the scheme's own error there - refining by one level moves the same field by
1.3e-02 - and the largest deviation sits at the same point in both comparisons, (0.688, 0.938),
where the scheme is most sensitive, not at the domain boundary.

455/455 C++ tests pass, 101 demo comparisons pass (1 skipped), pre-commit clean.

codacy-production · 2026-07-27T13:35:44Z

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 162 complexity · 33 duplication

Metric Results

Complexity 162

Duplication 33

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.}

…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.

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.

The prediction stencil now shifts inward near a boundary instead of reading the outer ghosts the boundary conditions filled. That is the numerical extension of the boundary rewrite, and it is where the wavelet stops being fed values it cannot reproduce. `compute_detail_op` (1d/2d/3d) and `prediction_op` (both arms) decompose their interval into the runs of constant shift that prediction_shifts.hpp hands back, and select the coefficient family per run and per direction. The hoists survive: the coefficient table and the stencil's storage offsets are rebuilt per run rather than per interval, and the inner loop is untouched, so a run with every shift zero runs exactly the code it ran before. Outside the domain, and where no shift fits, the stencil stays centred: that is today's behaviour, and those ghosts are what the boundary conditions write. The periodic wrap comes from the mesh, not from the caller. Threading it through the call chain was tried and abandoned: `prediction_fn` is a public customization point (`make_MRAdapt(prediction_fn, ...)`), so widening its arity is the user-visible API break that belongs to a later phase. `Mesh_base` therefore caches its domain's bounding box - a linear scan that `update_ghost_periodic` and MPI neighbour discovery each recompute today - and `prediction_period()` reads it per interval for a few integer operations. Two constants become `2r`, which is what they always meant: the `max_stencil_radius` floor (its comment read "2 is because prediction_stencil_radius=1, if >1 we don't know what to do" - the answer is `2r`), and the two `< 2` throws in `reconstruction()` and `transfer()`. At `r = 1`, the value everything ships with, nothing changes; at `r >= 2` the mesh now provisions the ghosts a clamped stencil needs, without which the roundtrip test at radius 2 in 3D reports a wrong value rather than a missing cell. Cost at `r = 1`, on `finite-volume-advection-2d` (max_level 10, interleaved runs): total runtime +3.6%, `detail computation` +8%, the prediction inside `ghost update` +9%, `update_sub_mesh` unchanged - the price of one position query per interval and one coefficient table per run, paid exactly in the two kernels that changed. The row lookups the query performs are repeated for every interval of a row and could be hoisted per transverse index, which is where that 8-9% would come back from. The two `s2_3D` roundtrip tolerances are raised from 1e-13 / 5e-13 to 3e-12, with the reason recorded next to the table: a shifted stencil has a larger l1 mass than a centred one (10/3 against 4/3 in 1D), so the roundoff floor is position dependent and worst where every direction is clamped at once - up to 2.5^dim, 15.6 in 3D. The identity the test asserts is unaffected, both operators reproducing the polynomial exactly; only its floor moved, and the observed errors (2.8e-13 and 5.7e-13) sit far inside what that bound allows. tests/test_prediction_boundary_reproduction.cpp carries both halves of the acceptance bar: - **the point of the change**: a polynomial of degree `2r` now has no detail anywhere, the boundary included, and one degree above still has one. It fails on main; - **the interior is bit-identical**, asserted with `EXPECT_EQ` on doubles against the centred formula computed in the test with the kernel's own summation order - a different order would give a different double, which is the whole point. On a non-polynomial field, so the details compared are large. 455/455 pass, pre-commit clean. Not in this commit, and needing a decision: the demo reference h5 files no longer match. The mismatch is structural, not a tolerance - adaptation decides differently near a boundary, so the mesh differs (27034 against 27094 connectivity rows on advection_2d). Regenerating them blesses whatever the new code emits unless something independent justifies it, which is why the interior bit-identity assertion above exists.

Adaptation decides differently near a boundary now, so the meshes differ and the stored `.h5` no longer match. The mismatch is structural, not a tolerance: `conftest.py` walks every dataset element-wise, and a changed mesh fails on shape - 27034 against 27094 connectivity rows on advection_2d. Regenerating a snapshot blesses whatever the new code emits unless something independent says the new code is right, so here is what was checked first, on a uniform raster that both meshes can be compared on (the cells are axis-aligned squares carrying one value, so rasterising is exact rather than an interpolation): | case | cells before / after | difference of the two solutions | error against the exact solution | |---|---|---|---| | advection_2d, `Tf = 0.01` | 35320 / 35380 | Linf 4.4e-16 | identical to 5 digits | | advection_2d, `Tf = 0.35`, the disc crossing the outflow boundary | 14722 / 15034 | Linf 5.7e-05, at the boundary, on values ~2e-4 | identical to 5 digits | | level-set MRA, `Tf = 0.1` | 5104 / 5188 | Linf 2.8e-02 on `phi`, `rho` identical | - | The first two are the cases where an answer is available: the solutions agree to roundoff, and where the boundary is actually crossed they agree to 6e-05 on values of 2e-4, with the accuracy against the advected disc unchanged. The level-set case is the honest one to report: 2.8e-02 on a field of order 0.35 is not roundoff. It is the scale of the scheme's own error at that resolution - refining by one level (5104 -> 10414 cells) moves the same field by 1.3e-02 - and the largest deviation sits at the same place in both comparisons, (0.688, 0.938), where the scheme is most sensitive, and not at the domain boundary. Two equally valid discretisations of a nonlinear advected interface, in other words, which is what a changed mesh produces. Alongside, in the library rather than on the outputs: the interior is asserted bit-identical (tests/test_prediction_boundary_reproduction.cpp), and a polynomial of degree 2r now has no detail at the boundary either, which is the property the references had no way of showing. 101 demo comparisons pass, 1 skipped.

gouarin force-pushed the bc-prediction-detail branch 2 times, most recently from 4331129 to 2904316 Compare July 27, 2026 15:10

gouarin mentioned this pull request Jul 27, 2026

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

Open

5 tasks

gouarin added 6 commits July 27, 2026 17:23

gouarin force-pushed the bc-prediction-detail branch from 2904316 to 523fb46 Compare July 27, 2026 15:25

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

feat(mr): predict and take details with the boundary-clamped prediction stencil - #494

feat(mr): predict and take details with the boundary-clamped prediction stencil#494
gouarin wants to merge 6 commits into
hpc-maths:mainfrom
gouarin:bc-prediction-detail

gouarin commented Jul 27, 2026

Uh oh!

codacy-production Bot commented Jul 27, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

gouarin commented Jul 27, 2026

Description

Related issue

How has this been tested?

Uh oh!

codacy-production Bot commented Jul 27, 2026

Up to standards ✅

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant