Running in Parallel

Monte Carlo transport is the rare scientific calculation that parallelizes almost for free. Each history is independent of every other: a neutron is born, scatters, and is absorbed or escapes without any knowledge of the other neutrons in the run. Split a million histories across sixteen processors and you have sixteen problems of sixty-odd thousand histories that never need to speak to each other until their tallies are summed at the end.

MCNP goes further and makes the split invisible in the answer. At the start of each history the random number generator advances by a fixed stride, so history number 731,004 draws the same sequence of random numbers regardless of which worker happens to run it. A fixed-source answer therefore does not depend on how you divided the work, which is what lets you develop a deck on a laptop and produce the identical result on a cluster.

What parallelism does not do is repeal the 1/√N law from output analysis. Sixteen cores let you run sixteen times the histories in the same wall-clock time, and sixteen times the histories cuts the relative error by a factor of four. That is a real gain, and it is a much smaller gain than the core count suggests. A calculation that is a factor of ten away from useful statistics needs a hundredfold increase in work, which is a variance-reduction problem wearing the costume of a hardware problem. Get the figure of merit up first, then buy cores.

On a single machine, threading is all you need and it is one keyword. MCNP calls its threads tasks, and tasks n starts n of them sharing one address space.

mcnp6 i=pin_cell.i o=pin_cell.out tasks 8

Because the threads share memory, the cross-section tables are loaded once rather than eight times, which is the practical reason to prefer threads over processes on one machine. A deck with a large FMESH is the case where this matters most: every algorithm except the remote-memory one uses memory more efficiently as you raise the thread count and lower the number of separate processes.

Set n to the number of physical cores, not the number of logical processors the operating system reports. MCNP transport is limited by memory bandwidth and by branch-heavy table lookups, so the second hardware thread on a core has little of its own to do and mostly contends with the first. On a machine advertising sixteen logical processors from eight physical cores, tasks 8 will usually beat tasks 16, and it will never lose by much.

Threads cannot cross a machine boundary. For that MCNP uses MPI, which needs a separate MPI installation and a different executable — and the name of that executable depends on your platform and MPI implementation. On Linux it is mcnp6.ompi for OpenMPI and mcnp6.mpich for MPICH; on macOS, mcnp6.mpi; on Windows, mcnp6.mpi.exe for Microsoft MPI.

The detail that costs people a core, and occasionally an entire node, is that MCNP's MPI mode is manager–worker. The count you give mpiexec is the total number of processes including the manager, and the manager tracks no particles. Ask for 16 and you get 15 workers.

# 17 processes: 1 manager + 16 workers actually tracking particles
mpiexec -n 17 mcnp6.mpi i=pin_cell.i o=pin_cell.out

# The same executable runs serially with no mpiexec at all — useful for
# checking a deck for input errors, or for plotting geometry
mcnp6.mpi i=pin_cell.i ip

The MPI executables include OpenMP as well, so the two models compose. Each worker can run its own threads, and the number of processes tracking particles is (m − 1) × n — one fewer than the process count, times the threads per process.

# 4 nodes: 1 manager + 3 workers, 6 threads each = 18 tracking particles.
# Not 24. The arithmetic here is (m - 1) x n, and forgetting the manager is
# the most common way to under-fill a cluster allocation.
mpiexec -n 4 mcnp6.mpi i=core.i o=core.out tasks 6

The usual arrangement on a cluster follows from the memory argument above: one MPI process per node, with threads equal to the cores on that node, plus one extra process for the manager. Four compute nodes of 24 cores is therefore -n 5 ... tasks 24, and the node hosting the manager will be nearly idle. On a large allocation that waste is negligible; on two nodes it is a quarter of your machine, which is an argument for threading alone until you genuinely need a second node.

For problems where workers finish at noticeably different rates, the BALANCE keyword turns on load balancing between MPI processes.

Threading required extensive changes to the MCNP source, and a number of features were never converted. If a deck uses one of them, MCNP prints a warning and runs the whole problem on a single thread — it does not refuse, and it does not stop. A run that is mysteriously eight times slower than the last one is almost always this, so the warning is worth reading rather than scrolling past.

The ones most likely to appear in reactor and shielding work are TMESH tallies, the legacy non-HDF5 PTRAC, surface source writing with SSR, delayed neutrons and photons through ACT's CINDER path, and the FMULT fission multiplicity models — CGMF, FREYA, and the LLNL models. DBCN event logging and HISTP file generation do it too. So does transporting any particle other than neutrons, photons, and electrons, and so does relying on model physics where a cross section is missing.

Two of these have modern alternatives that do thread. Prefer FMESH to TMESH, which is the current recommendation for other reasons as well, and prefer the HDF5 PTRAC format to the legacy one.

Run the same deck at several thread counts and compare. Keep nps fixed so that each run does identical work, and keep the runs long enough that startup and cross-section loading are not a measurable fraction of the total.

for n in 1 2 4 8 16; do
  /usr/bin/time -f "tasks=$n  wall=%e s" \
    mcnp6 i=scaling_test.i o=out_$n tasks $n
done

# Speedup S(n)   = T(1) / T(n)
# Efficiency E(n) = S(n) / n     -- 1.0 is perfect, 0.8 is a good result

There is a better measurement available, and MCNP has already computed it for you. The figure of merit is 1/(R²T), so it already contains the run time — at fixed statistics it is directly proportional to speedup, and it comes from inside MCNP rather than from a shell timer that also counted your file system. Compare the FOM printed by each run. If the FOM at eight threads is not close to eight times the FOM at one, the threads are not doing what you think, and the usual culprit is one of the features above.

Expect efficiency to hold well through a single socket and then soften. Beyond that the limits are hardware rather than MCNP: memory bandwidth shared across cores, and cross-section tables that stop fitting in cache. Deep variance reduction can also flatten scaling, because splitting turns one history into a tree of correlated work that lands on one thread.