Parcel-in-Cell Simulations

• Parcel-in-Cell simulations are Eulerian–Lagrangian methods that couple a mesh-based continuum with moving parcels carrying discrete state information.
• They enable efficient representation of multiphase flows, mantle convection, and plasma behavior by aggregating numerous similar physical particles into parcels.
• The approach features adaptive mesh refinement, two-way coupling, and scalable parallel execution with minimal particle loss.

Searching arXiv for the cited works to ground the article and confirm bibliographic details. Searching arXiv for "(Joshi, 2016) particle in cell code development high current ion beam transport plasma simulations". Parcel-in-Cell simulations are Eulerian–Lagrangian methods in which a mesh-based continuum discretization is coupled to moving parcels or particles that carry discrete state. In the broad particle-in-cell framing used for modern computational science, these Lagrangian entities may be passive tracers, material-interface trackers, carriers of history variables such as strain or damage, or weighted representatives of many physical particles. In CFD-oriented parcel formulations, a parcel is a computational representative of multiple similar physical particles that moves through the mesh, samples Eulerian fields, and deposits source terms back to cells. This usage is broader than electrostatic or electromagnetic plasma PIC, but it preserves the same central architecture: fields or continuum variables are represented on a mesh, particle or parcel state is advanced in physical space, and information is exchanged in both directions as required by the model (Gassmoeller et al., 2016, Marzouk et al., 2024, Vay, 2020).

1. Concept and scope

The broad definition adopted in large-scale finite-element work is that particle-in-cell methods “couple discrete (Lagrangian) particles and (Eulerian) mesh-based continuum methods.” In that sense, parcel-in-cell simulations include tracer transport, material-interface tracking, particles carrying constitutive history, and two-way coupled schemes in which particle-carried quantities are projected back onto quadrature points, cells, or field spaces. The reference implementation described in ASPECT on top of deal.II was motivated by mantle convection and adaptive unstructured finite-element meshes, but the algorithms were presented as generic for continuum solvers with dynamic mesh refinement and repartitioning (Gassmoeller et al., 2016).

In dispersed multiphase CFD, the same architecture appears in a more explicitly weighted form. The gas phase is solved on an Eulerian mesh, the dispersed phase is represented in Lagrangian form, and gas–particle interaction is handled through two-way momentum exchange. The “parcel” terminology emphasizes that one computational entity can stand for many physical particles with similar properties. The swirling-pipe study of Sommerfeld and Qiu’s experiment is explicit on this point: multiple similar particles are grouped in a single parcel, thereby reducing the amount of computation, while the parcel still obeys the same motion equation as an individual tracked particle (Marzouk et al., 2024).

This broader usage differs from plasma-physics PIC only in scope, not in principle. In electrostatic and electromagnetic plasma PIC, the discrete entities are macroparticles representing distribution functions, while the mesh carries electrostatic potential or electromagnetic fields. In parcel-in-cell formulations for continuum mechanics or two-phase flow, the parcels need not represent charged plasma particles, but they still participate in the same interpolate–advance–deposit cycle (Joshi, 2016, Vay, 2020).

2. Fundamental algorithmic structure

A generic parcel-in-cell workflow on adaptive finite-element meshes consists of particle generation, storage by containing cell, advection by a mesh-based velocity field, property update from sampled Eulerian fields, relocation when a particle changes cells, migration across MPI ranks when it crosses subdomain boundaries, transfer during refinement/coarsening or repartitioning, and optional reconstruction of particle-carried quantities back onto quadrature points, cells, or fields. In the notation of the finite-element formulation, particle trajectories satisfy

$$\frac{d}{dt}\mathbf x_k(t)=\mathbf u(\mathbf x_k(t),t),$$

or, more precisely, are advanced under the discrete field interpolation actually available from the mesh solver (Gassmoeller et al., 2016).

In CFD parcel simulations, the same loop appears in a force-based form. The gas phase satisfies the unsteady Favre-averaged Navier–Stokes equations, with a source term $S_p$ representing particle feedback in the momentum equation. The dispersed phase is advanced from

$$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$

where $U^*=U+U'$ combines interpolated resolved velocity and a stochastic turbulent fluctuation, and $\tau$ is the particle momentum relaxation time implied by the drag law. The practical cycle is to solve the Eulerian gas equations, interpolate gas velocity and turbulence quantities to parcel locations, sample or update the turbulent fluctuation, advance each parcel, and accumulate momentum exchange into crossed cells using parcel residence time (Marzouk et al., 2024).

