refactor(bc): fill the outer ghosts of a level from that level alone - #501
Open
gouarin wants to merge 33 commits into
Open
refactor(bc): fill the outer ghosts of a level from that level alone#501gouarin wants to merge 33 commits into
gouarin wants to merge 33 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.
5 tasks
Not up to standards ⛔🔴 Issues
|
| Category | Results |
|---|---|
| UnusedCode | 5 medium |
🟢 Metrics 399 complexity · 50 duplication
Metric Results Complexity 399 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-outer-ghosts
branch
13 times, most recently
from
September 3, 2026 13:45
5802bd5 to
80d7adc
Compare
…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.
The outer ghosts used to be filled by an orchestration from the fine levels to the coarse ones: `project_bc` averaged the finer level's outer ghosts into the coarse ones, `predict_bc` copied a coarse outer ghost into the fine ones under it at order 0, `project_corner_below` carried the corner ghosts two levels down with a hardcoded 2, a second polynomial extrapolation filled the far layers next to projection ghosts, and five branches on `is_box()` - a declared flag, false at every level of a plain box mesh but its reference level - decided which of these ran where. All of it existed for one reader, the prediction stencil, which read outside the domain wherever a cell touched the boundary and needed a coarse level's outer ghosts to hold what the fine level's did. That reader is gone: prediction shifts its stencil inward and reads only cells the domain holds, wherever the domain is wide enough for that. The outer ghosts of a level are now the physical extension of what that level holds, written once per level from its cells inside the domain, and nothing crosses levels: - the ghost update fills them in its bottom-up pass, after the prediction ghosts of the level, when every cell of the level inside the domain has a value - the real cells, the projection ghosts and the prediction ghosts - so a coarse level under a refined boundary gets its outer ghosts from the condition applied to its own projected values, at every level down to 0; - the top-down pass no longer touches them: the projection reads cells inside the domain only, and the prediction at the next level reads a level's outer ghosts only where the domain is too narrow for the stencil to shift inward, after they are filled; - the far layers and the corner blocks follow the same rule, and the corner extrapolation stays as the fallback for the schemes that read a diagonal ghost; - a condition is applied around the cells whose whole stencil the mesh holds at that level (`cells_holding_stencil`), per cell, rather than by a flag on the domain: a projection ghost next to an obstacle does not always have, at its own level, the inner neighbours a third-order condition reads - the obstacle demo threw on it. Where a level is `2r` cells wide or less - the coarsest level of a small box - no shift fits and the stencil stays centred, reading the outer ghosts written here at that level. That is the one place the numerical extension is the physical one, and it is where a `static_assert(prediction_stencil_radius <= 1)` used to stand: the assertion guarded `project_bc`'s two-level child search, which is gone with it. Cost on finite-volume-advection-2d: 2.53 s -> 2.31 s total, the outer ghosts from 0.36 s to 0.23 s. A first version applied the conditions in two passes and materialised the source set per direction and ran at 6.5 s: the cost is set traversal and per-evaluation overhead, not the values written. The adaptive von Karman LBM reference is regenerated: the cylinder is a hole in the domain and the coarse levels around it now take the condition applied to their projected values instead of an average of the fine outer ghosts, 6.5e-5 relative at most, same mesh.
…o runs The demo runs main_fct<0>() then main_fct<1>(). The options were registered, on the first call only, on that call's locals - Tf, the file name, the path - and parsed again by the second call: CLI11 then wrote through pointers into a stack frame that no longer existed. AddressSanitizer reports it as a stack-use-after-return in CLI::detail::lexical_cast during the second parse. What it corrupted depended on the second frame's layout: on the x86 CI runner, with this branch's inlining, it hit a std::string of main_fct<1>(), and the run ended with 'free(): invalid pointer' after its last iteration - at one rank, so the whole MPI job stopped there. Unaffected layouts kept the second run's Tf at its default of 0.1 instead of the value on the command line, which is why the pred_1 reference runs longer than pred_0. The parameters now live in one struct, registered and parsed once in main(), and copied into each run. Both runs see the same command line; the pred_1 reference follows.
…pdate The ghost update visits every level from 0 up; a mesh refined everywhere - the uniform max-level mesh burgers_mra carries next to its adapted one - holds one of them. Evaluating the boundary conditions' set expressions on the empty levels cost 0.02 s of the 0.04 s that its time-loop ghost update had gained; they return before building any expression now.
gouarin
force-pushed
the
bc-outer-ghosts
branch
from
September 3, 2026 17:52
747d3dd to
d8e69e9
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 is the phase the rest of the rewrite was for. The outer ghosts - the cells outside the
domain that the mesh holds next to a level's cells - used to be filled by an orchestration
from the fine levels to the coarse ones:
project_bcaveraged the finer level's outer ghostsinto the coarse ones,
predict_bccopied a coarse outer ghost into the fine ones under it atorder 0,
project_corner_belowcarried the corner ghosts two levels down with a hardcoded 2, asecond polynomial extrapolation filled the far layers next to projection ghosts, and five
branches on
is_box()- a declared flag, false at every level of a plain box mesh but itsreference level - decided which of these ran where. All of it existed for one reader, the
prediction stencil, which read outside the domain wherever a cell touched the boundary and
needed a coarse level's outer ghosts to hold what the fine level's did. Since #494 and #500 that
reader shifts its stencil inward and reads only cells the domain holds, wherever the domain is
wide enough for that.
The outer ghosts of a level are now the physical extension of what that level holds,
written once per level from the cells of that level inside the domain, and nothing crosses
levels:
when every cell of the level inside the domain has a value - the real cells, the projection
ghosts (the coarse cells under finer real cells, whose value is the average of theirs) and
the prediction ghosts - so a coarse level under a refined boundary gets its outer ghosts from
the condition applied to its own projected values, at every level down to 0, and a fine
layer under a coarse cell from its own predicted values;
The prediction at the next level reads a level's outer ghosts only where the domain is too
narrow for the stencil to shift inward, and they are filled before it;
stays as the fallback for the schemes that read a diagonal ghost - a lattice-Boltzmann stream
with diagonal velocities whose condition does not own the corner;
(
cells_holding_stencil, an intersection of the translated reference, one per stenciloffset), per cell rather than by a flag on the domain: a projection ghost next to an obstacle
does not always have, at its own level, the inner neighbours a third-order condition reads.
Where a level is
2rcells wide or less - the coarsest level of a small box - no shift fits andthe stencil stays centred, reading the outer ghosts written here at that level. That is the one
place the numerical extension is the physical one, and it is where a
static_assert(prediction_stencil_radius <= 1)used to stand: the assertion guardedproject_bc's two-level child search, which is gone with it.update_outer_ghost.hppgoes from 500 lines to 120, andapply_field_bc.hpplosestranslated_outer_neighbours,apply_extrapolation_bc_ghostsand the second extrapolationpass. The MPI synchronisation count per level is unchanged: the exchange that used to follow
the top-down boundary conditions now follows the bottom-up ones.
Cost
finite-volume-advection-2d, default options,--timers, sequential, three runs each, on thesame machine. The previous slice, #500, is the baseline; the plan measured the outer ghosts at
17.5 % of the runtime on
main.Two design points came out of the measurement, and are worth knowing for what follows:
the top-down pass, to the prediction ghosts in the bottom-up one - and materialised the
source set once per direction. It ran at 6.5 s: the cost of the boundary machinery is set
traversal and fixed per-evaluation overhead, not the values written, so the number of set
expressions evaluated per level is what counts. One pass per level, and no materialisation,
is what the table shows.
of evaluating a set expression. That is the lever for the plan's 5 % target, and it is
independent of this PR.
Against
mainitself (7fa7237), same protocol, three alternating runs, the whole chain up tothis PR:
finite-volume-advection-2d2.34 s against 2.34 s, with the detail computation 15 %cheaper, the mesh update 8 % dearer (the inward reach of #493) and the ghost update 5 % dearer;
burgers_mra(hat, levels 2 to 12) 2.89 s against 2.77 s, the adaptation 6 % cheaper.Regenerated reference
test_lbm_demo_d2q9_von_karman_adaptive-extra1: the adaptive von Karman case, whose cylinder isa hole in the domain. The outer ghosts of the coarse levels around the cylinder used to be
averaged from the fine ones and are now the condition applied to the projected values, so the
solution moves by 6.5e-5 (relative) at most; the mesh is unchanged. Nothing else in the
comparison suite moved, the obstacle demo included. The reference is generated from a build
configured with
-ffp-contract=off, for the reason given in #494: the default gcc/aarch64build fuses multiply-adds, and the x86 CI does not.
A demo bug this PR's stack layout exposed (last commit)
The MPI job failed on
finite-volume-advection-2dat one rank,free(): invalid pointerafterthe last iteration of its second run, on the x86 runner only: not reproducible on aarch64 with
mpich, nor on the runner under ASan or UBSan. The demo runs
main_fct<0>()thenmain_fct<1>(); its options were registered on the first call's locals and parsed again by thesecond, so CLI11 wrote through pointers into a dead stack frame (ASan does report the
stack-use-after-return in
CLI::detail::lexical_cast). What that corrupted depended on thesecond frame's layout, which this branch's inlining changed. The parameters are now parsed once
in
main()and copied into each run; a side effect is that the second run finally sees the--Tfof the command line instead of its default.Related issue
None. Phase 3 of the boundary-machinery rewrite; follows #492 to #500.
How has this been tested?
tests/test_ghost_update.cppandtests/test_bc_ghost_values.cpp: the existing checks of theghost values against the boundary conditions, on box meshes and on adapted meshes whose
refinement touches the boundary, at every order they anchor;
test_prediction_inward_reach.every_cell_the_mesh_holds_is_filled(from fix(mr): give the coarse levels the inward reach a clamped prediction stencil needs #493) is what pinsthe placement of the single pass: applied before the prediction ghosts exist, the outer
ghosts next to them under a refined boundary are held and never written;
tests/test_prediction_domain.cppandtests/test_reconstruction_boundary.cpp(from feat(mr): predict and take details with the boundary-clamped prediction stencil #494 andfeat: reconstruct, stream and assemble with the boundary-clamped prediction stencil #500), which poison every outer ghost with NaN, still pass: the predictions read none of them;
finite-volume-linear-convection-obstacle), whose domain has a hole, iswhat caught the stencil-holding condition, end to end through the demo comparison suite;
pass with the reference regenerated above.
Code of Conduct
By submitting this PR, you agree to follow our Code of Conduct