Halo Exchange and the Surface-to-Volume Tax on Domain Decomposition

A stencil update reads a point and its neighbours and writes a new value at the same point. On one processor the neighbour reads are ordinary loads that hit cache. On a cluster, the neighbours at the edge of a rank’s subdomain belong to another rank, so each iteration must fetch a thin layer of data from the neighbours before the update can proceed. That layer is the halo or ghost layer, and the cost of maintaining it is the single most important term in the scalability of a structured-grid code.

The cost is geometric. The amount of arithmetic a rank performs grows with the volume of its subdomain, while the amount of data it must exchange grows with the surface area. Splitting a fixed problem more finely increases the surface faster than the volume, so communication per unit of work rises with the process count. That is the surface-to-volume tax, and no amount of network bandwidth removes it.

The geometry of a ghost layer

Consider a three-dimensional grid of N points per side and a seven-point stencil: a point depends on itself and its six axis neighbours. Decompose the grid into P congruent blocks, with p = P^(1/3) blocks per side, so each rank owns a cube of N/p points per side.

The rank’s interior is (N/p)^3 points. The data it needs from outside is one layer on each of its six faces, or 6 (N/p)^2 points. The ratio of exchanged points to owned points is therefore 6p/N, or 6 P^(1/3) / N. Everything about halo cost follows from that expression: it grows with the cube root of the process count and falls with the problem size. Doubling the number of ranks increases the relative communication by about 26 percent, not by a factor of two, which is why three-dimensional decompositions scale better than intuition suggests.

In bytes, with eight bytes per double, a rank exchanges 48 (N/p)^2 bytes per iteration and performs roughly seven floating-point operations per interior point, so 7 (N/p)^3 flops. The ratio of bytes to flops is 48 / (7 N/p), which improves linearly as the subdomain grows. At a subdomain edge of 100 points the ratio is below 0.1 byte per flop, comfortably compute-bound on a node whose memory system delivers a few tenths of a byte per flop. At an edge of 10 points it is close to 0.7, and the exchange is comparable to the arithmetic.

Slabs, pencils, and blocks

The decomposition shape changes the surface term without changing the volume term.

A slab decomposition gives each rank a full N x N cross-section and N/P planes. Its interior is N^3/P, and it exchanges two faces of N^2 points each, independent of P. The ratio is 2P/N rather than 6P^(1/3)/N. At P = 1000 and N = 1000, that is 2.0 exchanged points per owned point against roughly 0.06 for the block decomposition, a difference of more than an order of magnitude. The slab is simple, its messages are large and contiguous, and it is a poor choice beyond a few dozen ranks.

A pencil decomposition sits between the two: each rank owns an N x N/p column of the grid and exchanges with four neighbours. It is the usual compromise when the code’s data layout makes a full block transpose expensive, and its ratio is 4P^(1/2)/N.

The block decomposition has the best asymptotics and the worst message structure. A cube has six neighbours, but a face is a two-dimensional slice with a stride, so the data is not contiguous in memory. Sending it requires either a pack into a contiguous buffer, which costs a copy of the face, or a derived datatype that describes the stride, which lets the implementation gather directly from the array. The derived datatype avoids the copy in principle; whether it does in practice depends on whether the implementation can hand the stride pattern to the adapter or falls back to packing internally.

Message count against message size

Six faces mean six messages per iteration in each direction, twelve sends and receives if every exchange is counted separately. At a thousand ranks the faces are small: an N = 1000 grid with p = 10 gives faces of 100 by 100 points, or 80 KB per face. That is large enough that bandwidth dominates the per-message latency, but only just.

The alternative is to exchange edges and corners as separate messages so that each face message is purely interior and no data is sent twice. A seven-point stencil does not need edges or corners, so the face messages can overlap at their boundaries without harm: the duplicated points are simply overwritten with identical values. A wider stencil, such as a 27-point operator or a higher-order compact scheme, does need the edge and corner data, and then the message count rises to 26 per iteration. The extra messages are tiny, and at that point latency, not bandwidth, sets the cost.

