Repository navigation
SC26 SCC: Parallelize the cell loops inside each IB patch's gang - #2097
aliencaocao wants to merge 3 commits into
Conversation
|
Claude Code Review Head SHA: 0290a49 Files changed:
Findings:
|
Lines of Code
|
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #2097 +/- ##
=======================================
Coverage 61.95% 61.95%
=======================================
Files 86 86
Lines 22909 22909
Branches 3379 3379
=======================================
Hits 14193 14193
Misses 6207 6207
Partials 2509 2509 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@danieljvickers can you comment on this before I/Spencer hit merge? |
|
our benchmarks in CI @aliencaocao can't replicate any of your speedups. maybe it's a failure of the benchmarking toolchain. not sure. but just one delta on one compiler on one case isn't convincing. moving this to an unapproved draft PR. |
|
@aliencaocao You ran this on a test with 20 IB particles. That parallelism was designed for tens of thousands of particles on the same device. Did you run any analysis on how this compares in the actual larger scale case. It was my experience that the speed up vs what you are sugesting is about 10x-50x using a vector parallelism for everything, which is the current default. My estimation is that this will be slower once your particle count saturates the number of streaming multi-processors on you device. @sbryngelson I don't think any of the benchmarks CI cover the many-ib-parallelism case. This goes back to my previous comment about wanting to add a moving IBM case to the benchmark suite at some point soon. |
|
@danieljvickers at that rate we may as well add 15 other tests to the benchmarks. perhaps we should think about how benchmarking can be better automated. |
|
Getting off topic, but yes that is why I haven't done that yet. Relevant to the conversation, you should not expect this optimization to appear in the CI. |
|
Hi @danieljvickers following your suggestion I ran some extra tests. Environment: 1x H200 SXM, NVHPC 26.9 OpenACC
So the gain shows up when there are few patches with large bounding boxes (20 particles), which is the case the PR targets. At 2,000 and 20,000 particles all three variants are within ~1%, and I see no evidence that it is "10-50x" slower as you mentioned, at 20k scale. nsys on the 20,000-particle case:
At 20k particle size the marking is under 0.4% of GPU time whichever path is used, but even in this case, the PR saves ~30% for this specific kernel. The difference is simply not visible at 20k scale due to it being a smaller part of the total runtime. However data suggests that there is no harm to it, and it is strictly a positive for cases with less particles. If we want CI to be able to see this path, add a benchmark that sets case.py (self-contained; derived from
|
|
@aliencaocao I will refer you to this PR when I merged this feature #1549. In it, I presented this plot. This was 8k particles on a single device with 10 grid cell diameters being placed in 3D. The performance difference was x150 when using many-ib parallelism. Your tables do not have units nor are they clear, so you are going to have to correct me if I interpret these incorrectly. Your results show that the grid cell parallelism is actually faster than many-ib, even with 20k particles on a device. Two separate subroutines with completely separate parallelism strategies that deviate this heavily should not produce identical results to 3 digits in all of these cases, and it makes your results appear questionable. To check, I ran the case file your provided on an A100 for the 2k particle case and measured the performance of the projection subroutine on master for each parallelism. I get this table
THis table shows that for your case, master with grid cell parallelism should be 78x slower. Since these numbers deviate by nearly two orders of magnitude from what you present, I don't really believe anything you are presenting here. |
|
@danieljvickers thanks, you were right on your accounts and I apologise for the confusion. In my previous comment the "master with However, my mistake did not change the conclusion of the PR, as we are comparing actually between this PR and I reran everything with the flag verified in the generated I ran on 2 cases, both 20 steps: the 2D case from my previous comment (5000x5000, 20 / 2,000 / 20,000 particles), and a 3D case in the regime of #1549: 8,000 spheres of ~300 grid cells each (radius 4.15 cells) in a periodic 320^3 box, same moving-IB/collision settings (case.py below). grind, ns/gp/eq/rhs (lower is better):
I believe the 3D case is a measurement noise. To further prove this, I ran nsys which give a much clearer breakdown: Marking cost per call of
Please try my PR on A100 and let me know if you cannot reproduce the PR's speedup. 3D case.py (self-contained)Run as #!/usr/bin/env python3
"""3D uniform flow past a bed of free-floating spheres (moving IBM), to reproduce the regime of
MFC PR #1549: thousands of small spheres (~300 grid cells each), timing s_apply_ib_patches.
Env (defaults): GRID cells per side (320), STEPS (20), NP particles (8000),
RC sphere radius in grid cells (4.15 -> ~300 cells/sphere), MANY many_ib_patch_parallelism (T).
"""
import json
import math
import os
GRID = int(os.environ.get("GRID", 320))
STEPS = int(os.environ.get("STEPS", 20))
NP = int(os.environ.get("NP", 8000))
RC = float(os.environ.get("RC", 4.15))
L = 1.0 # domain [-L/2, L/2]^3
dx = L / GRID
radius = RC * dx
bed = 0.8 * L # cubic bed, centred
dt = 5.0e-4 * (320.0 / GRID)
case = {
"run_time_info": "F",
"x_domain%beg": -0.5 * L,
"x_domain%end": 0.5 * L,
"y_domain%beg": -0.5 * L,
"y_domain%end": 0.5 * L,
"z_domain%beg": -0.5 * L,
"z_domain%end": 0.5 * L,
"cyl_coord": "F",
"m": GRID,
"n": GRID,
"p": GRID,
"dt": dt,
"t_step_start": 0,
"t_step_stop": STEPS,
"t_step_save": STEPS,
"num_patches": 1,
"model_eqns": "5eq",
"alt_soundspeed": "F",
"num_fluids": 1,
"mpp_lim": "F",
"mixture_err": "T",
"time_stepper": "rk3",
"weno_order": 5,
"weno_eps": 1e-16,
"weno_avg": "T",
"avg_state": "arithmetic",
"mapped_weno": "T",
"null_weights": "F",
"mp_weno": "T",
"riemann_solver": "hllc",
"wave_speeds": "direct",
"fd_order": 2,
"bc_x%beg": -1,
"bc_x%end": -1,
"bc_y%beg": -1,
"bc_y%end": -1,
"bc_z%beg": -1,
"bc_z%end": -1,
"ib": "T",
"num_ibs": 0,
"many_ib_patch_parallelism": os.environ.get("MANY", "T"),
"collision_model": 1,
"coefficient_of_restitution": 0.9,
"collision_time": 20.0 * dt,
"ib_coefficient_of_friction": 0.1,
"num_particle_clouds": 1,
"particle_cloud(1)%x_centroid": 0.0,
"particle_cloud(1)%y_centroid": 0.0,
"particle_cloud(1)%z_centroid": 0.0,
"particle_cloud(1)%length_x": bed,
"particle_cloud(1)%length_y": bed,
"particle_cloud(1)%length_z": bed,
"particle_cloud(1)%num_particles": NP,
"particle_cloud(1)%radius": radius,
"particle_cloud(1)%mass": 0.25,
"particle_cloud(1)%min_spacing": 0.5 * radius,
"particle_cloud(1)%moving_ibm": 2,
"particle_cloud(1)%seed": 42,
"particle_cloud(1)%packing_method": 1,
"format": "silo",
"precision": "double",
"prim_vars_wrt": "T",
"E_wrt": "T",
"ib_state_wrt": "F",
"parallel_io": "F",
"patch_icpp(1)%geometry": 9,
"patch_icpp(1)%x_centroid": 0.0,
"patch_icpp(1)%y_centroid": 0.0,
"patch_icpp(1)%z_centroid": 0.0,
"patch_icpp(1)%length_x": L,
"patch_icpp(1)%length_y": L,
"patch_icpp(1)%length_z": L,
"patch_icpp(1)%vel(1)": 0.5,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%vel(3)": 0.0,
"patch_icpp(1)%pres": 1.0,
"patch_icpp(1)%alpha_rho(1)": 1.4,
"patch_icpp(1)%alpha(1)": 1.0,
"fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
"fluid_pp(1)%pi_inf": 0.0,
}
print(json.dumps(case))The 2D case is the same as from my previous comment. |
|
For your 20k particle case, I am not getting the same answers as you on an A100, I get that these methods are identical, which makes sense because it utilizes the same amount of the GPU. I think your measurement is noise at this level. But in the smaller cases I am seeing a speedup of the kernel that about aligns with the increased theoretical usage of the device. I will run this on AMDFlang and make sure that it gives results that are not harmful to performance. If so I will return here and approve. |
Parallelize the cell loops inside each IB patch's gang
Problem
Currently,
s_apply_ib_patches_ib_parallelismhadpatch_idas its only parallel loop. For cases with few immersed bodies whose bounding boxes cover many cells, they leave the GPU nearly idle, because the code uses one thread per patch and every cell of every box is being processed serially by that thread.What this does
The bounding-box cell loops (k, j, i in 3D; j, i in 2D) become a collapsed vector loop inside each patch's gang.
grid_cell_parallelism(the default for static patches) already parallelizes the cell loops and is untouched.Measurements
Benchmark: https://github.com/MFlowCode/MFC/blob/master/examples/2D_mibm_particle_cloud/case.py
(
many_ib_patch_parallelism=T), scaled to 7000x7000 on 1xH200.Validation
On NVHPC 26.9 with OpenACC, taking the saved files on the 20th step and comparing against
master:ib_state_20.dat): 74/400 words differ by one ULP due to floating-point ordering, rest bit-identical.On AMD OpenMP and CPU this code change is a no-op. On NVHPC OpenMP, it just add a
!$omp loop bind(thread)which is also a no-op compared to no directive (current).Contribution Policy
All contributions are expected to demonstrate:
If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.
Acknowledgment
This PR is cherry-picked from #2093 and is contributed by Nanyang Technological University Team Supernova for SC26 Student Cluster Competition.