Moving EL Bubbles with MPI Decomposition

[wilfonba]

Description

This PR adds support for moving Euler-Lagrange bubbles and various force models. Tracer particles that simply follow the background flow are implemented, as are Newton's 2nd Law particles that respond to body forces, pressure gradients, and drag. This PR is marked as a draft for now as there are several things that I need to complete before a proper merge can occur, but I wanted these changes to be visible to other contributors working on related features.

Fixes #(issue)

Type of change

• New feature

Testing

16 rank CPU versus 16 rank GPU

This test was ran using the examples/3D_lagrange_bubblescreen test case and compares the results for a 16 rank CPU simulation to a 16 rank GPU simulation. The visualization shows the following things:

1. The difference in pressure in the vertical plane, with darker colors corresponding to larger errors. There seems to be a reflection of some sort from the subsonic buffer outflow condition, and the acoustic source introduces a small amount of difference, though neither of these features was modified by this PR. (I intend to add a video with no bubbles in the domain to see if these artifacts still exist)
2. The difference in void fraction in the horizontal plane. The difference is 1e-38 everywhere in ParaView, which corresponds to bitwise identical results.
3. The bubbles for each simulation overlayed on top of each other. The bubbles from the CPU result are colored in blue while the bubbles from the GPU result are colored in white. The small bit of blue that is observed in the video is a result of the rendering spheres with a fixed number of facets.
4. The pressure field in the lower front quadrant as a surface render with darker colors indicating higher pressure.
5. The MPI domain decomposition as a black lattice.

test.mp4

Checklist

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

AI code reviews

Reviews are not triggered automatically. To request a review, comment on the PR:

• @coderabbitai review — incremental review (new changes only)
• @coderabbitai full review — full review from scratch
• /review — Qodo review
• /improve — Qodo code suggestions

[github-actions]

Claude Code Review

Head SHA: 54d388d
Files changed: 47 (+5320 / -1238)

Changed files (key):

• src/simulation/m_bubbles_EL.fpp — Moving EL bubble dynamics, MPI migration, BC enforcement
• src/simulation/m_bubbles_EL_kernels.fpp — Cell-list smearing, Lagrange interpolation, force models
• src/simulation/m_mpi_proxy.fpp — Particle MPI send/recv infrastructure
• src/common/m_mpi_common.fpp — MPI reductions, beta-variable buffer exchange
• src/common/m_boundary_common.fpp — Boundary condition handling for beta fields
• src/simulation/m_bubbles.fpp — elemental → GPU_ROUTINE for device compatibility
• src/simulation/m_sim_helpers.fpp — max_dt array → scalar + GPU reduction
• src/pre_process/m_start_up.fpp — Namelist additions

Summary:

• Adds Newton's 2nd Law (vel_model 2/3) and tracer (vel_model 1) particle motion with drag, pressure gradient, and gravity forces
• Replaces GPU-atomic smearing with a cell-list approach (eliminates non-determinism for Gaussian/delta kernels)
• Adds full MPI particle migration: bubbles crossing subdomain boundaries are packed/unpacked and transferred to neighbor ranks
• Switches adap_dt_stop_max (MAX reduction) to adap_dt_stop_sum (GPU atomic add) in both EE and EL paths
• Simplifies max_dt from a per-cell array to a scalar via GPU loop reduction

---

Findings

🔴 Bug: Wrong argument passed to f_interpolate_velocity in f_advance_step (adap_dt + vel_model 1)

File: src/simulation/m_bubbles.fpp, inside f_advance_step, vel_model case(1)

vTemp(l) = f_interpolate_velocity(fR, cell, l, q_prim_vf)  ! ← fR is bubble RADIUS

f_interpolate_velocity(pos, cell, i, q_prim_vf) expects a spatial coordinate pos. fR is the current bubble radius, not a position. The non-adap_dt path (in s_compute_bubble_EL_dynamics) correctly passes myPos(l). This will silently produce wrong velocities for tracer particles when adap_dt = .true..

Fix: Replace fR with fPos(l):

vTemp(l) = f_interpolate_velocity(fPos(l), cell, l, q_prim_vf)

---

🔴 Bug: Wrong denominator in 4th-order Lagrange basis function L(1)

File: src/simulation/m_bubbles_EL_kernels.fpp, f_interpolate_velocity (fd_order == 4)

L(1) = ((pos - xi(2))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/ &
       ((xi(1) - xi(2))*(xi(1) - xi(3))*(xi(1) - xi(4))*(xi(2) - xi(5)))  ! ← last factor WRONG

The Lagrange basis requires (xi(1) - xi(5)) in the last factor, not (xi(2) - xi(5)).

Fix:

L(1) = ((pos - xi(2))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/ &
       ((xi(1) - xi(2))*(xi(1) - xi(3))*(xi(1) - xi(4))*(xi(1) - xi(5)))

---

🟡 Issue: Duplicate hyper_cleaning in pre_process namelist

File: src/pre_process/m_start_up.fpp, around line 148

g0_ic, p0_ic, hyper_cleaning, hyper_cleaning  ! ← 'hyper_cleaning' listed twice

Fortran namelist variables must be unique. Some compilers will silently ignore the duplicate; others will raise a compile-time or runtime error. This should be deduplicated.

---

🟡 Issue: Optional arguments in f_advance_step used without present() guards

File: src/simulation/m_bubbles.fpp

f_advance_step declares fPos, fVel, fRe, cell, q_prim_vf as optional. Inside the select case (lag_vel_model) block these are dereferenced without any if (present(...)) guard. The EE path calls the function without these optionals and relies on lag_vel_model == dflt_int reaching the case default branch. This is safe in practice but fragile — any future change to lag_vel_model initialization could cause an undefined-behavior dereference of absent optionals. Consider adding a GPU_ROUTINE-compatible guard (e.g., pass a sentinel moving_lag_bubbles logical) or use non-optional arguments.

---

🟡 Issue: Missing @:DEALLOCATE for MPI particle buffers

File: src/simulation/m_mpi_proxy.fpp

p_send_buff, p_recv_buff, and p_send_ids are allocated in s_initialize_particles_mpi but have no matching @:DEALLOCATE in any finalization routine (s_finalize_mpi_proxy_module or equivalent). Per project convention every @:ALLOCATE requires a matching @:DEALLOCATE. This is a memory leak for long or repeated runs (though the OS reclaims memory on exit).

---

🟢 Improvements (non-blocking)

• src/simulation/m_bubbles_EL.fpp line ~1987: 3*fV*fVel(l)/fR — bare integer 3 in floating-point arithmetic; prefer 3._wp for consistency with coding conventions (though standard Fortran promotes correctly).
• The Kahan summation introduced in delta/Gaussian smearing is a good accuracy improvement for bubble projection; the cell-list restructuring correctly eliminates atomic-write non-determinism.
• The max_dt scalar + GPU min reduction is cleaner than the previous array approach.

---

This is a draft PR; findings above are intended to help prioritize before merge. The GPU vs CPU comparison video is compelling evidence that the MPI decomposition is correct.