| title | lptlib: A parallel Lagrangian particle tracking library for compressible CFD data and tracer-response analysis in optical velocimetry | ||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| tags |
|
||||||||||||||||
| authors |
|
||||||||||||||||
| affiliations |
|
||||||||||||||||
| date | 24 August 2026 | ||||||||||||||||
| bibliography | paper.bib |
Laser-based velocity measurement never observes the fluid directly. Particle image velocimetry (PIV) seeds a flow with microscopic tracer particles, photographs them twice in quick succession, and infers the velocity field from how far they traveled in between [@raffel2018piv]. The inference assumes that the particles go where the fluid goes, and in high-speed flow that assumption fails. A tracer crossing a shock wave carries too much inertia to decelerate as abruptly as the gas around it, so the measured velocity jump is smeared and the measured gradients are weaker than the real ones [@melling1997tracer]. The size of that bias depends on the particle, not on the instrument.
lptlib quantifies that bias. Given a computed flow field, the library releases particles of a specified size and density into it and integrates their trajectories, so the motion a camera would have recorded can be compared against the motion of the fluid itself. The difference between the two is the measurement error attributable to the tracer, shown in \autoref{fig:lag} for three tracer sizes.
Working from a simulated flow rather than an experiment makes the comparison possible, because the true fluid velocity is known everywhere. The library reads structured, multi-block, curvilinear grids and flow solutions in the PLOT3D format [@walatka1990plot3d], interpolates the flow state to arbitrary points inside the curvilinear domain, and advances both fluid streamlines and inertial particle paths. Large ensembles are tracked in parallel through multiprocessing, thread pools, or MPI [@dalcin2021mpi4py], and the scattered tracks are reduced back onto a user-defined Eulerian mesh to produce PIV-like fields. The tracker itself is not specific to velocimetry, and advances inertial particles through any resolved compressible field.
Correcting inertial bias requires a drag law that stays valid across every regime a tracer visits while crossing a shock. A 1 micron alumina particle in a Mach 5 flow passes from continuum to slip to near-free-molecular conditions within a few particle diameters, and no single classical closure covers that path. The practical questions the community asks are consequently comparative. How much does the answer change if the tracer is 300 nanometers instead of 1 micron, or if Henderson's closure [@henderson1976drag] is used rather than Loth's [@loth2008compressibility]? Answering either requires a tool where particle properties and the drag closure are both free parameters over the same flow field and the same integrator.
Such a tool has not been openly available. Researchers in this area have generally written their own single-purpose trackers, which are rarely released, rarely tested, and impossible to compare against one another. lptlib was built to close that gap. It exposes nine drag closures through a single argument, from Stokes and Oseen [@oseen1927hydrodynamik] creeping flow through the Schiller-Naumann [@schiller1933drag] and Melling [@melling1997tracer] continuum corrections and a Cunningham-type slip correction [@cunningham1910slip] to the compressible and rarefied models of Henderson, Loth, Tedeschi [@tedeschi1999motion], and Subramaniam and Balachandar [@subramaniam2022particle], with temperature-dependent viscosity from Keyes' correlation or Sutherland's law [@sutherland1893viscosity]. All nine share one integrator, and the results reduce to a gridded field that downstream synthetic-imaging tools can consume. Optical velocimetry motivated the library and is what its shipped test cases cover, but nothing in the tracker is specific to it. The drag closures, the integrator, and the Eulerian reduction apply unchanged wherever inertial particles are advected through a resolved compressible field, from particle dispersion to erosion and impingement. The intended users are therefore experimentalists interpreting high-speed PIV data, and computational researchers studying one-way coupled particle-laden compressible flow.
Open-source Lagrangian particle tracking is mature, but it is partitioned into communities whose assumptions do not transfer to compressible aerodynamics. Geophysical transport frameworks are the largest group. Parcels [@delandmeter2019parcels], OpenDrift [@dagestad2018opendrift], and MPTRAC [@hoffmann2025mptrac] advect particles through ocean or atmosphere velocity fields with great sophistication in interpolation and scale, but their particles are effectively passive, and the compressible and rarefied drag physics that dominates a tracer crossing a shock has no counterpart in them.
Experimental particle tracking velocimetry tools such as MyPTV [@shnapp2022myptv] sit on the other side of the measurement. They consume camera images and produce trajectories, whereas lptlib consumes a CFD field and produces trajectories, then reduces them back to a grid. The two are complementary.
Closest in intent is ppiclF [@zwick2019ppiclf], a parallel particle-in-cell library for particle-laden flow that scales to very large processor counts. However, ppiclF is written in Fortran and designed to be linked into a host solver, so it tracks particles during a simulation rather than after one. Solver-embedded Lagrangian cloud libraries such as those in OpenFOAM [@weller1998openfoam] share that constraint and a dependence on the solver's own mesh format. lptlib instead operates on a written solution file, so tracer response can be studied in an archived result without rerunning it, and the drag closure can be varied without touching the solver.
General-purpose visualization reaches part of this space. VTK and ParaView [@ahrens2005paraview] read PLOT3D files and offer particle tracers, but those tracers integrate massless streamlines, so they answer where the fluid goes rather than where a tracer of a given size and density goes. Beyond visualization, open PLOT3D post-processing is dominated by legacy Fortran utilities and closed commercial tools.
Three decisions shape the library.
The first is adaptive integration keyed to trajectory curvature. Particle response near a shock is stiff, and a globally small time step is wasteful everywhere else, so lptlib adapts the step to the angle between successive displacements, refining where a path deflects sharply and coarsening again in smooth regions. Deflection was chosen over a conventional local error estimate because the quantity of interest is the tracer's spatial lag through a discontinuity, and deflection responds to the discontinuity directly rather than to the smoothness of the interpolated field. A blow-up detector guards the stiff branch, reduces the step, retries, and terminates the trajectory with a diagnostic if the step cannot be taken. The step is restored once integration succeeds again, so one shock crossing does not cripple the remainder of a path.
The second is one-way coupling. Particles do not act back on the fluid, which is accurate at the low seeding densities used in PIV and which makes every trajectory independent. Independence allows the same tracking code to be driven by multiprocessing, a thread pool, or MPI without changing the numerics.
The third concerns input throughput. An ensemble reads the flow once and shares it across every particle, so the cost that matters is unsteady tracking, where a run loads a whole directory of PLOT3D solution files before a single trajectory is advanced. The PLOT3D reader loads each file in one buffered numpy.fromfile call and reconstructs multi-block arrays by strided slicing and Fortran-order reshaping, avoiding Python-level iteration over grid points [@harris2020numpy; @virtanen2020scipy]. A reproducible benchmark under benchmarks/ measures the effect, comparing readers only when they perform matched work. On a 2.04 million-point grid the strided reader beats an otherwise identical naive Python reader by a ratio of 27.8, with a 95 percent bootstrap interval of 26.1 to 28.8. Against a compiled Fortran reader doing the identical transpose it is about 1.4 times slower, while the full read_grid, which also builds the padded double-precision array and the per-block bounds, is about 1.75 times faster than Fortran doing that same work. Ratios rather than absolute times are reported, since absolute times move by more than an order of magnitude with page-cache state. benchmarks/README.md records the machine and explains how to add another.
lptlib grew out of, and now underpins, a line of work on particle response and bias correction in supersonic PIV spanning 2020 to 2026. The methodology was first set out for turbomachinery CFD data in @kalagotla2020rotor37. The library then provided the tracking framework and PLOT3D post-processing interface of @kalagotla2023oop and @kalagotla2023api, produced the synthetic PIV data of @kalagotla2024synthetic, supplied the drag-model comparison and particle-size estimation of @kalagotla2024drag, and generated the training data behind the bias corrector of @kalagotla2025bias.
That thread has since appeared in peer-reviewed form. The neural bias corrector of @kalagotla2026bicsnet was trained on 600 synthetic cases generated by this library and reduces inertial bias by up to 87 percent within its training distribution. The solver described there as in-house is lptlib, under the name it now carries.
Use now extends past the group that wrote the library. Two 2026 conference papers from Florida State University apply it to a supersonic wall-mounted hemisphere flow [@kalagotla2026hemisphere] and to uncertainty quantification of particle-inertia effects [@kalagotla2026uq], and the co-authors on those papers took no part in writing it. Adoption by a wholly independent party is not claimed, since the library's first author is common to all of this work.
The library is installable from PyPI, documented, and covered by a test suite that runs on Linux and macOS across Python 3.10 to 3.13. Its verification tests check the implemented oblique-shock relations against an independently coded theta-beta-Mach solution and against published gas-dynamics tables. The particle integrator is checked against a high-accuracy reference integration of the same equation of motion, which it reproduces to a relative velocity error of about 0.1 percent under both Stokes and standard sphere drag, and a physics-matched cross-code comparison against OpenFOAM parcel tracking is provided under benchmarks/tracking/.
Generative AI coding assistants were used during the development of lptlib, in parts of the implementation, the test suite, the documentation, and the preparation of this paper. All such output was reviewed, corrected where necessary, and validated by the authors, who take full responsibility for the software and for the claims made here. The scientific formulation, the choice of models, the design decisions, and the interpretation of all results are the authors' own.
The authors thank Harpreet Chhabra for contributions to early test cases, and Daniel Cuppoletti for discussions on the experimental context. This work was carried out in the Department of Aerospace Engineering and Engineering Mechanics at the University of Cincinnati.