In-Situ Analysis and the Bandwidth That Is Not There

The traditional pipeline for simulation output is a file per step, written by the solver, read by a visualization tool in a separate process. It works because both sides are slow relative to the compute, and it fails when they are not. The storage system’s bandwidth and the solver’s rate of producing data grow at different rates, and once the solver can produce data faster than the parallel filesystem can absorb it, the cost of writing becomes a bottleneck that the simulation budget cannot cover. The response is to stop moving the data: run the analysis inside the simulation process, or in a process that receives data over the interconnect at a rate that matches production, and write only what the analysis reduces to.

The change is not merely an optimization. Coupling the analysis into the solver changes the memory footprint, the parallel decomposition, the scheduling, and the failure modes, and it makes the analysis a part of the simulation that must be reasoned about when either changes.

The arithmetic of not writing

Consider a volume of $512^3$ samples in single precision: 134 million values, 536 MB per field. A simulation with twenty material fields, saved every hundredth step over ten thousand steps, writes

$$100 \times 20 \times 536\ \mathrm{MB} \approx 1.07\ \mathrm{TB},$$

excluding halos, metadata, and any duplicate copies kept for restart. At a sustained parallel write rate of a few hundred MB/s per storage target, and with the last steps of a campaign often needing to be saved more frequently than the first, the write time becomes comparable to the compute time rather than a fraction of it.

The data that the analysis actually needs is much smaller. A volume rendering at $1024 \times 1024$ pixels is 4 MiB per frame in RGBA8, and an isosurface at this resolution is tens of megabytes if written as a mesh. Where the study’s scientific question really is about the images, the reduction from terabytes to megabytes is the whole argument. Where the question is about the field itself, the reduction is lossy in a way that must be designed rather than accepted.

The trap is that the reduction is chosen before the question is known. A pipeline that writes only renderings cannot answer a later question about a feature the images do not show, so the output is a fixed vocabulary of images rather than a dataset. This is the central design tension of in-situ work: the less data is kept, the more the analysis must anticipate.

Coupling modes and their synchronisation costs

In the weakest coupling, the analysis library is linked into the solver and called between timesteps. The solver hands over its arrays, typically through a description of the decomposition rather than by copying, and the analysis operates on the data in place. There is no data movement over the network, but there is a synchronization: the analysis runs while the solver waits, so its cost is added to the wall clock of the simulation, and every core not participating in the analysis is idle during that period.

The cost is therefore not zero even when the data does not move. A rendering pass over a field distributed across thousands of ranks requires each rank to contribute its portion of the image, which requires a compositing step; compositing is a collective operation and its cost grows with the number of participants. If the analysis runs on the solver’s ranks, the compositing cost is paid on the same cores that could be advancing the simulation, and its all-to-all character leaves most of them idle while it completes.

In-transit computing decouples the two: data is sent over the interconnect to a set of dedicated analysis ranks, which may be fewer than the solver ranks, so the analysis can be performed at a rate the solver can sustain while the simulation continues. The cost model changes from idle cores to network bytes. A field of 536 MB per step over a 100 Gb/s link takes roughly 43 ms at line rate before protocol overhead, which sets a ceiling on the sustainable output rate regardless of how fast the analysis is. The analysis ranks must also be able to consume a step before the next arrives, so the pipeline’s throughput is the maximum of the two stages and its latency is the sum.

A middle mode runs the analysis on a subset of the solver ranks and redistributes the data among them, which reduces the idle fraction without dedicating whole nodes. The redistribution is itself a communication pattern, and where the solver’s decomposition is a regular block grid and the analysis needs a different partitioning, the cost is a full data exchange.

Domain decomposition cuts the field where the analysis needs it whole

Every parallel simulation partitions its domain, and every analysis that looks across cell boundaries needs data that a single rank does not own. An isosurface extraction on a block boundary needs the neighbour’s samples to determine which cells are crossed. A gradient stencil at the edge needs values one layer out. A streamline that leaves the rank’s block needs the field on the rank next door.

The standard answer is a ghost layer: each rank stores a copy of its neighbour’s boundary values, exchanged before the analysis runs. One layer suffices for a first-order stencil and a single-cell-wide isosurface crossing; a streamline requires continuous tracking across ranks, which means either a distributed data structure that follows the curve or an exchange protocol that migrates the integration state to the owning rank when it crosses. The latter is the usual implementation and it turns the integration into a message-passing algorithm with a latency per crossing rather than a local computation.

Ghost exchange interacts with the timestep boundary. The solver typically exchanges ghosts only when it needs them for its own stencil, and the analysis may need a wider layer or a fresher copy. Requiring an extra exchange per analysis step adds a collective to every timestep, and its cost is independent of the analysis’s own cost, which makes it dominant for cheap analyses.

