Move Particle Cloud Generation to Pre-Process

[danieljvickers]

Simulation should always read the initial state of particles upon start/restart. The particle cloud code was added haphazardly in an attempt to test scaling, but was poorly designed because it performs CPU-only particle placement that is invarient under restart simultaneously on each processor. This compute would be more efficiently done once on CPUs in pre_process and then read in by simulation each time. It saves some placement time in simulation, allows users to not utilize GPU node time to place particles (or worse, to get an allocation just for the job to fail), and muddies the waters for development on particle clouds.

This PR moves that code to preprocessing, and has preprocessing produce IB state files for the first time step, matching the implementation of all other initial values of the code.

[github-actions]

Claude Code Review

Head SHA: 40af4f1

Files changed:

• 7
• src/pre_process/m_data_output.fpp
• src/pre_process/m_global_parameters.fpp
• src/pre_process/m_particle_cloud.fpp
• src/pre_process/m_start_up.fpp
• src/simulation/m_start_up.fpp
• toolchain/mfc/params/definitions.py
• toolchain/mfc/test/cases.py

Findings:

• Namelist patch_ib initial vel/angular_vel/angles are silently zeroed on every run, not just restarts. s_write_ib_state_0 (src/pre_process/m_data_output.fpp) zero-initializes ib_buf once and, for each of the num_ibs namelist patches, only ever populates fields 17–20 (centroid + radius); fields 8–16 (vel/angular_vel/angles) are written as 0 for every namelist patch. In src/simulation/m_start_up.fpp, the new call call s_read_ib_restart_data(0, particle_cloud_ibs, num_particle_cloud_ibs) is now unconditional — it runs on every start, including a fresh (non-restart) run — and applies s_apply_ib_buf to patch_ib(gbl_id) for every namelist IB, overwriting vel, angular_vel, and angles with the zeros read from ib_state_0.dat. Previously, a fresh start (t_step_start == 0, no cfl_dt restart) never called s_read_ib_restart_data at all, so namelist-set initial kinematics for a moving_ibm patch survived untouched. Now any case that sets a nonzero initial patch_ib(i)%vel, %angular_vel, or %angles in the case file will have that initial condition silently discarded at startup (it only self-corrects on an actual restart, where the second s_read_ib_restart_data(n_start/t_step_start, ...) call re-overwrites with the real saved state). This is a silent-wrong-answer regression, not a crash.
  • failure_scenario: A case defines ib = T with a namelist patch_ib(1) moving immersed boundary given a nonzero initial patch_ib(1)%vel(1) (e.g. a projectile launched with initial velocity). On a fresh run, s_write_ib_state_0 writes zeros for vel into restart_data/ib_state_0.dat, and the new unconditional s_read_ib_restart_data(0, ...) call in simulation overwrites patch_ib(1)%vel back to zero before the first timestep, so the IB starts at rest instead of at the specified initial velocity — no error or warning is produced.

[sbryngelson]

Reviewed e9f2158. The restart-format work is solid - I checked the pre_process writers field-for-field against s_write_parallel_ib_state/s_write_serial_ib_state and they agree on all 20 real(wp) fields, the file_per_process chunk layout, the directory and file names, and the (gbl_patch_id-1)*20*WP_MOK offset. Particle-id contiguity holds in both packers. And the underlying catch is real and valuable: the previous do i = 1, num_ibs read only namelist patches, so restart with particle clouds was silently dropping every particle.

Two things block it, though, and both are silent.

1. Every fresh IB run now zeroes patch_ib%vel, %angular_vel and %angles

s_read_ib_restart_data(0, ...) is now called unconditionally at src/simulation/m_start_up.fpp:924 - previously that path was reached only when t_step_start > 0. The helper applies all 20 fields:

entry%vel         = buf(8:10)
entry%angular_vel = buf(11:13)
entry%angles      = buf(14:16)

But all three pre_process writers do ib_buf = 0._wp and then populate only fields 17-20 (centroid and radius). Fields 8:16 go out as zeros, and the reader copies them over what the case file set.

Those are user-settable namelist parameters (definitions.py:1004-1010), and %angles is consumed right after by s_update_ib_rotation_matrix (m_ib_patches.fpp:427-443, via s_ibm_setup at m_start_up.fpp:933). So, silently:

• examples/2D_ibm_airfoil - angles(3) = -0.5236 - the airfoil runs at 0 degrees AoA instead of 30.
• examples/2D_tumbling_rectangle - vel(1) = -0.05, angles(3) = 0.785, angular_vel(3) = 100 - it neither translates nor tumbles.
• examples/3D_rotating_sphere - angular_vel(3) = 100 - it does not rotate.

Restarts are fine (the later s_read_ib_restart_data(t_step_start, ...) overlays the right kinematics); only fresh starts are wrong. No test sets those fields, and the new test uses num_ibs: 0, so the suite cannot see it.

Either have the pre_process writers fill 8:16 from patch_ib(i), or have the apply helper skip the kinematic fields when t_step == 0.

2. particle_cloud is never broadcast in pre_process

Only rank 0 reads the input file. particle_cloud is a derived-type array, and per .claude/rules/common-pitfalls.md those are not auto-broadcast - fortran_gen.py emits emitters only for lag_params/chem_params/rburn, and only for target == "sim". The simulation's own particle_cloud broadcast (src/simulation/m_mpi_proxy.fpp:200) is hand-written; this PR adds no pre_process counterpart.