Classical plasma PIC implements the same structure with field equations specialized to electrostatics or electromagnetics. In the Frankfurt electrostatic beam-transport code, the implemented loop is: initialize particles in phase space; deposit particle charge onto grid points; solve Poisson’s equation for the electrostatic potential; compute the electric field from the potential; interpolate fields back to particle positions; advance particles in time under electric and magnetic forces; and repeat with appropriate boundary and source/loss treatment. Jean-Luc Vay’s plasma-accelerator review states the corresponding electromagnetic cycle as particle push, charge or current deposition, field advance by Maxwell’s equations, and field gather, typically in a leapfrog arrangement (Joshi, 2016, Vay, 2020).

A plausible implication is that “parcel-in-cell” is best understood as a structural family of algorithms rather than a single physical model. What remains invariant across applications is the coupling of Lagrangian carriers to an Eulerian discretization through local interpolation and deposition.

3. Parcel representation and state variables

The large-scale finite-element formulation stores each particle with a unique global ID, global position, reference-cell coordinates, and a pointer to a run-time-sized array of scalar properties. The dynamic properties array is deliberate: the number of particle properties is a run-time constant, not a compile-time one, and vector or tensor state can be flattened into scalars. This allows the same particle container to carry composition fractions, age, accumulated strain, deformation-gradient components, damage, temperature history, or custom chemistry state without changing the underlying storage type (Gassmoeller et al., 2016).

A central implementation decision is to store locally owned particles in $\tau$3 that is, grouped by containing cell. This gives efficient cell-local iteration and makes relocation a delete/reinsert operation. The same work recommends a single particle schema for a given simulation together with a property manager that stores semantics and name-to-offset mappings. This avoids heterogeneous particle classes, simplifies serialization and memory pooling, and supports runtime selection of generators, integrators, and projection operators (Gassmoeller et al., 2016).

Reference-cell coordinates are especially important on finite-element meshes. If a particle lies in cell $K$, the stored coordinate

$\hat{\mathbf x}_k=\Phi_K^{-1}(\mathbf x_k)$

avoids repeated inversion of the physical-to-reference mapping during field interpolation. Because shape-function evaluation is naturally performed in reference space, caching reference_location is a substantial practical optimization for FE-based parcel methods (Gassmoeller et al., 2016).

In CFD parcel formulations, the parcel state is simpler but weighted. A parcel carries position, velocity, sampled particle diameter, and total represented mass or multiplicity. The swirling-flow study parameterized the weighting by particles-per-parcel, $\mathrm{ppp}=1,\ 10,\ 100$. When $\mathrm{ppp}=1$, the scheme is an individual-particle approach; when $\mathrm{ppp}>1$, each computational trajectory represents multiple physical particles assumed to share position, velocity, diameter, trajectory, and force history. The implementation notes that because diameters are sampled from a log-normal distribution, one cannot simultaneously impose exactly constant particle mass flow rate and constant particles-per-parcel in the OpenFOAM injection model (Marzouk et al., 2024).

4. Transport, cell crossing, and large-scale parallel execution

Finite-element parcel advection is constrained as much by field interpolation and mesh search as by the ODE integrator itself. The implemented integrators are explicit Forward Euler, RK2, and RK4, but the associated error analysis shows that increasing particle-integrator order beyond the accuracy of the underlying Eulerian field approximation often yields limited benefit. The particle-trajectory error relative to the exact velocity field is bounded by terms involving the particle time step, the spatial error $S_p$0, and the temporal interpolation error $S_p$1 (Gassmoeller et al., 2016).

After every particle step, or every substep for RK2 and RK4, the code checks whether the particle remains in its current cell by updating reference coordinates and testing whether $S_p$2. A global search over all local cells would be too expensive, so the implementation relies on a CFL-like condition: if CFL $S_p$3, a particle generally moves by at most one cell diameter per step. The search therefore first inspects neighboring cells of the old cell and only rarely falls back to a full search. In the common case, relocation is therefore $S_p$4 per moved particle. Candidate neighbors are further ordered by decreasing

$S_p$5

where $S_p$6 is the particle displacement relative to the closest old-cell vertex and $S_p$7 points from that vertex to the candidate-cell center; this directional heuristic often finds the correct new cell on the first try (Gassmoeller et al., 2016).

Parallel migration is neighbor-only under the same CFL assumption. If the new cell is locally owned, the particle is reinserted locally; if it is a ghost cell owned by another process, serialized particle data and destination-cell information are appended to a send buffer and transmitted by nonblocking point-to-point MPI; otherwise the particle is discarded. The discard path covers rare transitions into refined descendants not known locally or advection outside the domain. The reported losses quantify the tradeoff: with RK2 at CFL 1, less than $S_p$8 of particles per 1000 time steps are lost; Forward Euler at the same CFL loses about $S_p$9 per 1000 steps; reducing the time step to CFL 0.5 eliminates losses for RK2 and reduces them for Euler (Gassmoeller et al., 2016).

