Skip to content

Euler-Lagrange capabilities - Phase 1 - #1926

Draft
thierrydaoud wants to merge 1 commit into
MFlowCode:masterfrom
thierrydaoud:el-phase1
Draft

thierrydaoud wants to merge 1 commit into
MFlowCode:masterfrom
thierrydaoud:el-phase1

Conversation

@thierrydaoud

Copy link
Copy Markdown

Summary

This PR adds a Lagrangian solid-particle model (particles_lagrange) next to the existing Lagrangian bubble model. Rigid spherical particles are tracked individually. The gas acts on them through quasi-steady drag (Gidaspow, Parmar et al., or Osnes et al.), optional pressure-gradient and added-mass forces, and an optional stochastic drag-fluctuation force. With two-way coupling, the particles act back on the gas through momentum and energy source terms, projected onto the grid with the Gaussian kernel of Maeda & Colonius (2018).

The solver is ported from our research fork (MFC-EL) and adapted to current master (generated parameters, m_eos, the refactored Riemann-state reconstruction). Phase 1 has no particle–particle collisions.

Existing cases are unaffected: every non-particle code path, including the Lagrangian bubble model, behaves as before, and all 27 Lagrangian-bubble regression tests pass unchanged.

Usage and Applications

This code capability is essential to simulate complex flows in the Euler-Lagrange framework. Problems such as a hemi-sphere blast, a planar particle curtain shock, and many other were tested and can be simulated using the developed capability.

What's included

  • Solver: src/simulation/m_particles_EL.fpp (driver, dynamics, sources, RK update, boundaries, MPI handover, I/O) and src/simulation/m_particles_EL_kernels.fpp (GPU seq kernels: Gaussian projection, interpolation, force model, drag correlations). The module is private by default and exports 11 routines, which avoids name clashes with the bubble module that shares its ancestry.
  • Inputs: particles_lagrange, particle_pp%… (physical properties), particle_params%… (solver settings) and lag_voidfrac_wrt. They are registered in toolchain/mfc/params/definitions.py, so the namelists, broadcasts and declarations are generated.
  • Validation: check_particles_lagrange in case_validator.py, with 22 unit tests.
  • Docs: particle section and parameter table in docs/documentation/case.md.
  • Examples: examples/2D_particle_hemisphere (Mach 10 blast through a half-ring of particles) and examples/2D_particle_curtain (planar shock through a dense particle curtain between slip walls). Both generate their particles in case.py.
  • Regression tests: 8 golden tests (alter_particles in cases.py): 2D one-way, two-way, two-way on 2 ranks, Gidaspow, Parmar, and drag fluctuations; 3D two-way on 1 and 2 ranks. The goldens include every particle's position, velocity and force at every step.

Changes to shared code

