[sbryngelson]

IBM: parallelize the IB marker kernel over patch and bounding-box cell

One commit on origin/master 175ec9f (after #2207). It changes src/simulation/m_ib_patches.fpp (the marker kernel behind many_ib_patch_parallelism), one line in src/simulation/m_ibm.fpp, adds a test to toolchain/mfc/test/cases.py and the goldens of its four variants.

Summary

many_ib_patch_parallelism = T draws the immersed-boundary markers with s_apply_ib_patches_ib_parallelism. That kernel ran one GPU thread per patch, each walking its bounding box serially. A bed of a few hundred small bodies therefore used a few hundred threads, and one large body (a wall, a plate) held the whole kernel for as long as it took one thread to walk its box. With moving bodies the markers are rebuilt every Runge-Kutta stage, and in 3D with periodic faces the kernel is launched once per periodic image (27 launches), so the cost is paid 3 x 27 times a step.

Now:

• s_get_ib_boxes (new, one kernel) stores the bounding cells of every patch in every periodic image in ib_box (6 integers per patch and image, allocated on first use, freed in s_finalize_ibm_module) and returns the largest box in cells with a max reduction.
• The marker kernels collapse over (patch, box cell) instead of looping over patches only. A thread reads its box from ib_box, leaves at once if its cell is outside the patch's own box, and otherwise does what the old loop body did.
• The geometry tests and the atomic max on the marker are untouched. The atomic max already made the marker independent of thread order where bodies overlap (higher encoded id wins), so the markers do not change.
• Test: IBM -> Periodic Bodies -> many_ib_patch_parallelism (2D: circle, rectangle, ellipse; 3D: sphere, cuboid, cylinder; one body straddling each periodic face and moving across it), 1 and 2 fluids, so four goldens. No existing test used the flag in a periodic domain; the only user was the 2D_mibm_particle_cloud example.

AI disclosure: made with Claude Code.

The markers are unchanged

To compare master and this branch, a temporary wrapper around s_apply_ib_patches (not part of the PR) hashed the whole marker array, halo cells included, after every call, and the runs also wrote the restart data. Both builds use many_ib_patch_parallelism = T, prescribed motion (moving_ibm = 1, translation plus spin) so that markers are rebuilt at every stage, random positions, and a tenth of the bodies started on a periodic face.

case	calls compared	bodies
3D periodic, 104^3 cells, 120 steps	361	300 spheres
3D walls (no periodicity), 100 steps	301	300 spheres
3D static	1	100 spheres
3D periodic, 2 ranks	181 per rank	300 spheres
3D mixed sizes	25	30 spheres, sphere r = 3 D, 2 cuboids, cylinder
2D periodic, 673^2 cells, 250 steps	751	400 circles
2D periodic, 2 ranks	361 per rank	400 circles
2D mixed sizes, 60 steps	181	40 circles, rectangle, ellipse, airfoil
2D/3D, 10 to 2000 bodies, 15 steps (the speed cases below)	46 each	up to 2000
2D/3D with one plate spanning 90 % of the box, 500 or 2000 bodies	7 and 10	plate + spheres/circles

Result on both machines (below): the hash and the number of marked cells agree at every call; every restart file (conservative variables and the marker file ib.dat) is byte-identical (cmp). The only restart files that differ are the body-state files (ib_state_*.dat, the forces and torques), by at most 6e-17 relative: the forces are summed with atomics, and the same binary differs from itself by the same amount (master against a second run of master: 6e-17).

machine	compiler	GPU build	cases compared	result
Tuolumne, MI300A	LLVM flang 24.0.0 (AFAR 24.3.0)	--gpu mp --mpi	20	identical
A100 80GB node	nvfortran 25.11	--gpu acc --mpi	20	identical
NVIDIA workstation CPU	nvfortran 25.11	--no-gpu --mpi	4 (2D and 3D, 1 and 2 ranks)	identical

All IBM-labelled tests (the 105 listed by --list that contain "IBM", 93 of them non-chemistry) were also run with the flag forced on in both builds (a debug switch, not in the PR) and the goldens they generate compared. Master against this branch: the same 14 tests differ in a handful of entries each (1 to 73 entries out of 6 760 to 436 000), by at most 1e-14 absolute (4e-15 relative to the largest entry); these are two-way and force-coupled tests, and master run twice differs in the same 14 tests by the same amount. On the A100 node, 10 of them.

Speed

Time of one call of s_apply_ib_patches (system_clock around the call, synchronized by a one-element device-to-host update), one rank, prescribed motion with bodies moving every step, periodic box, random positions at about 3 % (3D) and 10 % (2D) volume fraction, 6 (3D) and 12 (2D) cells across a body. Median over calls 4 to 46 of a 15-step run (three calls a step); the mixed-size rows are 5 (3D) and 8 (2D) steps, the plate rows 2 (3D) and 3 (2D) steps, so they are medians of 13 to 22 and 4 to 7 calls. "branch" is this PR, "master" is 175ec9f.

NVIDIA A100, nvfortran 25.11, --gpu acc (two replicates agree to within 6 %; the first is shown):

case	grid	master ms	branch ms	speedup
3D, 10 spheres	34^3	6.52	1.80	3.6
3D, 100 spheres	72^3	7.82	1.82	4.3
3D, 500 spheres	124^3	8.77	2.07	4.2
3D, 2000 spheres	196^3	7.97	2.63	3.0
2D, 10 circles	106^2	1.21	0.61	2.0
2D, 100 circles	336^2	1.29	0.63	2.0
2D, 500 circles	752^2	1.24	0.65	1.9
2D, 2000 circles	1504^2	1.47	0.75	2.0
3D, 500 spheres + sphere r = 3 D, 2 cuboids, cylinder	88^3	176.5	30.3	5.8
2D, 500 circles + rectangle, ellipse, circle, airfoil	752^2	246.4	2.82	87
3D, 500 spheres + plate 0.9 L x 2 D x 0.9 L	124^3	5242	812	6.5
2D, 2000 circles + plate 0.9 L x 2 D	1504^2	1162	499	2.3

AMD MI300A (Tuolumne), LLVM flang 24.0.0 (AFAR 24.3.0), --gpu mp; the last column is the speedup of the median of the last 15 calls (see the second note):

case	grid	master ms	branch ms	speedup	speedup, last 15 calls
3D, 10 spheres	34^3	74.8	4.8	15	5.2
3D, 100 spheres	72^3	75.4	5.0	15	5.4
3D, 500 spheres	124^3	87.1	6.9	13	5.5
3D, 2000 spheres	196^3	87.7	9.2	9.6	5.3
2D, 10 circles	106^2	9.75	1.34	7.3	7.3
2D, 100 circles	336^2	8.61	1.42	6.0	6.0
2D, 500 circles	752^2	7.37	1.55	4.8	4.8
2D, 2000 circles	1504^2	10.6	2.88	3.7	3.7
3D, 500 spheres + sphere r = 3 D, 2 cuboids, cylinder	88^3	2471	12.8	190	190
2D, 500 circles + rectangle, ellipse, circle, airfoil	752^2	689	2.97	230	230
3D, 500 spheres + plate	124^3	121 162	260	470	460
2D, 2000 circles + plate	1504^2	27 246	145	190	190

Notes on these numbers:

• Small and uniform beds: the old kernel was bound by the serial walk of one box and the 27 (3D) or 9 (2D) launches per call, the new one by the launches. On the A100 that is a 2 to 4 times faster call; on the MI300A 4 to 15 times.
• On the MI300A every step, in master and in this branch alike, takes about 1.5 times longer from the tenth step on (a 3D case with 10 spheres: 0.205 s then 0.31 s per step in master, 0.054 s then 0.095 s in the branch; run without the timing wrapper). I did not find the reason and it is not caused by this change. In the branch's marker call of the 3D equal-body cases it shows as about 4.8 ms for the first 27 calls and 15 to 18 ms from then on, in master as little (the call is long, 75 to 96 ms), so the speedup in a longer run is 5.2 to 5.5 times (last column) rather than 9.6 to 15. The 2D and mixed rows do not show it. The A100 shows nothing like it.
• Mixed sizes and plates: the old kernel is as slow as its slowest body, and a plate takes seconds to minutes; the new one is limited by (number of bodies) x (largest box) threads that mostly leave at once. That is still much faster than master in every case tried, but it is wasted work that grows with both numbers.
• Whole-step times, 1 rank, same cases with the fluid solve on (median step):

case	machine	master s/step	branch s/step
3D, 300 spheres, periodic, 104^3	MI300A	0.339	0.126
3D, 300 spheres, periodic, 104^3	A100	0.347	0.328
3D, 300 spheres, walls, 104^3	MI300A	0.072	0.067
3D, 300 spheres, walls, 104^3	A100	0.322	0.324
2D, 400 circles, periodic, 673^2	MI300A	0.081	0.057
2D, 400 circles, periodic, 673^2	A100	0.128	0.126
2D, 40 circles + 3 large bodies, 240^2	MI300A	0.975	0.058
2D, 40 circles + 3 large bodies, 240^2	A100	0.260	0.018
3D, 30 spheres + 3 large bodies, 88^3	MI300A	0.218	0.108
3D, 30 spheres + 3 large bodies, 88^3	A100	0.649	0.243

So the step time moves where the marker update was a visible part of the step: on the A100 a bed of equal bodies gains 0 to 6 % of the step, a mixed-size case 2.7 to 14 times. On the MI300A the old kernel was a large part of the step and the gain is up to 2.7 times for equal bodies and 2 to 17 times for mixed sizes.

• Non-GPU builds: the collapsed loops run serially, and every (patch, cell) pair is visited, including the empty periodic images, so the marker update is slower than in master. nvfortran CPU, one core per rank: 0.91 ms to 13.7 ms (2D, 30 circles + 3 large bodies, 6 cells per diameter), 0.90 ms to 5.8 ms (3D, 60 spheres), medians over 118 and 34 calls. The option is documented as a GPU option; nothing else changes on CPU, but a CPU user of many_ib_patch_parallelism gets a slower marker update.

Tests

./mfc.sh test --no-chemistry --only <UUIDs> on the 93 non-chemistry tests with immersed boundaries (every test whose name contains "IBM", the IBM Example cases included), unmodified branch, goldens as on master:

machine	compiler	GPU build	passed / failed
Tuolumne, MI300A	LLVM flang 24.0.0 (AFAR 24.3.0)	--gpu mp --mpi	93 / 0
A100 node	nvfortran 25.11	--gpu acc --mpi	93 / 0

The same 93 on master: 93 / 0 on the MI300A, 93 / 0 on the A100.

The four new tests, goldens generated on master (175ec9f, MI300A, --gpu mp) so that they also check this branch against master, run on the branch:

machine	compiler	build	passed / failed
Tuolumne, MI300A	LLVM flang 24.0.0 (AFAR 24.3.0)	--gpu mp --mpi	4 / 0
A100 node	nvfortran 25.11	--gpu acc --mpi	4 / 0
NVIDIA workstation CPU	nvfortran 25.11	--no-gpu --mpi	4 / 0

CI covers Frontier (Cray CCE) and GNU on CPU; I did not run either by hand, and I did not run NVHPC 26.x.

Limitations

• The speedup depends on the machine: 2 to 4 times for the call on the A100 and no visible change of the step time when the fluid solve dominates (equal bodies), up to 15 times and 2.7 times of the step on the MI300A, and 2.3 to 470 times when body sizes differ greatly.
• Cost is proportional to (bodies) x (largest box) threads per periodic image, most of which leave immediately. It is large for a case with thousands of bodies and one box hundreds of cells across. It beat master in every case I tried, including the plates; I did not try a box comparable to the whole domain.
• ib_box is allocated on first use with num_local_ibs_max x 27 x 6 integers (19 MB at the current limit of 30 000); it is not allocated unless many_ib_patch_parallelism = T.
• On CPU builds the flag is slower than on master (numbers above).
• Not tested: Cray CCE, GNU, NVHPC 26.x, more than 2 ranks, STL (geometries 5 and 12) and 3D airfoil (11) bodies by hand. The STL and airfoil IBM tests ran with the flag forced on through the IBM test label (marker kernels unchanged for them), and their goldens agree with master's to the same 1e-14.
• Verified with prescribed motion (moving_ibm = 1) and fixed bodies; two-way bodies (moving_ibm = 2) were only in the test-label runs.

This PR was made with Claude Code (AI).

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_ib_patches.fpp	697	+57
src/simulation/m_ibm.fpp	2512	+1

Directory	Lines	Diff
simulation	31231	+58
total	50563	+58

[github-actions]

Claude Code Review

Head SHA: af59396

Files changed:

• 3
• src/simulation/m_ib_patches.fpp
• src/simulation/m_ibm.fpp
• toolchain/mfc/test/cases.py

Findings:
No issues found. Checked in step 4: src/simulation/m_ib_patches.fpp (diff hunks against s_get_bounding_indices and s_get_periodicities on disk), src/simulation/m_ibm.fpp and toolchain/mfc/test/cases.py (diff only).

The prepared context directory was empty. The checkout is master, not the PR head, so I could not read the post-change files. The review is based on the diff plus the master versions of the unchanged callees.

[codecov]

Codecov Report

❌ Patch coverage is 90.00000% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 63.99%. Comparing base (9649d50) to head (af59396).
⚠️ Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_ib_patches.fpp	89.79%	1 Missing and 4 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2214      +/-   ##
==========================================
+ Coverage   63.71%   63.99%   +0.28%     
==========================================
  Files          90       90              
  Lines       24569    24604      +35     
  Branches     3584     3588       +4     
==========================================
+ Hits        15653    15745      +92     
+ Misses       6135     6065      -70     
- Partials     2781     2794      +13     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.