The decomposition also determines what cannot be computed locally. Global reductions, a volume integral, a maximum over the whole domain, a histogram over all voxels, are collectives, and their cost grows with the logarithm of the rank count and with the size of the reduced quantity. A histogram is the expensive case because the reduced object is $O(\text{bins})$ rather than scalar: with 256 bins, a reduction over a thousand ranks moves 256 KB per step, and a 2D histogram with $256^2$ bins moves 256 MB before any compression. Per-rank pre-aggregation reduces this, since each rank’s local histogram can be sparse, but the sparse pattern differs per rank and the merge is irregular.

/* Ghost exchange before analysis: one layer, periodic in x.
   Plane layout: index 0 is the left ghost, 1..nx are the compute planes,
   nx+1 is the right ghost. */
static void exchange_x(const Grid *g)
{
    size_t plane = (size_t)g->ny * g->nz;
    float *d = g->data;

    /* own first compute plane goes left; the right neighbour's first plane
       arrives into the right ghost */
    MPI_Sendrecv(d + plane, plane, MPI_FLOAT, g->left, 0,
                 d + (size_t)(g->nx + 1) * plane, plane, MPI_FLOAT,
                 g->right, 0, g->comm, MPI_STATUS_IGNORE);

    /* own last compute plane goes right; the left neighbour's last plane
       arrives into the left ghost */
    MPI_Sendrecv(d + (size_t)g->nx * plane, plane, MPI_FLOAT, g->right, 1,
                 d, plane, MPI_FLOAT,
                 g->left, 1, g->comm, MPI_STATUS_IGNORE);
}

The ghost planes must be part of the same allocation as the compute planes, so that an analysis stencil can index i - 1 and i + 1 without a branch at the block boundary. Getting the send and receive directions transposed is the classic source of a seam at every block boundary in an extracted isosurface, because the extraction then reads a ghost value that was filled with the wrong neighbour’s data.

Scheduling, checkpointing, and the failure interaction

A solver that checkpoints periodically has an I/O pattern that in-situ analysis changes. If the analysis writes images at every step and the solver checkpoints occasionally, the two share the filesystem, and a burst of analysis output during a checkpoint window can extend the checkpoint’s duration. Scheduling analysis to avoid the checkpoint window is the obvious mitigation and it is not free: it means the analysis runs at a different step rate than the simulation, which for time-varying features means some steps are never looked at.

Failure modes are more consequential. A crash in the analysis library now crashes the simulation, and a job that has run for days loses its progress unless a checkpoint exists immediately before. The mitigation used in practice is to run the analysis in a separate process, even on the same node, so that a fault is contained; that reintroduces a copy of the data through shared memory, which costs bandwidth but not as much as the network, and it keeps the solver’s address space intact.

Restart introduces the inverse problem. A simulation restarted from a checkpoint has an analysis state, accumulator arrays, seeded positions, temporal filters, that must also be restored, and that state is not part of the solver’s checkpoint unless the analysis is integrated into the checkpoint format. An analysis that maintains a running average over timesteps and silently restarts from zero produces an image that represents the second run only, with a discontinuity at the restart point that is invisible in a single frame and obvious in a sequence.

Determinism is a related constraint. A composite of per-rank images depends on the order in which contributions are accumulated when the arithmetic is not associative in floating point, and the order is determined by the communication pattern, which can vary with the number of ranks. Running the same simulation on a different node count then produces an image that differs in the last bits, which is harmless for a picture and unacceptable when the analysis output is itself a scientific result. Fixed reduction trees, or reductions in double precision with a specified order, are the standard repairs.

What survives from the traditional pipeline

In-situ analysis does not replace the offline tools; it produces data for them. The useful patterns keep the reduction honest by emitting data that can be re-examined: a set of renderings from a fixed camera and transfer function, which is what the Cinema approach formalizes as a database of images that a viewer can query; a downsampled field at a resolution that fits in storage; or the extracted geometry rather than the volume.

Each of these sits at a different point on the same axis. A downsampled field preserves the possibility of a different transfer function later, at a resolution that has already discarded detail. Extracted geometry preserves exact positions of one level set and nothing about the rest of the field. Images preserve the appearance under one set of rendering parameters. Choosing among them is choosing which future questions the data will be able to answer.

The rendering itself is constrained by the same mechanisms as any other volume rendering: opacity correction and transfer function design in Direct Volume Rendering and the Opacity Correction Nobody Applies, and multi-dimensional classification in Multi-Dimensional Transfer Functions and the Gradient-Magnitude Histogram. The step that turns per-rank images into one image is compositing, and its cost at scale is the subject of Sort-Last Compositing and the Depth Order Across Ranks.