4. Parallel I/O fundamentals for crash codes

Crash simulations running on thousands of cores generate time-series output at a rate that quickly overwhelms naive file writing. A single explicit finite element model may produce 200–800 state variables per node per frame, and a 20-million-element model writing 500 frames can easily exceed tens of terabytes. When every rank opens the same file and writes its own portion, the metadata server becomes a bottleneck, lock contention serializes writes, and the job stalls in write() or MPI_File_write_at_all for minutes per frame. Parallel I/O middleware exists to solve exactly this problem: it aggregates thousands of small, non-contiguous rank-local writes into large, contiguous, collective transfers that the parallel file system (Lustre, GPFS, or DAOS) can absorb efficiently.

The Three-Layer I/O Stack

Understanding parallel I/O for crash codes requires separating three concerns that are often conflated:

• Transport layer (MPI-IO): Defines how a set of MPI ranks collectively read or write a shared file. It provides MPI_File_open, MPI_File_set_view, and collective operations such as MPI_File_write_at_all. MPI-IO handles the mapping from each rank's local buffer to a global file offset, and it is the layer that actually talks to the file system.
• Data model layer (HDF5, ADIOS2): Adds self-describing structure — datasets, groups, attributes, dimensions, and typed variables — on top of the byte-oriented file. This is what allows a restart file to be read back without external metadata, and what enables selective reads of a single field or a single time step.
• Access pattern layer: The application-level decision of how ranks partition the global mesh and how that partition maps to file offsets. This is where most performance is won or lost.

MPI-IO is the foundation; HDF5 and ADIOS2 are higher-level libraries that use MPI-IO (or their own transport) underneath. Choosing between them is not a matter of one being universally faster, but of matching the library's strengths to the code's output pattern.

How MPI-IO Aggregates Writes

The core mechanism is the file view. Each rank declares, via MPI_File_set_view, the region of the global file it owns. A typical crash-code pattern is a 1D block decomposition: rank r of P owns nodes [r*N/P, (r+1)*N/P), and each node contributes V contiguous bytes (the state vector). The rank's view is therefore a strided pattern — a contiguous block of (N/P)*V bytes at offset r*(N/P)*V.

When all ranks call MPI_File_write_at_all collectively, the MPI-IO implementation does not simply issue one write per rank. Instead, an aggregation strategy kicks in: a subset of ranks (the aggregators) are designated to perform the actual file system writes, and non-aggregator ranks send their data to aggregators via MPI messages. The aggregators then issue large, aligned, contiguous writes. On Lustre, this aligns writes to stripe boundaries and avoids the small-write penalty that plagues independent I/O.

HDF5: Self-Describing Time Series

HDF5 is the most common choice for crash-code restart and results files because it provides a hierarchical, self-describing container. A typical layout for a time-series crash result is:

/results
  /frame_0000
    /displacement   dataset [n_nodes, 3]  float64
    /stress         dataset [n_elements, 6] float64
    /energy         dataset [n_global]    float64
  /frame_0001
    ...
  /mesh
    /nodes          dataset [n_nodes, 3]  float64
    /connectivity   dataset [n_elements, 8] int32
  /metadata
    attributes: time_step, solver_version, material_model

Writing this with HDF5 in parallel requires the parallel-enabled build (--enable-parallel) and the use of H5Pset_fapl_mpio to attach an MPI communicator to the file access property list. Each rank then writes its hyperslab of each dataset. The critical performance knob is the collective metadata and collective data transfer mode:

hid_t dxpl = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(dxpl, H5FD_MPIO_COLLECTIVE);
H5Dwrite(dset, H5T_NATIVE_DOUBLE, mem_space, file_space, dxpl, buf);

With H5FD_MPIO_COLLECTIVE, HDF5 routes the write through MPI-IO's collective path, enabling the aggregation described above. With H5FD_MPIO_INDEPENDENT, each rank writes on its own — simpler, but typically 3–10× slower on Lustre for large frames.

HDF5's weakness for crash codes is metadata overhead at scale. Every dataset creation, every attribute write, and every H5Dclose touches the file's metadata structures. When thousands of ranks each create datasets, the metadata operations serialize. The standard mitigation is to have a single rank (or a small subset) create all datasets and groups, then use H5Dwrite collectively for the bulk data. Another is to use one file per frame rather than one file for the entire time series, trading file count for reduced metadata contention.

ADIOS2: Staging and Flexible Transport

ADIOS2 was designed for the extreme-scale case where even collective MPI-IO is not enough. Its distinguishing feature is the separation of the API from the engine. The same bpWriter code can write to a POSIX file, an HDF5 file, or a staged in-memory buffer, by changing only the engine name in the IO object:

adios2::IO io = adios.DeclareIO("crash_output");
io.SetEngine("BP5");          // or "HDF5", "SST", "POSIX"
io.SetParameter("Threads", "4");
auto var = io.DefineVariable<double>("displacement",
    {n_global, 3}, {offset, 0}, {n_local, 3});
engine.Write(var, local_buf);

The BP5 engine is ADIOS2's native format, optimized for many-rank writes with minimal metadata. The SST (Sustainable Staging Transport) engine writes to memory and allows a separate consumer process to read the data while the simulation is still running — useful for in-situ visualization of crash results without stopping the solver. For checkpoint/restart, BP5 to a parallel file system is the usual choice.

ADIOS2's DefineVariable takes the global shape, the local offset, and the local shape, which makes the rank-to-file mapping explicit and avoids the hyperslab selection boilerplate of HDF5. The trade-off is that ADIOS2's ecosystem is younger: fewer post-processing tools read BP files natively, and teams often convert to HDF5 or VTK after the run.

Choosing and Tuning

Regardless of library, three tuning levers dominate. First, aggregation: ensure collective writes are actually collective, and tune the number of aggregators (on Lustre, often one aggregator per OST or per stripe). Second, stripe alignment: set the Lustre stripe count and stripe size so that each frame's write spans an integer number of stripes; misaligned writes cause read-modify-write overhead. Third, frame granularity: writing one file per frame reduces metadata contention but increases file count; writing one file for all frames reduces file count but concentrates metadata operations. The right choice depends on frame size and file system limits.

The practical takeaway is that parallel I/O for crash codes is a layered problem, and the layer where you spend effort should match the symptom. If writes are slow because each rank writes a tiny non-contiguous chunk, fix the access pattern and use collective MPI-IO. If writes are slow because metadata operations serialize, reduce dataset creation frequency or switch to a format with lighter metadata. If the file system itself is the bottleneck, tune striping or move to a staged transport. The libraries are tools; the engineering judgment is in diagnosing which layer is actually limiting throughput.