These are the parts reviewers of other models may want to look at. Bubble call sites keep exactly their previous behavior.

  • Halo width (m_helper_basic.fpp): s_configure_coordinate_bounds takes particles_lagrange alongside bubbles_lagrange, and fd_number is set for particles.
  • Smoothed-field (beta) exchange: beta_vars holds up to num_beta_vars_max = 11 entries, and s_populate_beta_buffers has an optional vars argument. Without it (every bubble call), the list is reset to [1, 2, 5] as before.
  • Halo buffer (m_mpi_common.fpp): enlarged for the particle beta exchange (up to 11 fields over a 2*(mapCells + 1)-deep band). Particles only.
  • MPI particle handover (m_mpi_proxy.fpp): new solid-particle send/receive. The allocation and count exchange that bubbles and particles had in duplicate are factored into s_allocate_particle_comm and s_exchange_particle_counts.
  • Saved Riemann states (m_rhs.fpp): each direction's reconstructed states are copied at the end of s_reconstruct_riemann_states for the particle gradient fields.
  • Parallel-I/O slot for the volume-fraction field (MPI_IO_DATA sized sys_size + 1), and a t = 0 save for particle runs so post-process sees the real volume fraction.
  • Post-process: one lag_name replaces 14 hard-coded lag_bubbles names in the text and Silo writers. The Silo calls now pass the real name length (one existing call passed 16 for an 11-character name).
  • Misc: f_xorshift_rand in m_helper.fpp (the generator removed with m_model in Unify ICPP STL onto the shared IB model path #1546, restored for the drag fluctuations), _emit_struct_bcast in the Fortran generator, and a packer rule so lag_particle output files skip their header line like lag_bubble files.

Verification

Check Result
Particle regression tests, CPU (gfortran 15, macOS) 8/8 pass
Particle regression tests, GPU (LLNL Tuolumne, Cray CCE, OpenACC), after merging master 8/8 pass against the CPU goldens (tolerance 1e-10)
2D Mach 10 hemisphere, 1612 particles, two-way with pressure-gradient and added-mass forces: GPU (Tuolumne) vs CPU flow ≤ 5e-9 relative; particle positions ≤ 5e-18 m, velocities ≤ 1e-13 relative
Same case, one-way variants, GPU vs CPU flow ≈ 1e-9
Same case, full 1009-step run: HiPerGator 8 ranks (GCC 14) vs Mac 1 rank (GCC 15) flow ≤ 1.4e-8 relative, particle positions ≤ 2e-13 m at every save
Rank independence, hemisphere 1 vs 3 ranks bit-for-bit; 8 ranks (cloud split across two ranks) 5e-9
Particles moving between ranks (regression tests) 11 of 20 (2D) and 8 of 20 (3D) particles change ranks; 1 vs 2 ranks agree to 1e-16 (2D) and bit-for-bit (3D, same grid)
Debug vs release build (CPU) agree to 1e-13 (particles) and 1e-9 (flow)
Debug build with bounds checking, multi-rank clean
Full regression suite, HiPerGator CPU (at 26596a11) 721/723 pass; the 2 failures (3D multi-rank) were out-of-memory on the shared node and pass on macOS
Lagrangian-bubble regression tests (after all changes) 27/27 pass, unchanged
./mfc.sh precheck on HiPerGator (Linux), final commit after merging master 7/7 pass
Build: HiPerGator GPU (NVHPC 25.9, OpenACC) compiles; GPU runs were done on Tuolumne

Bugs found and fixed during verification

Running the GPU and debug builds against the CPU reference turned up several bugs, all fixed in this branch:

  1. Gradient weights computed on the GPU with run-time-sized private arrays. In two-way coupling the GPU gave particles about 5 % too slow. The weights are now computed once on the host and copied to the device.
  2. Unset z-component of the fluid velocity in 2D. It entered the added-mass force through a 3-component dot product: zero by luck in CPU release builds, NaN in debug builds (which silently dropped the added-mass force), something else on GPU.
  3. Zero-slip drag. The drag correlations are singular at zero slip; the resulting NaN was caught by a finiteness check that zeroed the whole particle force, including the pressure-gradient and added-mass parts. Drag is now skipped when the slip is zero.
  4. Halo buffer overflow in the multi-rank beta exchange with qs_fluct_force (found with bounds checking).
  5. GPU loop hygiene: four unclosed GPU_PARALLEL_LOOPs, and two scalars missing from private lists.
  6. Particle output on GPU was written from stale host copies; the state is now updated from the device first, in a format that cannot overflow.
  7. Parallel-I/O crash at the first save (MPI_IO_DATA too small for the volume-fraction slot).

Known limitations

These are features the phase-1 solver does not support yet. None of them is a regression, and the validator rejects the unsupported combinations where it can.

  • Periodic and reflective boundaries are rejected by the validator. The beta boundary routines are sized for bubble field counts, and particle wrapping is untested. This is the next PR.
  • Collisions are not included. The particle_pp%*_col constants are accepted but unused until then.
  • Axisymmetric (cyl_coord) has a code path in the projection kernel but has not been tested.
  • One-way coupling with added mass omits the fluid-acceleration term, as in the research code: the stored fluid RHS is only filled in two-way coupling.
  • Restarts and lag_db_wrt require parallel_io = T (the particle restart file); the validator enforces this.
  • Performance has not been profiled. The halo exchange uses the bubble mapCells = 3, although the particle kernel only reaches one cell.
  • OpenMP offload (--gpu mp) has not been verified. An earlier Cray build crashed at startup (Present Table Collision), and it has not been retested on a clean build.
  • Post-process labels: particle mass appears under the bubble label mg, and the per-particle volume fraction reuses the rvel column.

Notes for reviewers

  • The branch is up to date with master (1ececa34, merged in); the particle tests, bubble tests and precheck were rerun after the merge.
  • Commits were made with --no-verify on macOS, where four upstream toolchain tests (test_monitor_*, test_bench_preflight, test_submit_requeue) hang or fail without these changes too. Every precheck step was run manually for each commit, and the full precheck passes on Linux.

AI disclosure

This PR was developed with Claude Code (Anthropic). The model design, the choice of what to port, the test cases and every cluster run (HiPerGator and Tuolumne) were directed and checked by the author.


Contribution Policy

We do not accept pull requests generated primarily by AI without genuine understanding or real-world usage context.

All contributions are expected to demonstrate:

  • A clear understanding of the codebase
  • Alignment with product direction
  • Thoughtful reasoning behind changes
  • Evidence of real-world usage or hands-on experience with the problem

If these expectations are not met, we would prefer to implement the changes ourselves rather than spend time reviewing low-effort submissions.


Acknowledgement

-I confirm this PR meets the above expectations and reflects my own understanding and real-world context.


PR template credit: junegunn

🤖 Generated with Claude Code

@sbryngelson
sbryngelson marked this pull request as draft September 29, 2026 02:52

@sbryngelson sbryngelson left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this. It is a large, carefully tested port, and the verification table and the list of bugs found along the way are genuinely useful. I reviewed CI, the toolchain, the tests and the examples myself, and read the two new Fortran modules and the shared-code changes in full. Several issues need fixing before this can merge; inline comments carry the specifics. Nothing below was built or run locally except where stated; the CI observations come from the job logs of the current run.

Commit attribution. Please rewrite the branch so the commits carry no Co-Authored-By: Claude ... <noreply@anthropic.com> trailers (12 of the 13 non-merge commits have one), and drop the "Generated with Claude Code" footer from the description. We do not list AI tools as MFC contributors. The AI disclosure paragraph is fine to keep. While you are at it, the commits are authored as t.daoud <TDAOUD@MAE-RYX75DJ.local>, a machine-local address, so GitHub does not credit them to your account; squashing into commits authored with your GitHub-linked email fixes both.

CI: 38 failing jobs, three causes

  1. Both new examples fail case validation on every lane: fluid_pp(1)%eos = 'ideal_gas' has no stiffness; do not set fluid_pp(1)%pi_inf. That rule was already in the master merged into this branch, so the examples were not rerun after the merge. They also have no goldens (tests 242025F6, 45F842CA) and are not in the Example skip list, so they will fail again once they validate.
  2. 2D -> Lagrange Particles -> Two-way Coupling -> qs_fluct_force (EDD3540A) fails on NVHPC CPU (e.g. 23.11, 26.3) with a relative error of 8e3 in beta.6. See the inline notes on the seed overflow and the per-stage OU update; either would make this non-portable.
  3. The NVHPC GPU builds crash the compiler on m_particles_EL.fpp (fort2 TERMINATED by signal 11, on 26.3 OpenACC and 25.5 OpenMP on Phoenix). The GPU verification in the description was done with CCE only; NVHPC is our main GPU compiler, so this has to build there.

Blockers (physics and correctness)

  1. Momentum is not conserved in two-way coupling: the source weights are normalized by sum(G V alpha_f) and the RHS divides by alpha_f again, so the fluid receives about F/<alpha_f> for a particle force F (inline).
  2. Energy is not conserved: SE is always zero and the fluid energy source is S.u_f, so the drag dissipation beta|u_f - u_p|^2 disappears instead of heating the gas (inline).
  3. A quiescent cloud (uniform p, u = 0, fixed particles, nonuniform alpha_f) is not an equilibrium: the -(p/alpha_f) grad(alpha_f) momentum term is unbalanced (inline).
  4. fd_order unset or 1 is accepted for particle runs and gives a negative or empty gradient stencil (inline).
  5. One-way coupling with added mass runs with du/dt = 0 because rhs_old is filled only in the two-way path. This is listed under known limitations; the validator should reject the combination rather than run it (inline).
  6. The new t = 0 save also runs on restarts and rewrites the checkpoint being restarted from (inline).
  7. Restarts assign particles by rank index (part_id = proc_particle_counts(proc_rank + 1)), so a restart on a different rank count reads out of range or loses particles, and particle_seed/fqs_fluct are not saved. Either support it or abort when the rank count changes.

Should fix

  • The drag correlation details flagged inline (Ma and Ahmadi coefficient, Osnes Mach floor, undocumented Loth modifications, stochastic force advanced every RK stage).
  • Non-finite forces are silently zeroed (kernels ~648-667), and interpolated rho/p/alpha_f are unbounded near shocks. Please count or abort instead of hiding them.
  • Cost and memory: six full sys_size arrays of reconstructed states are allocated and copied every stage (inline); particle arrays and p_send_ids are sized by nparticles_glb on every rank; the particle state moves between host and device several times per stage. None of this is a problem at 20 particles; all of it is at 1e6.
  • lag_voidfrac_wrt is not restricted to particle runs; on a bubble run it writes the bubble rvel column labelled as void fraction.
  • cyl_coord has a partial code path but no 1/r terms or velocity conversion; please prohibit it in Phase 1.
  • Two small behavior changes to non-particle runs should be mentioned in the description: the nt <= 1 timing fix, and the bubble Silo multimesh name length (16 -> 11).

Cleanliness

  • Registered but never read: particle_pp%cp_particle, ksp_col, nu_col, E_col, cor_col (the tests even set cp_particle). Please register them with the collision PR instead.
  • Dead code: the periodic and reflective wrap paths (forbidden by the validator), particle_in_domain, s_check_celloutside, s_get_cell, the never-written kahan_comp arrays, and the bubble leftovers that mean nothing for rigid particles (particle_draddt, the radius RK update, Rmax/Rmin_stats_part, the stats file). Commented-out code should be deleted (e.g. kernels 826-877, m_particles_EL 1899-1901).
  • About 330 lines of m_particles_EL are verbatim copies of m_bubbles_EL, and s_mpi_sendrecv_solid_particles duplicates the bubble routine's structure. A shared Lagrangian layer would remove most of this.
  • Gidaspow, Parmar et al. and Osnes et al. are not in docs/references.bib; please cite them with the equation numbers each kernel implements.

What would make this excellent

  1. A conservation regression (one particle in a closed box, total momentum and energy to round-off) and a quiescent-cloud regression. Either would have caught blockers 1-3.
  2. Validation of the curtain example against Wagner et al. (2012): Mach 1.66, a 2 mm curtain at 21 percent volume fraction of 115 um glass, 82.7 kPa ambient is that experiment. Plotting the upstream and downstream curtain fronts against the measured ones would be a strong result, and the example should cite it.
  3. Unit tests of the drag correlations against published tables, including the limits (Stokes, M -> 0, phi -> 0) and continuity at every switch. For the record, Parmar 2010 checks out, including continuity at M = 0.6, 1.0 and 1.75, and the Gidaspow dilute limit matches Schiller-Naumann.
  4. A restart regression: 2N steps must equal N steps, a restart, then N more.
  5. Longer term, one Lagrangian container shared with the bubble model (handover, compaction, restart I/O, RK through rk_coef) that stays on the device.

Comment thread src/simulation/m_particles_EL_kernels.fpp Outdated
Comment thread src/simulation/m_particles_EL_kernels.fpp Outdated
Comment thread src/simulation/m_particles_EL.fpp
Comment thread src/simulation/m_particles_EL.fpp
Comment thread src/simulation/m_particles_EL_kernels.fpp Outdated
Comment thread toolchain/mfc/params/definitions.py Outdated
Comment thread toolchain/mfc/case_validator.py Outdated
Comment thread src/post_process/m_data_output.fpp
Comment thread examples/2D_particle_curtain/case.py Outdated
Comment thread examples/2D_particle_curtain/case.py Outdated

@wilfonba wilfonba left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR needs to introduce a shared m_euler_lagrange.fpp to house shared utilities for bubbles and particles. As it stands, there are hundreds of lines of almost verbatim duplicate code, and several hundred more of very similar code.

Comment thread examples/2D_particle_curtain/case.py Outdated
Comment thread examples/2D_particle_curtain/case.py
Comment thread examples/2D_particle_hemisphere/case.py Outdated
Comment thread examples/2D_particle_curtain/README.md
Comment thread examples/2D_particle_hemisphere/README.md
Comment thread src/simulation/m_particles_EL.fpp Outdated
Comment thread src/simulation/m_particles_EL.fpp Outdated
Comment thread src/simulation/m_particles_EL.fpp Outdated
Comment thread src/simulation/m_particles_EL.fpp Outdated
Comment thread src/simulation/m_particles_EL.fpp
@thierrydaoud

thierrydaoud commented Sep 29, 2026 •

Copy link
Copy Markdown
Author

Hello all,

I appreciate you taking the time for this pull request, and provide detailed comments and feedback.

I am going over them and will answer/resolve them one by one.

This might take a day or more, but hopefully will be resolved soon.

After resolving most/all comments, I will push the code changes/commits and squash them here, instead of bugging users with notifications about the various commits. I'll be doing them instead on my local branch

Thanks for your patience and assistance.

@github-actions

github-actions Bot commented Oct 1, 2026

Copy link
Copy Markdown

Lines of Code

File Lines Diff
src/simulation/m_particles_EL.fpp 1236 +1236
src/simulation/m_particles_EL_kernels.fpp 679 +679
src/simulation/m_euler_lagrange.fpp 391 +391
src/simulation/m_bubbles_EL.fpp 1314 -325
src/simulation/m_mpi_proxy.fpp 718 +186
src/common/m_derived_types.fpp 500 +21
src/simulation/m_global_parameters.fpp 811 +18
src/simulation/m_bubbles_EL_kernels.fpp 502 -14
src/simulation/m_rhs.fpp 1980 +12
src/common/m_helper.fpp 499 +9
src/common/m_mpi_common.fpp 1493 +9
src/common/m_boundary_common.fpp 455 +7
src/simulation/m_start_up.fpp 1251 +7
src/post_process/m_data_output.fpp 1262 +5
src/simulation/p_main.fpp 79 +5
src/simulation/m_time_steppers.fpp 883 +2
src/common/m_global_parameters_common.fpp 251 +1
src/post_process/m_global_parameters.fpp 415 +1
Directory Lines Diff
common 10463 +47
simulation 30192 +2197
post_process 3505 +6
total 49189 +2250

Lagrangian solid-particle model (particles_lagrange) next to the Lagrangian
bubbles. Rigid spheres are tracked individually; the gas acts on them through
quasi-steady drag (Gidaspow, Parmar et al., Osnes et al.), optional
pressure-gradient and added-mass forces and an optional stochastic drag
fluctuation. With two-way coupling the particles act back on the gas through
momentum and energy sources projected with a Gaussian kernel. No collisions in
phase 1.

- Solver: src/simulation/m_particles_EL.fpp and m_particles_EL_kernels.fpp.
- Shared module m_euler_lagrange.fpp for bubbles and particles: cell location,
  domain and cell-volume helpers, input parsing and start-up, void-fraction and
  evolution files, MPI-IO restart read/write. Bubble results are unchanged.
- Coupling: force deposit normalized by sum(G V); the gas receives the particle
  momentum change m du/dt (including the added-mass reaction) and the work
  -F.u_p; pressure terms -(alpha_p/alpha_f) grad p and
  -(alpha_p/alpha_f) div(p u) as in the bubble solver, so a quiescent cloud
  stays at rest.
- Inputs registered in toolchain/mfc/params/definitions.py, validated by
  check_particles_lagrange (unsupported combinations rejected), documented in
  docs/documentation/case.md with references in docs/references.bib.
- Examples 2D_particle_curtain and 2D_particle_hemisphere with goldens; 12
  particle regression tests, including a quiescent cloud and a closed box.
- Shared-code changes: halo width and buffers, beta exchange, solid-particle
  MPI handover, particle gradients inside the RHS direction loop, parallel-I/O
  volume-fraction slot, post-process naming and rank-count handling.
@thierrydaoud

Copy link
Copy Markdown
Author

Summary

This PR adds a Lagrangian solid-particle model (particles_lagrange) next to the existing Lagrangian bubble model. Rigid spherical particles are tracked individually. The gas acts on them through quasi-steady drag (Gidaspow, Parmar et al., or Osnes et al.), optional pressure-gradient and added-mass forces, and an optional stochastic drag-fluctuation force. With two-way coupling, the particles act back on the gas through momentum and energy source terms, projected onto the grid with a Gaussian kernel.

The solver is ported from our research fork (MFC-EL) and adapted to current master. Phase 1 has no particle–particle collisions. Utilities the bubble and particle solvers had in duplicate now live in a shared module, m_euler_lagrange.

Existing cases are unaffected: every Lagrangian-bubble golden passes unchanged, bit for bit.

Usage and Applications

The model simulates particle-laden compressible flows in the Euler–Lagrange framework, such as shock–particle-curtain interaction and blast waves through particle clouds. Two examples are included and run in CI.

What's included

  • Solver: src/simulation/m_particles_EL.fpp (driver, dynamics, sources, RK update, boundaries, I/O) and m_particles_EL_kernels.fpp (GPU seq kernels: Gaussian projection, interpolation, force model, drag correlations).
  • Shared Euler–Lagrange module: src/simulation/m_euler_lagrange.fpp, used by bubbles and particles: cell location, domain and cell-volume helpers, input-file parsing, start-up (restart point, rank bounds, output directory), the void-fraction and evolution files, and the MPI-IO restart read/write. Each solver keeps only what touches its own state.
  • Inputs: particles_lagrange, particle_pp%rho0ref_particle, particle_params%… and lag_voidfrac_wrt, registered in toolchain/mfc/params/definitions.py.
  • Validation: check_particles_lagrange in case_validator.py, with 19 unit tests. It rejects the combinations phase 1 does not support (see Known limitations).
  • Docs: particle section and parameter table in docs/documentation/case.md; drag-law and model references in docs/references.bib.
  • Examples: examples/2D_particle_curtain (Mach 1.66 shock through a 2 mm glass-particle curtain; the setup of Wagner et al., 2012) and examples/2D_particle_hemisphere (Mach 10 blast through a half-ring of particles). Shock states come from a Mach number via the normal-shock relations; particles come from a separate gen_particles.py that case.py calls only when ./mfc.sh run runs pre_process. READMEs show result images.
  • Regression tests: 12 golden tests: 2D one-way, two-way, two-way on 2 ranks, Gidaspow, Parmar, drag fluctuations, a quiescent cloud and a closed box; 3D two-way on 1 and 2 ranks; and the two examples.

Changes to shared code

  • Shared module m_euler_lagrange (above). The bubble module now uses it; bubble results are unchanged. For bubbles this also adds error checks on every restart-file open, an abort on a rank-count mismatch at restart (it used to read the wrong part of the file), and appending to the evolution file on a restart.
  • Particle field gradients (m_rhs.fpp): computed inside the RHS direction loop from that sweep's reconstructed face states, only when a gradient force is on. No extra copies of the face states are kept.
  • Halo width (m_helper_basic.fpp): s_configure_coordinate_bounds takes particles_lagrange alongside bubbles_lagrange.
  • Smoothed-field (beta) exchange: beta_vars holds up to 11 entries and s_populate_beta_buffers has an optional vars argument; bubble calls are unchanged.
  • Halo buffer (m_mpi_common.fpp): sized for the particle beta exchange, with the bubble case included.
  • MPI particle handover (m_mpi_proxy.fpp): solid-particle send/receive. Allocation and count exchange are shared with bubbles (s_allocate_particle_comm, s_exchange_particle_counts).
  • Parallel I/O slot for the volume-fraction field, and a t = 0 save on fresh particle starts so post-process sees the initial volume fraction.
  • Post-process: one lag_name replaces 14 hard-coded lag_bubbles names, and particles are written for any post-process rank count.
  • Misc: f_xorshift_rand in m_helper.fpp (for the drag fluctuations), _emit_struct_bcast in the Fortran generator, a packer rule for lag_particle output headers, and the Example test harness (particle examples run their own generator; lag_db_wrt is off and t_stop is capped at 5e-5 in that suite).

Physics changes made during review

These change two-way results; all particle goldens were regenerated.

  1. Momentum exchange. The force deposit was normalized by Σ(G V α_f) and then divided by α_f again, so the gas received F/⟨α_f⟩ (1.27 F at the curtain's packing). It is now normalized by Σ(G V).
  2. Added-mass reaction. The particle equation is solved implicitly, (m + m_a) du/dt = f_p, but the gas received −f_p and missed the reaction to the −m_a du/dt part. It now receives −m du/dt.
  3. Energy exchange. The energy source was zero and the gas received S·u_f, so drag dissipation left the system. The gas now receives −F·u_p.
  4. Pressure terms. −(p/α_f)∇α_f and −(p/α_f)u·∇α_f (the expansion of a ∇·(α_f p) form whose balancing interphase term was missing) are replaced by the bubble solver's −(α_p/α_f)∇p and −(α_p/α_f)∇·(p u). A cloud at rest in uniform pressure now stays exactly at rest; before, the same setup reached 22% pressure errors.
  5. Drag floors. The volume-fraction floors and the Gidaspow and Osnes Reynolds floors are removed. The Gidaspow floor cut the Stokes drag by Re/0.1 below Re = 0.1. Parmar keeps its floor, which bounds its log(Re) fits.
  6. Drag-fluctuation model. Advanced once per time step with the exact Ornstein–Uhlenbeck update, a portable seed, and corrected Ma–Ahmadi coefficients.

In a closed box with slip walls (moving particle cloud, gas at rest), total momentum drift went from +14.9% to +1.9% after 100 steps, and from +25.7% to +2.4% after 1000 steps in a longer box, where it no longer grows once the exchange is over. The energy lost early in the exchange went from 78% of the particles' kinetic energy to a few percent. The remainder comes from discretizing the α_f terms and shrinks as the α_f field gets smoother.

Verification

Check Result
Particle goldens (12), CPU, gfortran 15 (macOS) 12/12 pass
Lagrangian-bubble goldens (28), CPU, after the shared-module refactor and the master merge 28/28 pass, bit for bit
Particle (12) and Lagrangian-bubble (28) goldens, GPU, NVHPC 25.9 OpenACC (HiPerGator), at 7b829d67 40/40 pass against the CPU goldens
Particle restart: 20 steps straight vs 10 + restart + 10 flow 4e-16; positions, seeds and drag-fluctuation state identical
Bubble restart (2D two-way, 2 ranks), same check flow 8e-16, bubble state 3e-18
Quiescent cloud (fixed particles, uniform pressure) gas stays exactly at rest

Runs before the review's physics fixes, kept for reference: GPU vs CPU on LLNL Tuolumne (Cray CCE, OpenACC) agreed to 5e-9 in the flow on the hemisphere case; NVHPC 25.9 OpenACC on a HiPerGator GPU node passed the 8 particle tests then present; 1 vs 2 ranks agreed to 1e-16 (2D) and bit for bit (3D) with particles crossing ranks; and a debug build with bounds checking ran clean on multiple ranks.

Bugs found and fixed

During the port: gradient weights computed on the GPU with run-time-sized private arrays (particles about 5% slow on GPU); an unset z-velocity entering the added-mass force in 2D; zero-slip drag NaNs that zeroed the whole particle force; a halo-buffer overflow in the multi-rank beta exchange; unclosed GPU loops and missing private scalars; particle output from stale host copies; and a parallel-I/O crash at the first save.

During review: the four coupling errors above; interpolation overshoot (the barycentric interpolant is now clamped to its stencil range, which removed negative interpolated pressures at blast fronts); a non-finite force now aborts with a report of the particle, its Reynolds and Mach numbers and the force term responsible, instead of being zeroed silently; seed overflow (undefined behaviour in Fortran); restart files that lost the drag-fluctuation state and the seed; and post-process missing particles when run on a different rank count.

Known limitations

The validator rejects the unsupported combinations.

  • Collisions are not included (phase 2).
  • Periodic and reflective boundaries, and axisymmetric coordinates, are rejected.
  • One-way coupling with added mass is rejected: the fluid acceleration it needs is only available with two-way coupling.
  • Restarts need parallel_io = T and the same number of ranks as the run that wrote the file (inherited from the bubble restart format).
  • Added mass across a restart: the fluid acceleration from the previous stage is not saved, so the first stage after a restart uses zero (a difference of about 5e-5 in the tests).
  • Multi-fluid added mass: dρ/dt uses fluid 1's partial density only.
  • Shared container: adding, copying, the RK update, boundary enforcement and the MPI send/receive packing are still per solver. Sharing them needs one Lagrangian state container for both, which is planned with the collision work.
  • Grid-aligned artifact: with HLLC, the hemisphere blast shows a carbuncle along the x = 0 grid line; the example uses HLL.
  • Performance has not been profiled.

Acknowledgement

  • I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

@thierrydaoud
thierrydaoud marked this pull request as ready for review October 2, 2026 03:50
@codecov

codecov Bot commented Oct 2, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 65.88670% with 554 lines in your changes missing coverage. Please review.
✅ Project coverage is 63.83%. Comparing base (c7feed2) to head (2c148cd).

Files with missing lines Patch % Lines
src/simulation/m_particles_EL.fpp 63.77% 202 Missing and 57 partials ⚠️
src/simulation/m_euler_lagrange.fpp 48.31% 97 Missing and 26 partials ⚠️
src/simulation/m_particles_EL_kernels.fpp 79.45% 28 Missing and 48 partials ⚠️
src/simulation/m_bubbles_EL.fpp 49.25% 30 Missing and 4 partials ⚠️
src/simulation/m_mpi_proxy.fpp 78.57% 18 Missing and 12 partials ⚠️
src/post_process/m_data_output.fpp 36.36% 13 Missing and 1 partial ⚠️
src/pre_process/m_data_output.fpp 0.00% 1 Missing and 2 partials ⚠️
src/simulation/m_rhs.fpp 70.00% 0 Missing and 3 partials ⚠️
src/common/m_boundary_common.fpp 60.00% 0 Missing and 2 partials ⚠️
src/common/m_helper_basic.fpp 0.00% 0 Missing and 1 partial ⚠️
... and 9 more
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1926      +/-   ##
==========================================
+ Coverage   62.80%   63.83%   +1.02%     
==========================================
  Files          86       89       +3     
  Lines       22385    23649    +1264     
  Branches     3304     3457     +153     
==========================================
+ Hits        14060    15096    +1036     
- Misses       6073     6166      +93     
- Partials     2252     2387     +135     

☔ 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.

@sbryngelson
sbryngelson marked this pull request as draft October 2, 2026 16:28

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

4 participants