feat: reconstruct, stream and assemble with the boundary-clamped prediction stencil - #500
Open
gouarin wants to merge 30 commits into
Open
feat: reconstruct, stream and assemble with the boundary-clamped prediction stencil#500gouarin wants to merge 30 commits into
gouarin wants to merge 30 commits into
Conversation
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.
The consumers apply the 1D coefficient family as a tensor product, so the cells a shifted stencil reads are a whole box, mixed terms included. Availability read one direction at a time cannot see a cell that is missing only diagonally: at the cell diagonally off the corner of a hole every direction reads cells the domain has, yet the corner of the box is inside the hole, and in the new design nothing fills it. The query now asks the domain about the box. A shift is admissible when the domain holds all of it, and the answer is the most centred admissible one: least shifted overall, then shifting x least (a transverse shift only picks a different row, an x shift moves the innermost loop's reads), then negative, so that the answer never depends on how the domain is stored. On a box domain the two rules agree at every cell, corners included, so this is a statement about re-entrant corners only, whether they belong to a hole or to an L-shaped domain. It also makes fits the joint condition, which is strictly stronger: a plus-shaped domain is three cells wide in each direction and still holds no 3x3 box. Neighbouring runs carrying the same shift are merged, so the runs are maximal in the strict sense the docstring claims and a consumer launches one kernel per genuine change of shift. The slow reference the tests cross-check against was rewritten to the box rule, and it is independent of the machinery: reversing the tie-break in the query alone fails five tests. The 3D sweep now also runs periodic in all three directions, the only case where two transverse directions step off the end of the domain at once.
Same answers, same tests. The query used to be one loop interleaving the row cursors, the availability counts, the periodic top-up, the breakpoint computation and the shift search; it is now named stages: - DomainRow::around() returns how far one row of the periodically extended domain covers around a cell, wrap along the row included, so the driver never mentions the period; - displaced_row() keeps the transverse half of the wrap and enumerates only the directions that can have stepped off the end; - most_centred_fit() picks the shift from the covers alone; - TransverseRows owns both the enumeration of the reachable rows and the index of an offset in it, so the shift table and the row array agree by construction rather than by convention. Small diagrams document the re-entrant corner, the reach and the run decomposition over a hole.
The black-box suite survived the refactor unchanged, which is what it is for; two properties live below what it can observe: - index_of inverts offset, checked at compile time: the shift table and the row array index rows the same way by construction, and a change of enumeration order on one side alone stops compiling; - DomainRow::around reports exactly where its answer stops holding. A cover cut too early only fragments the sweep before the merge glues the pieces back together, so the one-query-per-interval cost rests on a property the public run decomposition cannot show - including that a periodic wrap tops the counts up without moving the breakpoints.
The box-rule search tabulated, for every candidate shift, the indices of the rows its stencil box covers: `(2r+1)^dim x (2r+1)^(dim-1)` entries, which is 16 GB of static storage at radius 3 in six dimensions - the instantiation the prediction roundtrip test makes as soon as a consumer uses the query, and the binary then fails to link. The candidates stay tabulated, on the heap; the rows of a candidate are a few integer operations and are recomputed where they are needed. The per-query arrays of row cursors and row covers had the same shape, `(4r+1)^(dim-1)` entries: five in 2D at radius 1, 371293 in the same six dimensions, which is not a stack object either. They now live on the stack while small and on the heap otherwise. TransverseRows is parameterised by its reach rather than by the radius it is `2r` of, which is what it always was and what a later consumer will need at another reach.
…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.
…s identically The owner of a shared unknown was decided from local heuristics - "the minimum rank of the children *this rank* holds", "the closer of the two gravity centres, pairwise" - whose answers differ from rank to rank, and the disagreements were then negotiated over a bounded number of correction passes. Pairwise closeness is not transitive, so three ranks sharing a coarse cell could each name a different owner and the passes never converged. The wider prediction margin the coarse levels now carry shares more cells and made that failure the norm: `finite-volume-heat` at 3 ranks stopped with "Maximum number of correction passes reached". Ownership is now a function of data every holder has, evaluated identically everywhere, and the exchange only checks that neighbours agree. The principle: a row is assembled by the rank that classifies the cell - real cell, projection ghost, prediction ghost, boundary ghost - in its own mesh, since that classification is what the assembly visits and what guarantees the rank holds every cell the row reads. So the owner is the lowest rank among those that classify the cell, and a cell no rank's real cells give a meaning to goes to the holder whose gravity centre is the closest, ties to the lowest rank. A disagreement is now an error with the rule that fired on each side, not something to repair: it can only mean two ranks hold a common cell without being neighbours. Two more things the same runs turned up: - `Mesh_base::swap` did not swap the gravity centre, so a rank that had adapted kept the centre of its previous mesh while its neighbours held the centre of its current one - the one input to the ownership rule that could differ between ranks, and it did. - the projection rows were preallocated for half their children off-process; the owner is the lowest rank holding one child as a real cell, so every other child may be remote. The check exchange also walks a level range derived from both meshes and an intersection whose operands are ordered by rank, as the numbering exchange already does, since a positional exchange is only as good as the agreement on what is walked. The numbering report goes to stderr so that every rank is heard, the demos silencing stdout on ranks other than 0. Measured with the wider ghost band of the boundary rewrite (`finite-volume-heat --init-sol crenel --min-level 3 --max-level 8`): 1 to 3 ranks pass, with and without load balancing at 1 and 2 ranks, as do the lid-driven cavity at 1 and 2 ranks and the parallel ghost tests. Not fixed here: at 4 ranks, and at 3 with load balancing, petsc still refuses a column of a scheme row (`New nonzero at (335,21735)`), with the same row before and after this change, so it is a preallocation defect of the scheme rows and not an ownership one.
…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.
Not up to standards ⛔🔴 Issues
|
| Category | Results |
|---|---|
| UnusedCode | 5 medium |
🟢 Metrics 403 complexity · 50 duplication
Metric Results Complexity 403 Duplication 50
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
force-pushed
the
bc-prediction-reconstruction
branch
from
September 2, 2026 16:09
fa2917c to
65d7c53
Compare
This was referenced Sep 2, 2026
gouarin
force-pushed
the
bc-prediction-reconstruction
branch
5 times, most recently
from
September 3, 2026 06:06
699babc to
0daac0e
Compare
4 tasks
…iodic Without a periodic direction the two contributions - the r margin, and the 2r band clipped to the domain - are one set: the 2r expansion of the cells, clipped to the domain expanded by r. The two add_prediction_ghosts arms per level become one, with two operands each instead of one and three, and the cell list merges one contribution instead of two. update_sub_mesh goes from 0.229 s to 0.201 s on finite-volume-advection-2d (0.184 s before this PR) and from 3.67 s to 3.06 s on finite-volume-advection-3d (1.75 s before). The periodic case keeps its two arms: there the band extends beyond a periodic boundary and the margin does not, and no single expansion of the domain says both.
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.
os.listdir() hands back the entries in filesystem order, which is not the same on every filesystem; the reference and the generated lists were compared as read, so the same set of files failed the comparison on a filesystem that lists them differently.
…d cavity moved Three references, for three different reasons: - `obstacle_linear_convection` is the one holed domain among the demos. The shift rule now asks whether the whole stencil box lies in the domain rather than each direction on its own, which is the same rule on a box domain and clamps one more cell at a re-entrant corner - exactly the cells around the obstacle. The mesh adapts differently there, so the reference is structurally new (352856 bytes against 500664). - `diff_heated_cavity`'s previous reference held 16 cells against the 3310 the demo produces: it had been regenerated from a run that stopped at its first output. It is now the demo's actual result. - `level_set_from_scratch` differs from its reference by up to 4.5e-3 on about 100 cells around the interface, not at the boundary (4 of 108 within two cells of it), and the previous head of this branch already differed from that reference by the same kind of amount on the same cells while producing bit-identical velocity fields and mesh. The reinitialisation loop of the level set amplifies last-bit differences, and the reference is regenerated here so that it matches the code as it stands; if the comparison turns out to be platform dependent it is this demo's tolerance that needs a look, not the boundary machinery.
The obstacle reference regenerated in the previous commit came from a gcc/aarch64 build, where -ffp-contract=fast is the default: the cell corners of this domain, whose origin is not a multiple of the cell length, round differently under a fused multiply-add, so the point list of the file and its connectivity differ from what the x86 CI runners and clang produce. The other demos have an origin at zero and are unaffected. Regenerated with -DCMAKE_CXX_FLAGS=-ffp-contract=off, which on the code before this PR reproduces the previous reference exactly in mesh and to 2e-12 in the field, and whose connectivity mismatch against the fused build is the one the CI reported (15792 of 15840 entries).
The adapted geometry of the catalog is produced by a multiresolution adaptation, whose mesh near the boundary changes with the clamped detail stencil, so the contact sheet of suite B (geometry x decomposition, coloured by rank) no longer matched its baseline: RMS 22.6 against a tolerance of 15, on the CI as locally. Regenerated at 4 ranks with pytest test_ghost_cases.py --mpl-generate-path=reference/ghost_cases; suite A is unchanged and keeps its baseline.
…he inward reach The clamped prediction stencil reads cells up to 2r inside the domain at the level it predicts from. MRMesh holds them and fills them - that is what hpc-maths#493's inward reach is for - but nothing else does: the from-scratch AMR mesh of level_set_from_scratch holds, with its one-cell ghost layer, coarse cells two in from the boundary under a finer region, because they neighbour a coarse real cell, yet its projection cannot write them, their children not being held. The centred stencil never read them; the clamped one does, and predicts from whatever the allocation left there. The demo's result then differed from one run to the next (3.8e-3 in phi), and from one compiler to the next on the CI, which is how it was found: poisoning every unwritten ghost with NaN puts NaN in 48 real cells after the first time step on this branch, and in none on main, where the same 234 ghosts are unwritten but unread. A mesh now declares that it provides the reach with `static constexpr bool holds_inward_prediction_reach = true;` - MRMesh does - and prediction_domain() records it. Without it every query answers the centred stencil, so a mesh that does not provide the reach keeps reading its outer ghosts exactly as it did before the stencil could shift: level_set_from_scratch is bit for bit main again, and matches its original reference to 2e-16. The test pins both answers on the same boundary cell of an MR and an AMR mesh.
…he AMR meshes are centred again Every reference this PR had regenerated is redone from the state before it, against the demos built with -ffp-contract=off (see the obstacle commit for why): the ones the demos still match keep their original bytes, the others are regenerated. Three come back to what they were, with the AMR meshes reading the centred stencil again: amr-burgers-hat, level-set-amr and level-set-from-scratch, the last of which was the nondeterministic one. Twenty-eight are regenerated - the MR demos the clamped detail and prediction stencils move (advection_2d, level-set-mra, the four heat cases, the lid-driven cavity, the heated cavity, the obstacle, the six adaptive LBM cases, the four reconstruction tutorials), plus the initial-state files that come with them. The heated cavity test now prints the demo's output: it is the only demo with a non-linear block solver, its result has been seen to differ between CI runners, and nothing showed why.
…t prediction uses The matrix rows of the prediction ghosts were still the centred family while the ghost update filled those same ghosts with the clamped one: for an implicit scheme the Jacobian no longer matched the residual near the boundary. The heated cavity showed it on the CI and nowhere else: its Newton iterations went from 2 to 31 between t = 0.3 and t = 0.6, then diverged, the time step collapsed to 5e-5 and the run ended on the initial 16-cell mesh on macos, gcc-11 and clang, while gcc-13 and aarch64 happened to converge to the 3310-cell mesh the reference held. The nonlinear solve of the only demo with a non-linear block solver is what made the inconsistency visible; the linear implicit demos absorbed it as a slightly different operator. assemble_prediction now visits the same shifted nodes as the explicit operator, through prediction_shifts_at() on the coarse level's prediction domain, in a dimension-generic loop that replaces three hand-written copies. This is the assembly half of the consumer switch, brought forward from the reconstruction slice because it belongs with the operators it has to agree with.
…ion rows move The implicit heterogeneous heat, the lid-driven cavity and the heated cavity go through the petsc prediction rows, which now carry the clamped stencil (previous commit). Regenerated from the build without FMA; nothing else moves.
…e node loop In 3D the ghost update predicts through sixteen million runs of two cells, so what the operator does once per run costs as much as what it does per cell. Three things were paid per run and are not any more: - the node loop computed `index + view(stencil) + view(stencil_start)` for each of the 27 nodes, one xtensor sum more than before the stencil could shift; the stencil's origin is now added once per run and the loop is back to the single sum it always had. That was 5.8 s of the 8.0 s the operator spent, measured with per-section timers; 2.5 s of it go away; - the 1D coefficient tables were served from an unordered_map, six lookups per run; they come from a flat array indexed by (shift, parity); - the box query decomposed every interval, when an interval at least `reach` away from every non-periodic face is one centred run - all but one percent of them here. That check comes first and returns at once. finite-volume-advection-3d: ghost update inside the adaptation 6.7 s -> 4.9 s (4.3 s before this PR), the time loop's 3.3 s -> 2.4 s (2.1 s); total 22.9 s -> 20.2 s (17.9 s).
…iction stencil The four consumers of the prediction coefficients now take the same stencil near a boundary. `compute_detail_op` and `prediction_op` already shifted it inward; the composed maps of reconstruction.hpp - `reconstruction()`, `portion()`, `transfer()` and the LBM stream that is built on them - and the prediction rows of the petsc assembly still used the centred family everywhere, so a field predicted one way was reconstructed another, and the matrix a solver assembled was not the operator the explicit path applied. The difficulty is the memo. `prediction(delta_l, indices)` is a mesh-free, global memo of the composed map, and it has to stay one: the LBM stream builds one tap list per (level, component) from it. What the composed map depends on near a boundary is the domain within a fixed reach of the reference cell, and only that, so the memo is now keyed on a **position class** - the cover of every transverse row within `3r + 1` of the cell, capped at that reach - from which the shift rule of prediction_shifts.hpp can be asked at every step of the composition, for any parent and any gap, without a mesh. `3r + 1` is derived, not chosen: the support of the composed map stays within `2r`, but the last step asks the rule about parents within `r` and their candidate boxes reach `2r` further, plus one coarse cell for the children a stream slice places next door. The class is read off the mesh once per run of constant class, by the same box arithmetic or row scan as the shift, and a mesh with no domain to be positioned against - a uniform mesh - is one interior class. The interior is unchanged: the recursion visits the nodes in the order it always did and a centred parent carries the coefficient 1 exactly, so away from every boundary the maps are bit for bit those of the centred recursion, which the test asserts against a copy of it. `reconstruction_op_` and `portion` decompose their coarse row into runs of constant class; the LBM stream builds its tap lists per class and looks one up per run, one per strip in the bulk; the petsc prediction rows visit the same shifted nodes as the explicit operator, in a dimension-generic loop that replaced three hand-written copies. `tests/test_reconstruction_boundary.cpp` states the property the way this whole rewrite does: a field holding the cell averages of a polynomial the operator reproduces, every cell outside the domain poisoned with NaN, the ghost update told it has nothing to do - and the reconstruction is exact everywhere, portion on every child of every boundary cell is exact, one degree above is not, and no NaN appears, so no outer ghost was read. The two ways of reading the class, box and row scan, are held to the same runs cell for cell.
cppcheck (shadowFunction) flags a local named 'row' in box_position_runs and in the portion reconstruction, both of which sit in scopes where the free function row() is visible.
…flat table reconstruction_op_ looked the composed map up once per child of every run - a hash of a key that now carries the position class, into a memo nine times larger than the centred one on a 1D domain, whose every coarse level lies within reach of both boundaries - and then iterated the map's own hash table. On burgers_mra (hat, levels 2 to 12) that cost 1.82 s of reconstruction against 1.42 s on main, and a run 12 % slower overall while every timed section had become cheaper. The maps of all the children of one (class, delta_l) are now built once, from the memo, into one contiguous table indexed by the child's x-major linear index, each map flattened into a sorted vector of (offset, weight) terms. The operator pays one lookup per run and iterates contiguous terms. The sums are unchanged term for term; the sort fixes their order so the accumulation is deterministic. burgers_mra goes to 2.79 s against 2.86 s on main.
…ward reach The position-class query follows the shift query: on a mesh that does not declare prediction_inward_reach (see hpc-maths#494), prediction_domain() carries clamp = false and every run is the interior class, so reconstruction(), portion(), transfer(), the LBM stream and the petsc prediction rows use the centred maps there, as they did before the stencil could shift.
…moves, without FMA The composed maps of reconstruction(), portion(), the LBM stream and the petsc prediction rows take the boundary-clamped stencil (previous commits), so the demos that go through them move. Against the references of hpc-maths#494, built with -ffp-contract=off as they were, nine comparisons fail and are regenerated from the same kind of build: the implicit heterogeneous heat and the lid-driven cavity (petsc rows), the heated cavity (both), and the six adaptive LBM cases (stream and reconstruction). Everything else keeps its bytes.
clang (16 to 19, and Apple's) rejected prediction<r, value_t>(cls, delta_l, idx...) where the trailing indices come from a generic lambda: the first index type explicit, the others deduced, and the class argument's dimension depending on the whole pack. The calls now cast every index to value_t and let the pack be deduced uniformly. gcc accepted the previous form, which is why the local builds did not see it.
…maps move as well With the petsc rows now consistent from hpc-maths#494 on, the lid-driven cavity still moves in this slice: its ink is transported with the reconstruction-based flux at the finest level, which takes the composed maps. Regenerated from the build without FMA.
…al away from every face The same shortcut as the shift query's (see hpc-maths#494): an interval at least `reach` away from every non-periodic face of the box is one interior run, and the decomposition into runs of constant class need not walk it.
gouarin
force-pushed
the
bc-prediction-reconstruction
branch
from
September 3, 2026 17:51
a6bf1c6 to
f9b16f7
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Description
This closes slice 2c of the boundary rewrite: the four consumers of the prediction coefficients
now take the same stencil near a boundary.
compute_detail_opandprediction_opalreadyshifted it inward (#494); the composed maps of
reconstruction.hpp-reconstruction(),portion(),transfer()and the LBM stream built on them - and the prediction rows of thepetsc assembly still used the centred family everywhere, so a field predicted one way was
reconstructed another, and the matrix a solver assembled was not the operator the explicit
path applied.
The difficulty is the memo.
prediction(delta_l, indices)is a mesh-free, global memo of thecomposed map, and it has to stay one: the LBM stream builds one tap list per (level, component)
from it. What the composed map depends on near a boundary is the domain within a fixed reach of
the reference cell, and only that, so the memo is now keyed on a position class - the cover
of every transverse row within
3r + 1of the cell, capped at that reach - from which the shiftrule of #492 can be asked at every step of the composition, for any parent and any gap, without
a mesh.
3r + 1is derived, not chosen: the support of the composed map stays within2r, butthe last step asks the rule about parents within
rand their candidate boxes reach2rfurther, plus one coarse cell for the children a stream slice places next door. The class is
read off the mesh once per run of constant class, by the same box arithmetic or row scan as the
shift, and a mesh with no domain to be positioned against - a uniform mesh - is one interior
class.
The interior is unchanged: the recursion visits the nodes in the order it always did and a
centred parent carries the coefficient 1 exactly, so away from every boundary the maps are bit
for bit those of the centred recursion, which the test asserts against a copy of it kept in the
test.
reconstruction_op_andportiondecompose their coarse row into runs of constant class;the LBM stream builds its tap lists per class and looks one up per run, one per strip in the
bulk; the petsc prediction rows visit the same shifted nodes as the explicit operator, in a
dimension-generic loop that replaced three hand-written copies.
The regenerated references are listed and quantified in the second commit: the two implicit
cases move where the petsc rows changed (1.9e-5 on the heterogeneous heat; the lid-driven
cavity's ink front, sharp and next to the refined lid corners, moves by about one fine cell
while its total is preserved), the adaptive LBM cases by 3e-7 to 3.8e-4, two of them adapting
to a slightly different mesh.
Cost, and the flat table (perf commit)
On
burgers_mra(hat, levels 2 to 12,--timers), the first version of this PR ran 12 %slower than
main(3.21 s against 2.86 s) while every timed section had become cheaper:reconstruction()had gone from 1.42 s to 1.82 s. Not the position query - 16 ns perinterval, 7 ms in total - but one hash lookup per child of every run, into a memo whose key now
carries the class and which is nine times larger on a 1D domain where every coarse level lies
within reach of both boundaries, followed by the iteration of the map's own hash table. The
maps of all the children of one (class, delta_l) are now built once into a contiguous table
indexed by the child's linear index, each map flattened into a sorted vector of terms; the
operator pays one lookup per run and iterates contiguous terms.
burgers_mragoes to 2.79 s,under
main. The sums are unchanged term for term; the sort fixes their order.Where the mesh has no inward reach (fix commit)
#494 found that the clamped stencil is only legitimate on a mesh that holds and fills what it
reads, and made MRMesh declare it. The position-class query follows the shift query: on any
other mesh,
prediction_domain()says so and every run is the interior class, soreconstruction(),portion(),transfer(), the LBM stream and the petsc prediction rows usethe centred maps there, as before.
The references (last commit)
Redone against the references of #494 with demos built with
-ffp-contract=off(see #494 forwhy the sandbox's default build cannot produce a reference the x86 CI matches): nine
comparisons move and are regenerated - the implicit heterogeneous heat and the lid-driven
cavity (petsc rows), the heated cavity (both), the six adaptive LBM cases (stream and
reconstruction). Everything else keeps its bytes. The earlier paragraph on the regenerated
references, and its figures, stands.
Related issue
None. Part of the boundary-machinery rewrite; follows #492 to #494.
How has this been tested?
tests/test_reconstruction_boundary.cppstates the property the way the rest of this rewritedoes: a field holding the cell averages of a polynomial the operator reproduces, every cell
outside the domain poisoned with NaN and the ghost update told it has nothing to do - so a
single read of an outer ghost turns a value into NaN. Then:
reconstruction()is exact everywhere on an adapted mesh whose refinement touches theboundary, in 1D and 2D, and not exact one degree above;
portion()on every child of every coarse real cell, scalar form and slice form, is exact;3 in 2D, children of the neighbouring cells included;
cell for cell.
485/485 tests and the 41 petsc explicit-versus-implicit consistency tests pass with PETSc on
(the latter are the I2 check: the assembled matrix and the explicit operator agree), and the
101 demo comparisons pass with the regenerated references.
Code of Conduct
By submitting this PR, you agree to follow our Code of Conduct