Fix slow CPU halo update for the middle (y) decomposed dimension - #124
Merged
Merged
Conversation
memcopy! flattened both operands with view(.,:) before copying. For the middle-dimension face slice view(A,:,iy,:) (IndexCartesian) this forces a per-element linear->Cartesian index conversion in the copy, making the y-face halo 25-100x slower than the x/z faces on the CPU (the cost grows with the face size). The x- and z-face slices are IndexLinear, so they were unaffected. Branch on IndexStyle(dst, src): keep the flat-vector memcopy for the linear case (x/z faces, and all the flat send/recv buffers), and copy the shaped views directly when either operand is IndexCartesian (the y face) so iteration stays contiguous along the inner dimension. Output is bit-identical. Measured on update_halo! (np=2, Float32, system OpenMPI), y-face per call: N=66: 882.7 -> 25.9 us (34x) N=130: 2836 -> 41.3 us (69x) N=258: 13794 -> 119.7 us (115x) x and z unchanged. CPU path only; the GPU (CUDAExt/AMDGPUExt) packing kernels and CUDA-aware MPI buffers do not use this function. The Polyester path already iterates with eachindex (Cartesian-aware) and is unaffected. Assisted-by: Claude:claude-opus-4-8
Collaborator
|
Thanks for catching that one! LGTM |
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #124 +/- ##
=========================================
+ Coverage 2.91% 3.56% +0.64%
=========================================
Files 23 25 +2
Lines 754 1038 +284
=========================================
+ Hits 22 37 +15
- Misses 732 1001 +269 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
omlins
approved these changes
Jul 2, 2026
Collaborator
|
Thanks a lot for the PR @b-fg ! |
b-fg
added a commit
to WaterLily-jl/WaterLily.jl
that referenced
this pull request
Jul 2, 2026
Brings the grid-independent Poisson stopping criterion to mpi-igg. tol (default 2e-3) is the max-norm knob max|r| < tol; the bulk term is the mean residual Σ|r|/N < tol/10 (L1/N, in max-norm units). Convergence requires BOTH. Supersedes the earlier RMS-only (L2/N) form and its superseded MPI port c4cec26. MPI adaptations (dispatch-free, via p.inslen): - src/parallel.jl: add local_sumabs(a)=sum(abs,a), the L1 reduction (didn't exist; only local_dot/local_sum were present). - src/Poisson.jl: L₁(p)=global_allreduce(local_sumabs(p.r)) (ghosts≡0, so the whole-array sum = global interior, no double-count); l1n_tol(p::Poisson,tol)=(tol/10)*p.inslen uses the GLOBAL interior count (not length(inside), which is rank-local under MPI); L∞ keeps mpi-igg's global_max form; both solver!s take PR 307's combined break and keep mpi-igg's pin_pressure!+comm! ending. - test/test_poisson.jl: PR 307's bounds + the new L∞ < 2e-3 checks. FP64 accumulation — investigated, NOT needed. Julia sum(abs,·) is pairwise (err ~log₂N·eps): rel.err ~1e-7 at 1e9 cells for both uniform and peaked residuals, four orders below the ~1% the criterion tolerates. The earlier Σr² saturation (~290³) was the naive BLAS-dot accumulator (p.r⋅p.r), not the summation — and L₂ is no longer in the stopping path. GPU reductions are tree-based; the MPI Allreduce sums only nprocs partials. Kept plain FP32. (bench/tol_percell/fp64_l1_investigate.jl) Also bump ImplicitGlobalGrid compat to 0.17.1 (registered 2026-07-02, carries the y-halo fix, upstream PR eth-cscs/ImplicitGlobalGrid.jl#124). Validated: serial grid-independence mean_iters=1.0 at 32/64/128 (L∞ 3.6e-6→8e-8); MPI np=1 ≡ serial bit-identical; np=4 parity max|u| −0.13%, ke −0.17%; criterion lowers np=4 iters 97→80 vs the old tol. KA suite 293/293 (incl. the 2 new L∞<2e-3 checks); MPI suite 23/23. Assisted-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
21 of 25 tasks
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
4 participants
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.
On CPU arrays,
update_halo!is ~25–115× slower (growing with face size) when the middle (y) dimension is the decomposed one than for x or z. This PR closes #123.The issue is that
memcopy!flattens both operands withview(.,:)before copying. The y-face sliceview(A,:,iy,:)isIndexCartesian, while x/z faces areIndexLinear. So the flat copy forces a per-element linear→Cartesian index conversion. And this is not a stride effect since the x-face is the most scattered yet fast.The fix is to branch
memcopy!onIndexStyle(dst, src): keep the flatview(.,:)memcopy for theIndexLinearcase (x/z faces + all flat buffers), elsecopyto!(dst, src)on the shaped views so iteration stays contiguous along the inner dimension. Using the combinedIndexStyle(dst,src)matters:read_h2h!passes the y-slice as the destination.I have conducted the tests only on CPU backend, and I think GPU backends use a different function for this. Also I have verified this in WaterLily.jl and in the linked issue MRWE. below is the WaterLily weak scalability plot using master and this PR. In the MRWE, full
update_halo!(np=2, Float32) results in y-face 882→26 µs (66³), 13794→120 µs (258³) with x/z unchanged.test/test_update_halo.jlpasses 922/922.PR and debugging assisted by Claude Opus 4.8.