This is the trade the decomposition shape controls. Blocks minimize bytes and maximize message count. Slabs minimize message count and maximize bytes. The right choice depends on the machine’s latency-to-bandwidth ratio and on the process count, which is why production codes often switch decomposition at a threshold rather than fixing one for all job sizes.

#include <math.h>

/* Bytes exchanged per iteration by one rank of a block decomposition,
   relative to the flops it performs for a 7-point stencil. */
double block_comm_per_flop(long n, long p, int halo)
{
    double edge = (double)n / cbrt((double)p);
    double face_bytes = 6.0 * edge * edge * halo * sizeof(double);
    double flops = 7.0 * edge * edge * edge;
    return face_bytes / flops;      /* falls as 1/edge */
}

/* The same quantity for a slab decomposition: two full faces per rank. */
double slab_comm_per_flop(long n, long p, int halo)
{
    double face_bytes = 2.0 * (double)n * (double)n * halo * sizeof(double);
    double flops = 7.0 * (double)n * (double)n * (double)n / (double)p;
    return face_bytes / flops;      /* grows linearly with p, unlike the block */
}

The two expressions differ in the exponent of p, and that exponent is the whole argument for three-dimensional decomposition.

Overlap and the packing cost

The halo exchange has a natural overlap structure. The interior points that do not touch a boundary can be updated while the boundary messages are in flight, and the boundary update waits for the receives. The non-blocking interface exists for this pattern: post the receives, post the sends, update the interior, then wait and update the boundary.

The fraction of work that can overlap is the interior volume divided by the total, which is (edge - 2)^3 / edge^3 for a halo of width one. At an edge of 100 that is over 94 percent; at an edge of 10 it is about 51 percent. Overlap therefore works well exactly where the communication ratio is already small, and poorly where communication matters most. The tax is not just the bytes; it is that the bytes cannot be fully hidden in the regime where they are expensive.

Packing adds a cost that the byte count does not show. Copying a strided face into a contiguous buffer touches edge^2 points and writes them again, which is memory traffic on the same order as the arithmetic for a cheap stencil. A seven-point update reads seven values and writes one, so a pack that reads and writes one value per boundary point is a noticeable addition. Derived datatypes avoid the explicit copy but may not avoid an internal one, and the only way to know which happens is to inspect the implementation or measure. Since measurement is the only honest answer, the design guidance is to keep the face layout regular enough that a derived type can describe it with a small number of extents.

Where the decomposition breaks

Four boundary conditions cause trouble.

The first is a subdomain thinner than the halo. If N/p is one or two points and the halo is one point wide, the interior vanishes and every point is boundary. The update still works, but the communication ratio is meaningless and the code spends all its time exchanging. A decomposition that assigns fewer than a few halo widths per rank in any dimension is not a decomposition.

The second is a process count that does not factor into the grid. A block decomposition wants P to be a product of three integers close to P^(1/3). A prime process count forces a slab or a very uneven block. Implementations handle this by factoring and allowing the last block in each dimension to be smaller, which introduces imbalance that grows with the remainder.

The third is a rank that is its own neighbour. With periodic boundaries and a process grid of size one in a dimension, the rank to the left and the rank to the right are the same rank, and in the extreme case both are the caller. A blocking send to self deadlocks; the standard’s self-send rules and MPI_Sendrecv handle it, but only if the code does not special-case the neighbour arithmetic incorrectly.

The fourth is load imbalance. Equal subdomains assume equal work per point, which fails for adaptive refinement, for immersed boundaries, and for variable material properties. The surface-to-volume analysis then describes the average rank and not the slowest one, and the slowest rank sets the iteration time because the exchange synchronizes neighbours. The strong scaling discussion explains why that maximum, rather than the mean, is the quantity that matters.