num_particle_clouds is auto-broadcast, so on ranks > 0 the loop in s_generate_particle_clouds runs with packing_method == dflt_int and hits case default -> s_mpi_abort("Particle cloud packing method is not a known packing method of MFC"). Every other cloud field is equally garbage there, so even without that abort each rank would generate a different bed.

Repro: ./mfc.sh run <any cloud case>.py -n 2. The pitfalls file asks for a ppn=2 test with exactly this change; the new test at cases.py:1312 is single rank.

3. f_local_rank_owns_location changed meaning during the move

The original in m_collisions.fpp read ib_bc_${X}$%beg; the copy now in src/common/m_helper.fpp:727 reads bc_${X}$%beg. Not the same variable - m_start_up.fpp:1031-1034 saves ib_bc_x = bc_x with the comment "Save original BCs before decomposition overwrites them with MPI neighbor ranks", and m_mpi_common.fpp:1650-1663 does exactly that, so in a periodic decomposed direction bc_x%beg is a neighbour rank id and never equals BC_PERIODIC.

Net effect: the out-of-domain projection now fires where it did not, changing IB ownership at m_collisions.fpp:131/210, m_ibm.fpp:1448 and m_start_up.fpp:1405/1421. An IB drifting past a periodic boundary gets clamped and claimed by the top-edge rank rather than wrapping. ib_bc_* does not exist in pre_process, which is presumably why it was dropped - passing the pre-decomposition BCs in alongside glb_bounds_in would keep both callers honest.

Smaller

• Startup cost. m_start_up.fpp:1306-1314 issues 7 MPI_BCAST calls per particle in a serial loop over the global bed. Placement used to be replicated with zero communication; a 1e5-particle bed is now 7e5 blocking collectives at startup. One packed buffer would do it.
• parallel_io = F + file_per_process = T writes one file and reads another. The writer dispatches on parallel_io first (m_data_output.fpp:756-768), the reader on file_per_process first (m_start_up.fpp:774). Inherited rather than introduced, but this PR makes it fire on every IB run instead of only on restart. Worth a validator rule.
• ib_state_0.dat is now a hard prerequisite for every IB run - a case directory produced before this change aborts on the missing file.
• s_assign_particle_cloud_ib_defaults is now the routine the s_add_cloud_particle pitfall applies to, since it builds the struct that reaches the solver. It sets the same set (all but airfoil_id/model_id/length_*), so no regression - but the two are in different executables now and want a comment saying they must stay in lockstep.
• use m_constants was added to m_helper unqualified for one symbol; only: BC_PERIODIC would avoid widening the import surface of a module every target compiles.

CI

Three failures are yours to clear: 1E6FD465 and 1D4832F8 (the new Particle Cloud -> Box -> Restart cases) have no goldens, and 4E0FBE72 (3D -> IBM -> STL) is out of tolerance, which looks like a real regression rather than a portability issue. I left the goldens alone deliberately - generating them now would bake in the behaviour from finding 1.

[danieljvickers]

@sbryngelson, if this branch passes tests, it will now be ready for review.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_particle_cloud.fpp	0	-411
src/pre_process/m_particle_cloud.fpp	384	+384
src/pre_process/m_data_output.fpp	674	+50
src/common/m_helper.fpp	490	+28
src/simulation/m_start_up.fpp	1231	-27
src/simulation/m_collisions.fpp	312	-25
src/pre_process/m_global_parameters.fpp	471	+21
src/pre_process/m_start_up.fpp	497	+16
src/pre_process/m_mpi_proxy.fpp	142	+12
src/post_process/m_start_up.fpp	789	+3
src/post_process/m_global_parameters.fpp	400	+1
src/pre_process/p_main.f90	19	+1

Directory	Lines	Diff
common	10390	+28
pre_process	4997	+484
simulation	27626	-463
post_process	3403	+4
total	46416	+53

[sbryngelson]

Went through the failures here while triaging CI across the open PRs. They are three different things, not one:

1. Two missing golden files — 2D -> 1 Fluid(s) -> IBM -> Particle Cloud -> Box -> Restart (1E6FD465) and its 2-fluid counterpart (1D4832F8) both fail with The golden file does not exist!. These look like new cases the branch adds without generating goldens for them; ./mfc.sh test --generate --only 1E6FD465 1D4832F8 and committing tests/<UUID>/ should settle both. (#1869 had the same gap for a different case.)

2. One real numerical difference — 3D -> IBM -> STL (4E0FBE72):

Variable n°121225 in D/cons.2.00.000050.dat
  Candidate:  0.00119
  Golden:     0.001189974589
  Error:      abs 2.54E-08, rel 2.14E-05
  Tolerance:  abs 1.00E-10, rel 1.00E-10

Five orders of magnitude outside tolerance, so not roundoff. Given the branch moves particle-cloud generation into pre_process, an initial-condition difference feeding the STL case seems the likelier explanation than a solver change, but that is a guess and worth confirming from the case itself.

3. 3D -> Example -> IGR_jet_1fluid (43B5FEBD) fails to execute rather than mismatching, so it needs its own look.

None of these are the 23 model_eqns=3 failures every branch is currently inheriting from master — those are separate and #1878 clears them.