1. Per-Step Particle Pipeline

Order inside the coupled step, as implemented in AdvanceSimulation (src/runloop.c):

refresh Eulerian transport properties when a turbulence model is active: ComputeEulerianDiffusivity and ComputeEulerianDiffusivityGradient
advect particle positions using the velocity interpolated at the previous step: UpdateAllParticlePositions
settle particles - locate new host cells and migrate across ranks: LocateAllParticlesInGrid
remove particles that are now lost or outside the global domain: CheckAndRemoveLostParticles
interpolate the newly computed Eulerian velocity onto the settled positions: InterpolateAllFieldsToSwarm
update particle physics fields: UpdateAllParticleFields
scatter particle fields back to Eulerian storage: CalculateParticleCountPerCell and ScatterAllParticleFieldsToEulerFields

The advect-then-interpolate ordering is deliberate and is the single most important property of this loop. A particle carries the state interpolated at \(t_n\) into step 2, and UpdateParticlePosition displaces it by three contributions - convective velocity, drift along the diffusivity gradient, and a stochastic Brownian kick:

\[ \mathbf{X}^{n+1} = \mathbf{X}^{n} + \left( \mathbf{u}_p^{\,n} + \nabla D^{\,n} \right) \Delta t + \Delta \mathbf{X}_{\mathrm{brownian}}, \]

where \(\nabla D\) is the interpolated diffusivity gradient carried on the swarm and \(\Delta\mathbf{X}_{\mathrm{brownian}}\) comes from CalculateBrownianDisplacement, scaled by the particle's interpolated diffusivity. Step 5 then supplies the state at \(t_{n+1}\) from the flow field FlowSolver has just computed, ready for the next step. Interpolation therefore closes each step rather than opening it.

Particles are not advected without a valid interpolated state: both PerformInitializedParticleSetup and PerformLoadedParticleSetup settle the swarm and run one interpolation at \(t=0\) before the step loop begins, so the first physical advection already has interpolated velocity, diffusivity, and diffusivity gradient available.

Each stage relies on valid DMSwarm_CellID, interpolation weights, and synchronized ghost data.

Step 5 dispatches through InterpolateEulerFieldToSwarm, which selects between two interpolation paths based on SimCtx.interpolationMethod:

Trilinear (default): direct interpolation from the 8 nearest cell centers, second-order on curvilinear grids.
CornerAveraged (legacy): two-stage center-to-corner averaging then trilinear from corners.

See Trilinear Interpolation and Particle-Grid Projection for method details and configuration.

2. Core Particle Fields In Use

Commonly used swarm fields include:

position (position)
velocity (velocity)
cell ID (DMSwarm_CellID)
location status (DMSwarm_location_status)
diffusivity and scalar (Diffusivity, Psi)

Status transitions (NEEDS_LOCATION, ACTIVE_AND_LOCATED, MIGRATING_OUT, LOST) determine whether particles are walked, migrated, or skipped in settlement passes.

3. Current Built-In Physics

Current scalar model path:

This presently implements IEM-style relaxation for Psi, with the mixing constant C_IEM set by scalar_transport.iem_constant (default 2.0; 10. scalar_transport). Psi starts at zero unless models.physics.particles.fields gives it an initial value (1.1 Configured Particle Values (fields)); without one, the relaxation has nothing to act on (1. IEM Mixing Update In Current Code).

4. Statistics and Diagnostics

Post statistics currently include global kernels such as:

ComputeParticleMSD

Additional health indicators are available from migration counters and settlement-pass counts stored in SimCtx fields updated by location logic. For particle-enabled runs, the runtime also writes <run.analysis.metrics>/search_metrics.csv with per-step and cumulative loss plus aggregated search-attempt, traversal, tie-break, boundary-clamp, bbox-guess, and pass-depth signals. The examples/search_robustness/ bundle is the dedicated end-to-end reference for interpreting that artifact.

5. Limitations

Particles are not wrapped across periodic boundaries. Neither the motion update nor the walking search handles periodicity, so a particle that crosses a periodic face leaves the domain and is removed as lost.
A lost particle is removed, not reflected. Steady loss through a boundary you did not intend usually means the carrier velocity or the diffusion step carries particles out; the per-step lost count reports it.
Diffusion, advection, the gradient drift and interpolation are verified. A 10,000-particle cloud in zero flow reproduced the Einstein relation to 0.17% in the MSD slope; a cloud in uniform flow drifted at exactly the carrier velocity; paired runs with the same seed resolved the drift a linear diffusivity gradient adds to 1.1% of a t (diffusivity-gradient-drift-paired-2026-09-18); and Trilinear interpolation of the TGV3D field converged at order 1.97. Only one gradient, and zero or uniform carriers, were measured. The examples are listed in 3. Particle Verification Family.
Brownian trajectories repeat only for the same rank count. Placement and Brownian draws are seeded from models.physics.particles.random_seed (default 12345) plus the MPI rank, so identical inputs on the same number of ranks give identical trajectories, while a different rank count gives a different, equally valid, realization. A restart does not continue the sequence it stopped in: it draws a fresh sequence determined by the seed and the restart step.

6. Extending To New Closures

Recommended extension pattern:

add particle model selector to YAML/flags,
branch inside particle physics orchestrator,
isolate model-specific kernels from motion/location kernels,
keep scatter/interpolation interfaces unchanged,
validate with smoke tests and deterministic small cases.

For configuration contract changes, update:

Table of Contents