Large-Scale Startup Optimizations - #1686
Conversation
…nstantiation, saving memroy.
|
Claude Code Review Head SHA: f35421e Files changed:
Findings:
|
…ated memory to pair down collision list
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1686 +/- ##
==========================================
+ Coverage 60.77% 60.83% +0.05%
==========================================
Files 83 83
Lines 20872 20924 +52
Branches 3101 3108 +7
==========================================
+ Hits 12685 12729 +44
+ Misses 6121 6118 -3
- Partials 2066 2077 +11 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
sbryngelson
left a comment
sbryngelson
left a comment
There was a problem hiding this comment.
Overview
Three separable changes bundled together:
- Particle-cloud generation becomes neighborhood-aware (
m_particle_cloud.fpp,m_start_up.fpp) — the real content.get_neighbor_bounds()moves ahead of generation; each packing method now allocates its own working array, filters to the rank's IB neighborhood, and returns a count.gbl_patch_idis assigned at placement time (glbl_idx) instead of being reassigned ins_reduce_ib_patch_array. For lattice packing this makes per-rank memory O(local) instead of O(global). - Capacity constant bumps (
m_constants.fpp) — 15× across the board. - Two drive-bys —
DelayFileAccess→s_delay_file_accessrename, and re-enablings_detect_ib_collisionsovers_detect_ib_collisions_n2.
The core design is sound. The global-ID scheme is deterministic across ranks (glbl_idx advances for every candidate site whether kept or not, and both packers place exactly num_particles), so IDs stay dense and consistent. Moving the moving_ibm scan to particle_cloud(:) specs is required by the filtering and is correct. And moving the num_local_ibs > num_local_ibs_max PROHIBIT inside the loop fixes a genuine out-of-bounds write to local_ib_patch_ids that could previously happen before the check fired.
The problems are mostly in what the constant bumps drag along with them.
High priority
1. The constant bumps have large unconditional memory and per-timestep costs.
ib_patch_parameters is 60 reals + 6 ints + 1 logical ≈ 512 B. With num_ib_patches_max_namelist: 54000 → 810000:
| Structure | Before | After |
|---|---|---|
patch_ib (fixed-size, GPU_DECLARE) |
27.6 MB | 415 MB |
collision_lookup (num_local_ibs_max*27*8, 4) |
6.9 MB | 104 MB |
wall_overlap_distances |
2.6 MB | 39 MB |
recv_bufs in s_handoff_ib_ownership |
26.6 MB | 399 MB |
patch_ib is a fixed-size module array in the generated decls with the GPU flag set, so the 415 MB is resident on device for every simulation run, IB or not — and pre_process carries the same 415 MB in host BSS, which matters at high ranks-per-node on CPU partitions.
Worse, m_ibm.fpp:1429 and :1454 do $:GPU_UPDATE(host='[patch_ib]') / device='[patch_ib]' on the whole array, and s_handoff_ib_ownership runs every time step when moving_immersed_boundary_flag is set (m_time_steppers.fpp:591). That's 415 MB D2H + 415 MB H2D per step, up from 27.6 MB — in a PR whose purpose is performance at scale.
The cheap fix already exists elsewhere in the file: m_ibm.fpp:83/:94 and m_data_output.fpp:923 use patch_ib(1:num_ibs). Scoping lines 1429/1454/1541 the same way costs nothing and removes the regression. Separately, worth considering making patch_ib allocatable and sized from the actual case rather than paying 415 MB statically.
2. The collision-detection swap is an unexplained behavior change.
call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions)
! call s_detect_ib_collisions_n2(num_considered_collisions)git log shows master deliberately runs the _n2 variant and this PR flips it back, with no mention in the description. The motivation is obvious at scale — _n2 is O(num_ibs²) over local patches, hopeless at 810k — but these are not equivalent detectors: the ghost-point path only sees pairs whose IB markers are cell-adjacent, and it dedups serially on the host in O(num_raw × num_collisions), then pushes the full collision_lookup (now 104 MB) to device on every call. That's a different accuracy/cost tradeoff, not a bug fix.
Please state the reasoning in the PR body, and make it a case-file switch or delete the loser — leaving commented-out dead code contradicts #1699, which just landed.
3. The p == 0 → num_dims < 3 fix is only half applied.
Correct and valuable in s_particle_cloud_lattice: in simulation p is the local z-extent, so a rank owning a single z-cell has p == 0 in a 3D run and would silently take the 2D branch — very reachable at 78k ranks. But s_particle_cloud_random_box still has four p == 0 tests (m_particle_cloud.fpp:116, 147, 156, 169) with exactly the same latent bug. Fix both or neither.
4. ib_gbl_idx_lookup is still O(N_global).
m_start_up.fpp:1332 allocates ib_gbl_idx_lookup(1:num_gbl_ibs) on device unconditionally — 2 GB/rank at 500M particles — and s_update_ib_lookup() zeroes it on host and transfers it both directions inside every handoff. Not introduced here, but it's the next wall for exactly the workload this PR targets, and it means the memory win applies to static beds only. Worth a note in the PR so it doesn't get lost.
Medium
5. Argument shadowing. s_reduce_particle_cloud_ibs(cloud_ibs, num_ibs) — the dummy num_ibs shadows the module global num_ibs from m_global_parameters, which the surrounding code reads. Rename to num_cloud_ibs.
6. Fragile single-rank invariant. The num_procs == 1 and no-MPI paths copy num_particle_cloud_ibs entries but then set num_ibs = num_gbl_ibs. That's only correct because f_neighborhood_ranks_own_location short-circuits to .true. for num_procs <= 2 — a cross-module invariant with nothing enforcing it. Use num_namelist_ibs + num_particle_cloud_ibs, or add a PROHIBIT on the mismatch.
7. No case default on packing_method. If it ever falls through, cloud_ibs is unallocated and num_cloud_ibs undefined, so deallocate(cloud_ibs) is a runtime error. m_checker.fpp:194-198 makes this unreachable today, but the failure mode went from benign to a crash — add the default abort.
8. Latent integer overflow, now reachable. hash_size = max(16, 4*particle_cloud(i)%num_particles) overflows default integer above ~537M particles, and max(16, <negative>) degrades silently to a 16-bucket table (O(N²) placement) rather than aborting. Raising num_ib_patches_max to 1.5e9 puts this inside the advertised range.
9. The delay before the collective MPI_FILE_OPEN. MPI_FILE_OPEN(MPI_COMM_WORLD, ...) is collective — rank-staggered arrival can't reduce contention on a single shared-file open; every rank just waits for the largest delay (rank 78000 ⇒ ~6.1M random_number calls). The existing s_delay_file_access call sites all precede per-rank file opens, where staggering does help. If you measured a real win here, please put the numbers in the PR; otherwise I'd drop it.
10. num_ib_patches_max = 1500000000 is dead code. It's referenced only by its own definition — case_validator.py:632 binds a same-named local to num_ib_patches_max_namelist. Delete it rather than bumping it; leaving a 1.5e9 constant named like an array bound is a trap.
11. Undocumented coupling between the two constants. num_local_ibs_max * 27 == num_ib_patches_max_namelist holds before (2000×27 = 54000) and after (30000×27 = 810000), and it has to: wall_overlap_distances(patch_id, i) is indexed by patch_id ∈ [1, num_ibs] with num_ibs bounded by num_ib_patches_max_namelist, while the array's first dimension is num_local_ibs_max*27. Nothing states or enforces this. Add a comment and a compile-time check so a future bump of one constant doesn't silently create an out-of-bounds write.
Minor / style
- Half-finished rename: the subroutine is now
s_delay_file_access, butProcessRank,iDelay,nFileAccessDelayIterations,Dummy,Numberremain CamelCase. glbl_idxvs. the repo's establishedgbl_prefix (gbl_patch_id,ib_gbl_idx_lookup).stat=alloc_stathandling appears only on the random-box allocate; the lattice allocate and the caller'sparticle_cloud_ibsallocate have none.- The PROHIBIT in
s_generate_particle_cloudscompares againstnum_ib_patches_max_namelist, but the array is sizedmin(num_ib_patches_max_namelist, n_total_particles). Safe today (the sum can't exceedn_total_particles), but checkingsize(particle_cloud_ibs)would be self-evidently correct. - Continuation-line alignment at
m_particle_cloud.fpp:296-297.
Test coverage
The tests checkbox is unchecked and nothing in tests/ exercises particle_cloud — the only coverage is examples/2D_mibm_particle_cloud/case.py. More to the point, f_neighborhood_ranks_own_location returns .true. unconditionally for num_procs <= 2, so the entire new filtering path is inert at typical CI rank counts — even adding a 1- or 2-rank test wouldn't touch it.
Worth adding: a ≥3-rank lattice-packing case asserting (a) the union of per-rank particles reconstructs the full set, (b) gbl_patch_id values are unique and dense over 1..num_gbl_ibs, and (c) results match a single-rank run.
Verdict
The generation rework is good and the constants clearly need to move for the target scale. The blockers before merge are the unscoped full-patch_ib GPU transfers in the per-timestep handoff path (#1), a justification and cleanup for the collision-detector swap (#2), and finishing the p == 0 fix in the random-box packer (#3). Everything else can be follow-up, though #10 and #11 are one-liners worth doing now.
sbryngelson
left a comment
sbryngelson
left a comment
There was a problem hiding this comment.
Following up on my earlier comment-only review with a blocking one, and moving the actionable items inline.
Overall: the neighborhood-aware generation rework is good, and the design is sound. The global-ID scheme is deterministic across ranks (glbl_idx advances for every candidate site whether kept or not, and both packers place exactly num_particles), so IDs stay dense and consistent; moving the moving_ibm scan to particle_cloud(:) specs is required by the filtering and correct; and moving the num_local_ibs > num_local_ibs_max PROHIBIT inside the loop fixes a real out-of-bounds write to local_ib_patch_ids.
The blockers are (1) the constant bumps turning the per-timestep patch_ib GPU transfers into 415 MB each way, (2) the unexplained collision-detector swap plus commented-out dead code, and (3) the p == 0 -> num_dims < 3 fix being applied in only one of the two packers. Details inline.
Two items with no diff anchor
ib_gbl_idx_lookup is still O(N_global). m_start_up.fpp:1332 allocates ib_gbl_idx_lookup(1:num_gbl_ibs) on device unconditionally -- 2 GB/rank at 500M particles -- and s_update_ib_lookup() zeroes it on host and transfers it both directions inside every handoff. Not introduced here, but it is the next wall for exactly the workload this PR targets, and it means the memory win applies to static beds only. Worth a note in the PR body so it does not get lost.
Test coverage. The tests checkbox is unchecked and nothing in tests/ exercises particle_cloud -- the only coverage is examples/2D_mibm_particle_cloud/case.py. More to the point, f_neighborhood_ranks_own_location returns .true. unconditionally for num_procs <= 2, so the entire new filtering path is inert at typical CI rank counts -- even adding a 1- or 2-rank test would not touch it. Worth adding a >=3-rank lattice-packing case asserting (a) the union of per-rank particles reconstructs the full set, (b) gbl_patch_id values are unique and dense over 1..num_gbl_ibs, and (c) results match a single-rank run.
Minor / style
glbl_idxvs. the repo's establishedgbl_prefix (gbl_patch_id,ib_gbl_idx_lookup).stat=alloc_stathandling appears only on the random-box allocate; the lattice allocate and the caller'sparticle_cloud_ibsallocate have none.- The PROHIBIT in
s_generate_particle_cloudscompares againstnum_ib_patches_max_namelist, but the array is sizedmin(num_ib_patches_max_namelist, n_total_particles). Safe today (the sum cannot exceedn_total_particles), but checkingsize(particle_cloud_ibs)would be self-evidently correct. - Continuation-line alignment at
m_particle_cloud.fpp:296-297.
| integer, parameter :: num_stl_models_max = 10 | ||
| !> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing) | ||
| integer, parameter :: num_ib_patches_max = 2050000 | ||
| integer, parameter :: num_ib_patches_max = 1500000000 |
There was a problem hiding this comment.
Dead constant -- delete rather than bump.
num_ib_patches_max is referenced only by its own definition. toolchain/mfc/case_validator.py:632 binds a same-named local to num_ib_patches_max_namelist, not to this. Raising it to 1.5e9 has no effect, and leaving a 1.5e9 constant named like an array bound is a trap for whoever reaches for it next.
Related: it also puts hash_size = max(16, 4*particle_cloud(i)%num_particles) in s_particle_cloud_random_box inside the advertised range -- that overflows default integer above ~537M particles, and max(16, <negative>) degrades silently to a 16-bucket table (O(N^2) placement) instead of aborting.
There was a problem hiding this comment.
Constant is gone in 75c7c61, thanks -- and the hash_size overflow goes with it in practice, since s_particle_cloud_random_box allocates placed(3, n_total_particles) and would run out of memory long before 537M particles.
One leftover from the deletion though: the doc comment above it was not removed, so m_constants.fpp:29-30 now reads
!> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing)
!> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction)
integer, parameter :: num_ib_patches_max_namelist = 54000Two stacked !> blocks with the orphaned one binding to num_ib_patches_max_namelist -- which it now describes exactly backwards, since that constant is what sizes patch_ib. Drop the first line.
Reopening just for that one-line fix.
| !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) | ||
| integer, parameter :: num_ib_patches_max_namelist = 54000 | ||
| integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) | ||
| integer, parameter :: num_ib_patches_max_namelist = 810000 |
There was a problem hiding this comment.
Blocking: this 15x bump has large unconditional memory and per-timestep costs.
ib_patch_parameters is 60 reals + 6 ints + 1 logical ~= 512 B. With num_ib_patches_max_namelist: 54000 -> 810000:
| Structure | Before | After |
|---|---|---|
patch_ib (fixed-size, GPU_DECLARE) |
27.6 MB | 415 MB |
collision_lookup (num_local_ibs_max*27*8, 4) |
6.9 MB | 104 MB |
wall_overlap_distances |
2.6 MB | 39 MB |
recv_bufs in s_handoff_ib_ownership |
26.6 MB | 399 MB |
patch_ib is a fixed-size module array in the generated decls with the GPU flag set, so the 415 MB is resident on device for every simulation run, IB or not -- and pre_process carries the same 415 MB in host BSS, which matters at high ranks-per-node on CPU partitions.
Worse: m_ibm.fpp:1429 and :1454 do $:GPU_UPDATE(host='[patch_ib]') / device='[patch_ib]' on the whole array, and s_handoff_ib_ownership runs every time step when moving_immersed_boundary_flag is set (m_time_steppers.fpp:591). That is 415 MB D2H + 415 MB H2D per step, up from 27.6 MB -- in a PR whose purpose is performance at scale.
The cheap fix already exists elsewhere in the same file: m_ibm.fpp:83/:94 and m_data_output.fpp:923 use patch_ib(1:num_ibs). Scoping lines 1429/1454/1541 the same way costs nothing and removes the regression. Separately, worth considering making patch_ib allocatable and sized from the actual case rather than paying 415 MB statically.
There was a problem hiding this comment.
Reverting the bumps removes the memory concern, but it also removes the capability this PR exists for, so I do not think this is settled.
s_reduce_ib_patch_array PROHIBITs on num_local_ibs > num_local_ibs_max (m_start_up.fpp:1296 and :1312). With num_local_ibs_max back at 2000, the branch aborts on its own advertised workload: 500M IBs over 78k ranks is ~6.4k locally-owned IBs per rank, and your note on the collision thread says 81k particles per rank. Either way, 40x over the cap. Did the Frontier run use locally-patched constants? If so, the branch as it stands does not reproduce it, and the constants need to come back.
What I was asking for was the other half of the tradeoff, not the revert: keep the bumps and stop moving the whole fixed-size array.
m_ibm.fpp:1429,:1454,:1541still doGPU_UPDATE(host='[patch_ib]')/GPU_UPDATE(device='[patch_ib]')on the entire array, insides_handoff_ib_ownership, which runs every time step whenmoving_immersed_boundary_flagis set (m_time_steppers.fpp:591). Scoping them topatch_ib(1:num_ibs)-- exactly whatm_ibm.fpp:83/:94andm_data_output.fpp:923already do -- makes the cost proportional to the actual IB count rather than the compile-time cap. That is the change that makes 810000 affordable; without it the cap is priced into every time step whether or not the case has any IBs.recv_bufsin the same routine is sizedbuf_size = storage_size(0)/8 + patch_bytes*num_local_ibs_maxper neighbor over 26 neighbors (m_ibm.fpp:1459). Atnum_local_ibs_max = 30000that is ~399 MB allocated per handoff regardless of how many patches actually migrate. Sizing it from the realnew_count(exchange counts first, then allocate) removes that one too.
With both of those, I have no objection to 810000/30000 -- the point was never that the constants should not move, it was that the cost of moving them was being paid unconditionally by every run.
Reopening because the issue is live in either direction: either the constants come back up and the transfers get scoped, or they stay down and the PR description should say the full-system run is not reachable with the merged code.
| integer, parameter :: num_ib_patches_max_namelist = 54000 | ||
| integer, parameter :: num_local_ibs_max = 2000 !< Maximum number of immersed boundary patches (patch_ib) | ||
| integer, parameter :: num_ib_patches_max_namelist = 810000 | ||
| integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) |
There was a problem hiding this comment.
Undocumented coupling between these two constants.
num_local_ibs_max * 27 == num_ib_patches_max_namelist holds before (2000x27 = 54000) and after (30000x27 = 810000), and it has to: wall_overlap_distances(patch_id, i) is indexed by patch_id in [1, num_ibs] with num_ibs bounded by num_ib_patches_max_namelist, while the array's first dimension is num_local_ibs_max*27.
Nothing states or enforces this. Please add a comment plus a compile-time check so a future bump of one constant does not silently create an out-of-bounds write.
There was a problem hiding this comment.
Still open after 75c7c61 -- the constants went back down, but the relation num_local_ibs_max * 27 == num_ib_patches_max_namelist is exactly as load-bearing at 2000/54000 as it was at 30000/810000, and still unstated. wall_overlap_distances is allocated (num_local_ibs_max*27, 6) at m_collisions.fpp:46 but indexed by patch_id in [1, num_ibs], and num_ibs is bounded by num_ib_patches_max_namelist (m_start_up.fpp:1307). Bump either constant alone and you get an out-of-bounds write with nothing to catch it.
Cheapest fix is to make the relation structural rather than documented:
integer, parameter :: num_local_ibs_max = 2000 !< Max IBs owned by one rank
integer, parameter :: num_ib_patches_max_namelist = num_local_ibs_max*27 !< Fixed capacity of patch_ib: one rank's 27-neighborhoodThen the coupling cannot be broken by editing one line, and the *27 documents where it comes from.
| ! call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) | ||
| call s_detect_ib_collisions_n2(num_considered_collisions) | ||
| call s_detect_ib_collisions(ghost_points, ib_markers, num_gps, num_considered_collisions) | ||
| ! call s_detect_ib_collisions_n2(num_considered_collisions) |
There was a problem hiding this comment.
Blocking: unexplained behavior change, and dead code left behind.
git log shows master deliberately runs the _n2 variant and this PR flips it back, with no mention in the description. The motivation is obvious at scale -- _n2 is O(num_ibs^2) over local patches, hopeless at 810k -- but these are not equivalent detectors: the ghost-point path only sees pairs whose IB markers are cell-adjacent, dedups serially on the host in O(num_raw x num_collisions), and then pushes the full collision_lookup (now 104 MB, per the m_constants.fpp comment) to device on every call. That is a different accuracy/cost tradeoff, not a bug fix.
Please state the reasoning in the PR body, and make this a case-file switch or delete the loser -- leaving commented-out code contradicts #1699, which just landed.
There was a problem hiding this comment.
It is a somewhat straightforward exercise to show that the cell-adjacency is a necessary requirement of spherical/circular collisions. There is no accuracy tradeoff. Previously when I was running at hundreds of local IBs, the cost of the two subroutines were similar, but one of more clear how it worked (_n2). Now at scale, with 81k particles per rank, it is clear that the performance is not similar in this limit, and the tradeoff is justified. I will follow deletion of ht eunused code, as I was not sure which I would use for a while. But this is not an accuracy loss.
There was a problem hiding this comment.
Thanks -- dead code is gone, and I take the point that at 81k particles per rank the O(N^2) variant is not a real option. Two things still worth doing.
Put the reasoning in the PR description. Right now the body does not mention the swap at all, so the next person reading git log sees master deliberately choose _n2, then this PR silently choose the other one. One sentence ("_n2 is O(num_ibs^2) per rank and untenable above ~1e4 local particles") closes that.
On "no accuracy tradeoff": the continuum statement is fine -- overlapping spheres necessarily have surfaces in contact. But the implemented test is not the continuum one. It is cell-centered ib_markers adjacency over a +/-1 stencil (m_collisions.fpp:266-276), which is a discrete proxy, and it differs from the geometric test in three ways that are worth knowing about even if none of them is a reason to keep _n2:
- The test is one-sided.
if (gp_patch_id < neighbor_patch_id)records a pair only from the lower-encoded-id side, so the pair is seen only if the lower-id particle has a ghost point adjacent to the higher-id marker. Ghost points require a fluid cell withingp_layers(m_ibm.fpp:632-644); in a dense bed a contact-region cell can be enclosed by solid markers of both particles and fail that. You already canonicalize the ordering in the dedup loop (the swap at:306-313), so changing<to/=removes the asymmetry for free -- it doublesnum_raw, which is still bounded bynum_gps. exit neighbor_searchat:281stops after the first differing marker, so a ghost point touching two other particles contributes one pair. Fine when contacts span many ghost points; marginal for a one-cell contact patch in a tightly packed bed.- Detection is grid-resolution-limited. At zero overlap the nearest A-marked and B-marked cell centers can be ~2dx apart, outside a +/-1 stencil, so a pair is registered only once overlap reaches O(dx) rather than at
r1 + r2. For particles resolved with a few cells per radius that is a real slice of the soft-sphere force curve. That is a resolution requirement, not a bug -- but it is a resolution requirement that_n2did not have, so it should be stated somewhere rather than described as equivalent.
Missing bound check. num_considered_collisions increments at :325 with no check against size(collision_lookup, 1) = num_local_ibs_max*27*8. That was tolerable when this was one of two paths; now that it is the only detector, an overflow here is a silent out-of-bounds write into a GPU-resident array. Please add a PROHIBIT.
Reopening for the PR-description note, the < -> /= change, and the bound check.
| @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & | ||
| & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") | ||
| do i = 1, num_bed_ibs | ||
| do i = 1, num_particle_cloud_ibs |
There was a problem hiding this comment.
Fragile cross-module invariant.
This copies num_particle_cloud_ibs entries but then sets num_ibs = num_gbl_ibs two lines down. That is only correct because f_neighborhood_ranks_own_location short-circuits to .true. for num_procs <= 2, so the generation-time filter is a no-op here -- an invariant living in m_collisions.fpp with nothing enforcing it from this side. If that short-circuit ever changes, this silently publishes uninitialized patch_ib entries.
Use num_ibs = num_namelist_ibs + num_particle_cloud_ibs, or add a PROHIBIT on the mismatch. Same applies to the no-MPI branch at line 1322.
There was a problem hiding this comment.
Still open. To be concrete about the failure mode, since it is not obvious from the diff:
The num_procs == 1 branch (m_start_up.fpp:1269-1283) writes num_particle_cloud_ibs entries into patch_ib, then sets num_ibs = num_gbl_ibs, where num_gbl_ibs = num_namelist_ibs + sum(particle_cloud(:)%num_particles). Those two counts agree only if generation kept every particle it placed -- which is true today solely because f_neighborhood_ranks_own_location short-circuits to .true. for num_procs <= 2. That is an invariant owned by m_collisions.fpp, asserted nowhere, and consumed here.
If that short-circuit ever gains a real check (say someone makes the 2-rank case exact), this branch publishes patch_ib(num_particle_cloud_ibs+1 : num_gbl_ibs) uninitialized, local_ib_patch_ids indexes into it, and ib_gbl_idx_lookup is allocated to a count that no longer matches -- all without a crash. Same code at the no-MPI branch, :1320-1331.
Either of these closes it:
num_ibs = num_namelist_ibs + num_particle_cloud_ibsor keep the current line and add
@:PROHIBIT(num_namelist_ibs + num_particle_cloud_ibs /= num_gbl_ibs, &
& "Particle cloud generation dropped patches on a single rank.")The first is better -- it makes the count come from what was actually written.
| file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc) | ||
|
|
||
| call s_mpi_barrier() | ||
| call s_delay_file_access(proc_rank) |
There was a problem hiding this comment.
Does this actually help?
MPI_FILE_OPEN(MPI_COMM_WORLD, ...) is collective -- staggering rank arrival cannot reduce contention on a single shared-file open; every rank just ends up waiting for the largest delay (rank 78000 => ~6.1M random_number calls). Every other s_delay_file_access call site precedes a per-rank file open, where the staggering does do something.
If you measured a real win here at full-system scale, please put the numbers in the PR -- otherwise I would drop the barrier and the delay.
There was a problem hiding this comment.
Still open, and I want to be precise about why this one differs from the other call sites, because the rename made them all look alike.
Every pre-existing s_delay_file_access call precedes a per-rank open:
m_data_output.fpp:716->MPI_FILE_OPEN(MPI_COMM_SELF, ...)on<t_step>_<rank>.datm_data_output.fpp:964-> per-rankib_state_<t_step>_<rank>.datm_boundary_io.fpp:133/:236-> per-rankbc_<rank>.datpre_process/m_data_output.fpp:457-> same pattern
There, staggering spreads 78k independent metadata operations across time, which is the thing that helps.
The new site is before MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ...) on a single shared ib.dat (:905). That call is collective: no rank returns from it until all ranks have entered, so staggering arrival cannot spread anything -- it just idles the early ranks. The added wall time per write is the maximum delay over all ranks, (num_procs/128)*10000 random_number calls, which at 78k ranks is ~6.1M iterations on the last rank, paid on every ib.dat write. The s_mpi_barrier() on the preceding line makes it strictly worse: it forces every rank to wait for the slowest before any of them starts its delay.
If you measured a win at full scale I will take the numbers over my reasoning -- MPI-IO on Lustre does surprising things, and it is possible you were working around something in the file-creation path specifically. But as written I would expect this to cost ~seconds per output step at 78k ranks and buy nothing, so please either post the timings or drop both lines.
|
|
||
| !> Introduce a rank-dependent busy-wait delay to stagger parallel file access and reduce I/O contention. | ||
| impure subroutine DelayFileAccess(ProcessRank) | ||
| impure subroutine s_delay_file_access(ProcessRank) |
There was a problem hiding this comment.
Rename is the right direction, but it stopped at the subroutine name -- ProcessRank, iDelay, nFileAccessDelayIterations, Dummy and Number are all still CamelCase. Worth finishing while you are in here.
There was a problem hiding this comment.
Still open after 75c7c61 -- ProcessRank, iDelay, nFileAccessDelayIterations, Number and Dummy are all still CamelCase in m_delay_file_access.f90:21-33. Five identifiers in one 12-line subroutine; worth finishing since the rename already touched every call site.
| impure subroutine s_generate_particle_clouds(particle_cloud_ibs) | ||
| !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Each packing | ||
| !! method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / s_particle_cloud_random_box) and | ||
| !! hands back only the entries that fall within this rank's IB neighborhood |
There was a problem hiding this comment.
75c7c61 trimmed this docstring down, and it dropped the one sentence that was actually load-bearing:
Only the first num_particle_cloud_ibs of them are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is left uninitialized.
That is the invariant that makes size(particle_cloud_ibs) a bug rather than a synonym, and it belongs on the producing side where someone changing the allocation will see it -- the array is allocated to min(num_ib_patches_max_namelist, n_total_particles) at :48, so capacity and valid-entry count genuinely differ. s_reduce_ib_patch_array still states it on the consuming side (m_start_up.fpp:1233-1235), but that only helps the caller who already knows to look.
The gbl_patch_id sentence is worth keeping too -- "already the final, absolute global patch id, so s_reduce_ib_patch_array copies it as-is" is not deducible from the signature, and it is the thing a reader has to trust to believe the IDs stay dense across ranks.
Happy to see the rest go, just not those two.
|
Went through 75c7c61 and your two replies. Status roundup so nothing gets lost. Closed:
Reopened, because the commit did not address them:
Still open, no response yet: the Blockers for merge, in order: the constants question (it determines whether half of this list is live), the collision bound check, and the |
2 participants
Description
I performed a weak scaling test to the full of OLCF Frontier. To get the simulation to run correctly, the initialization process in simulation had to be modified to accommodate the number of IBs in the simulation as well as specific slowdowns I encountered at the full system size. This branch contains those modifications.
Some initialization optimizations are exclusive to the lattice packing method.
Type of change (delete unused ones)
Testing
Ran 500 million immersed boundaries on 78k MPI ranks.
Checklist
Check these like this
[x]to indicate which of the below applies.See the developer guide for full coding standards.
GPU changes (expand if you modified
src/simulation/)AI code reviews
Reviews are not retriggered automatically. To request a review, comment on the PR:
@claude full review— Claude full review (also triggers on PR open/reopen/ready)claude-full-review— Claude full review via label