Adaptive mesh refinement and repartitioning are handled by treating particles as cell-attached data. On local refinement, parent-cell particles are redistributed to children; on coarsening, child particles are merged to the parent; and on repartitioning, serialized particle data piggyback on the mesh framework’s cell-transfer mechanism. Because particle counts per cell may vary strongly, load balancing by cell count alone can fail badly. The implemented compromise assigns each cell the composite cost

$$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$0

with $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$1 reported as adequate for realistic runs (Gassmoeller et al., 2016).

The scaling data show that this is not a merely conceptual framework. On a Cray XC-40, a uniform-mesh benchmark with $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$2 degrees of freedom and $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$3 particles showed almost perfect strong scaling from 12 to 12,288 processes, while weak scaling with about $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$4 DoFs and $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$5 particles per process remained essentially constant from 6 to 6,144 processes. On adaptive meshes, strong scaling was again nearly linear from 12 to 1,536 cores on a case with about $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$6 DoFs and about 15 million particles. The most practically important adaptive-mesh result is that combined cell-plus-particle weighting reduced weak-scaling runtime growth from about a factor of 60 to about a factor of 3 (Gassmoeller et al., 2016).

5. Coupling, projection, and weighted deposition

Parcel-in-cell simulations support both one-way and two-way coupling. In the passive-tracer case, particles sample velocity or other fields but do not affect the Eulerian solve. In two-way coupling, particle properties are reconstructed back to the mesh and influence constitutive behavior or source terms. The geodynamics examples in ASPECT use arithmetic averaging of particle composition within a cell to influence material density and therefore the Stokes solve; more generally, the framework allows nearest-neighbor, arithmetic/geometric/harmonic averaging, distance-weighted averaging, shape-function-weighted local reconstruction, or least-squares projection onto a finite-dimensional space (Gassmoeller et al., 2016).

The finite-element work is deliberately nonprescriptive about the optimal reconstruction operator, but it emphasizes locality. Because particles are stored grouped by cell, any projection method using only particles in a cell or immediate neighbors can be implemented with optimal complexity and local communication only. If projection requires particles beyond one layer of neighboring cells, scalability becomes difficult. This is a foundational design constraint for large-scale parcel-in-cell codes on distributed-memory systems (Gassmoeller et al., 2016).

In two-phase CFD, the corresponding coupling is through cellwise momentum exchange. The gas momentum equation is

$$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$7

where $$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$8 is the source due to particles. The force model includes drag, gravity, and buoyancy; turbulent dispersion is modeled by a discrete random walk, with interaction time

$$\frac{dx}{dt}=u, \qquad \frac{du}{dt}=\frac{U^*-u}{\tau}+\left(1-\frac{\rho}{\rho_s}\right)g,$$9

and the instantaneous fluid velocity seen by the parcel is $U^*=U+U'$0. The implementation computes the resultant force due to the particle on the fluid phase while accounting for particle residence time in each crossed cell (Marzouk et al., 2024).

In weighted-parcel practice, the source-term scaling is by represented multiplicity. The same motion law is solved for a parcel as for an individual tracked particle, but the deposited momentum exchange is multiplied by the number of physical particles represented. The swirling-flow study makes clear that this is a sampling or compression strategy rather than a new force model: parcels preserve the same local physics while reducing the number of tracked Lagrangian trajectories (Marzouk et al., 2024).

6. Representative application domains

In mantle convection and related continuum problems, parcel-in-cell methods are used to carry composition-like or history-like state on adaptive unstructured finite-element meshes. The 3D entrainment benchmark reported that increasing from 40 to 500 particles per cell only slightly changed a bulk entrainment metric but substantially improved resolution of fine-scale interface structures. This indicates that modest particles-per-cell can be sufficient for averaged quantities, while sharply deformed interfaces and history fields benefit from much higher parcel density (Gassmoeller et al., 2016).

In dispersed two-phase CFD, the vertical swirling-pipe case provides a concrete benchmark for parcel weighting. The gas phase was solved on a $U^*=U+U'$1 3D wedge with 40,080 cells using PISO, backward Euler time integration, upwind convection, linear interpolation to face centers, linear diffusion terms, and a variable time step chosen to keep $U^*=U+U'$2. Three $U^*=U+U'$3-$U^*=U+U'$4 variants were compared: standard, RNG, and realizable. The standard model achieved the best overall performance; the realizable model was unable to satisfactorily predict the radial velocity and was the most computationally expensive; the RNG model predicted additional recirculation zones. For 0.5 s simulated time, the CPU times were 17.13 hr for the standard model, 19.44 hr for RNG, and 24.81 hr for realizable (Marzouk et al., 2024).

The parcel comparison in the same study is especially instructive. For 0.5 s of simulation, $U^*=U+U'$5 produced 12,099 injected parcels and 12,017 parcels remaining at the end, with a CPU time of 31.15 hr. At $U^*=U+U'$6, only 143 parcels were injected and 142 remained, but CPU time decreased only to 29.84 hr, a reduction of about $U^*=U+U'$7. Mean gas-phase velocity profiles at the sampled locations showed no discernible effect of $U^*=U+U'$8, but parcel snapshots indicated degraded particle-phase structure at $U^*=U+U'$9, including shortened convection of heavier parcels and weakened entrainment of smaller parcels outside the central bubble (Marzouk et al., 2024).

In plasma simulation, the same general coupling idea extends to electrostatic and electromagnetic PIC. The Frankfurt code developed for F8SR and FRANZ uses an electrostatic PIC formulation with first-order weighting, finite-difference Poisson solves, irregular embedded electrostatic boundaries, multiple coordinate systems, and an added Monte Carlo Collision module for simple helium plasma in a volume ion source. Jean-Luc Vay’s review covers the corresponding electromagnetic regime, including the Yee FDTD solver, Boris, Vay, and Higuera–Cary pushers, charge-conserving deposition, PSATD, moving windows, boosted frames, numerical Cherenkov instability, and mesh refinement. These are not parcel methods in the CFD sense, but they occupy the same mesh-particle methodological lineage (Joshi, 2016, Vay, 2020).

A further specialized branch is the narrowband envelope-tracked electromagnetic finite-element PIC method, ET-FEMPIC. It is designed for high-frequency, narrowband devices with self-consistent particle or plasma interaction, and combines compatible finite-element spaces, implicit Gauss-law-preserving time stepping, and envelope tracking of the electromagnetic carrier. In the beam-loaded klystron example, ET-FEMPIC reproduced the loaded cavity quality factor while using approximately one tenth the number of timesteps of regular EM-FEMPIC, and the discrete Gauss-law error remained at machine tolerance (Ramachandran et al., 2022).

7. Relation to classical PIC and principal limitations

A recurring source of confusion is terminology. In plasma physics, “PIC” often refers specifically to electrostatic or electromagnetic solvers advancing Maxwell’s equations or Poisson’s equation together with charged macroparticles. In the broader computational-science usage, however, PIC includes any method coupling discrete Lagrangian carriers to Eulerian mesh-based continuum methods. Parcel-in-cell simulations belong to this broader category, especially when parcels carry arbitrary state, are updated from Eulerian fields, and are projected back to influence continuum constitutive laws or source terms (Gassmoeller et al., 2016, Vay, 2020).

Several limitations are intrinsic to the current state of the art. On adaptive finite-element meshes, prescribed-location seeding can be much more expensive than cell-wise random generation; rejection sampling in distorted or curved cells can be inefficient; and particle loss under unusual transitions is accepted as negligible rather than eliminated. The communication strategy also leans on CFL $\tau$0, so applications with very large parcel jumps, ballistic particles, or nonlocal transport would need more general routing and search. Likewise, the framework handles compact-support particle-to-mesh coupling well, but does not address broad-kernel interactions or long-range neighbor searches (Gassmoeller et al., 2016).

In CFD parcel methods, computational savings are not guaranteed merely by reducing the number of tracked entities. The swirling-flow study reduced the tracked parcel count by nearly two orders of magnitude between $\tau$1 and $\tau$2, yet total CPU time changed only modestly because the case had low particle loading and the Eulerian gas-phase solve dominated the cost. The same study also notes a practical modeling restriction: the parcel must remain smaller than the cell under the point-force and point-mass treatment, so particles-per-parcel cannot be increased arbitrarily without compromising source localization (Marzouk et al., 2024).

Specialized PIC extensions impose additional assumptions. The Frankfurt electrostatic code does not specify the exact discrete particle pusher and does not include self-consistent magnetic-field evolution. ET-FEMPIC is not a universal replacement for time-domain EM-PIC; it assumes a narrowband electromagnetic system around a known carrier, and the presented particle pusher is nonrelativistic. Vay’s plasma-accelerator review similarly stresses that convergence in one numerical parameter is insufficient: longitudinal and transverse mesh size, time step, particles per cell, and, where relevant, boosted-frame choice and filtering strength all require systematic assessment (Joshi, 2016, Ramachandran et al., 2022, Vay, 2020).

Taken together, these results suggest that parcel-in-cell simulation is best viewed as a general computational pattern for mesh-particle coupling whose concrete realization depends strongly on the physics, the intended feedback path, and the dominant cost center. In some regimes the central challenge is scalable storage, relocation, and repartitioning on adaptive meshes; in others it is conservative weighted deposition, turbulence–particle interaction, or charge-preserving electromagnetic field advance. The common lesson is that the fidelity of parcel-in-cell simulations depends at least as much on coupling design, boundary treatment, and load balancing as on the particle or parcel ODE itself.