The Monte Carlo Transport Project at Lawrence Livermore National Laboratory develops two Monte Carlo transport codes used in production by a sizable internal user community. Mercury is a Monte Carlo particle transport code used to model the interaction of neutrons, gammas, and light ions with a material. Imp is an implicit Monte Carlo thermal x-ray photon transport code used to model the interaction of x-ray photons with a material. This paper describes the two codes and highlights recent developments.

Mathematics and computational methods are of fundamental importance in the design and analysis of nuclear systems. The “Math & Comp” research community is actively engaged in advancing numerical algorithms and computational methods to simulate nuclear systems with increasingly higher fidelity. The American Nuclear Society (ANS) 2021 International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2021) continued the rich tradition of this biennial conference series, showcasing preeminent computational science research for nuclear applications. M&C 2021 was an M&C conference like none before. Much of the planning for the conference was conducted amid the uncertainty of the COVID-19 pandemic via e-mail and Zoom sessions in 2020 and 2021. The conference was originally scheduled for April 11–15, 2021, in Raleigh, North Carolina. The pandemic resulted in the delay of the conference to October 3–7, 2021, and the COVID-19 Delta surge led to the last-minute decision to go fully virtual for the first time in the history of the conference. A constant focus in planning the conference was a strong commitment to excellence in the technical program, a hallmark of this series of conferences. The papers in this special issue of Nuclear Science and Engineering represent selected highlights of the research presented at M&C 2021. The themes of M&C 2021 were mathematical and computational methods, numerical analysis, computer codes, and high-performance computer architectures for solving problems in nuclear science and engineering. The M&C 2021 technical program was based on 14 technical tracks:

This paper considers Monte Carlo photon transport applications on heterogenous compute architectures with both CPUs and GPUs. Previous work on this problem has considered only meshes that can fully fit within the memory of a GPU, which is a significant limitation: many important problems require meshes that exceed memory size. We address this gap by introducing a new dynamic replication algorithm that adapts assignments based on the computational ability of a resource. We then demonstrate our algorithm’s efficacy on a variety of workloads, and find that incorporating the CPUs provides speedups of up to 20% over the GPUs alone. Further, these speedups are well beyond the FLOPS contribution from the CPUs, which provide further justification for continuing to include CPUs even when powerful GPUs are available. In all, the contribution of this work is an algorithm that can be applied in real-world settings to make more efficient use of heterogeneous architectures.

A graphics processing unit (GPU) has become a core technology for modern supercomputers. Applications that once ran on supercomputers are being forced to make significant changes to their designs to utilize these new machines. This paper introduces the concept of Thin-Threads as a method for history-based Monte Carlo transport applications on GPUs. The key principles behind Thin-Threads are light memory usage and communication and managing data race issues via atomics. We show that we can achieve a 10x speedup when moving from the traditional method to Thin-Threads on GPUs. Additionally, we demonstrate the viability of the Thin-Threads model at scale for GPU and CPU platforms.

The Chord Length Sampling (CLS) algorithm is a powerful Monte Carlo method that models the effects of stochastic media on particle transport by generating on-the-fly the material interfaces seen by the random walkers during their trajectories. This annealed disorder approach, which formally consists of solving the approximate Levermore-Pomraning equations for linear particle transport, enables a considerable speed-up with respect to transport in quenched disorder, where ensemble-averaging of the Boltzmann equation with respect to all possible realizations is needed. However, CLS intrinsically neglects the correlations induced by the spatial disorder, so that the accuracy of the solutions obtained by using this algorithm must be carefully verified with respect to reference solutions based on quenched disorder realizations. When the disorder is described by Markov mixing statistics, such comparisons have been attempted so far only for one-dimensional geometries, of the rod or slab type. In this work we extend these results to Markov media in two-dimensional (extruded) and three-dimensional geometries, by revisiting the classical set of benchmark configurations originally proposed by Adams, Larsen and Pomraning [1] and extended by Brantley [2]. In particular, we examine the discrepancies between CLS and reference solutions for scalar particle flux and transmission/reflection coefficients as a function of the material properties of the benchmark specifications and of the system dimensionality. (C) 2017 Elsevier Ltd. All rights reserved.

This paper presents ongoing research efforts at Lawrence Livermore National Laboratory to enable the Mercury Monte Carlo particle transport code to run efficiently on current and upcoming advanced computing architectures. We briefly describe the Quicksilver proxy application that we have developed to enable more rapid prototyping of new algorithms and to engage the external computer vendor communities. We present research and development efforts with Quicksilver and Mercury focused toward the Trinity machine at Los Alamos National Laboratory that uses both Intel Xeon Haswell processors and Intel Xeon Phi Knights Landing many integrated core processors. Finally, we describe research into Monte Carlo event-based and history-based algorithms for the Lawrence Livermore National Laboratory Sierra machine that will use IBM Power processors along with Nvidia Volta graphics processing unit architecture accelerators.

The material properties at a given spatial location are known only statistically in a stochastic medium [1]. The most common approach to solving particle transport problems involving binary stochastic media is to use the atomic mix (AM) approximation [1] in which the transport problem is solved using ensemble-averaged (homogenized) material properties. A common deterministic model developed for solving particle transport problems in binary stochastic media is the Levermore-Pomraning (LP) model [1, 2]. Adams et al. [3] investigated the accuracy of the LP model in onedimensional (1D) planar geometry by comparing LP solutions to benchmark solutions obtained by sampling large numbers of material realizations using Markovian statistics, solving a transport problem for each realization, and ensemble averaging the results. Zimmerman and Adams [4] proposed a Monte Carlo LP algorithm that produces the Levermore-Pomraning solution and another local realization preserving (LRP) algorithm that typically improves accuracy by locally preserving the material realization encountered by the Monte Carlo particle. These algorithms have been extensively studied using 1D planar geometry benchmark investigations [4, 5]. The Monte Carlo LP algorithm is generally more accurate than the AM approximation, and the LRP algorithm is generally more accurate than both the AM and LP approximations. The Monte Carlo LP and LRP algorithms have been extended for use in thermal photon implicit Monte Carlo (IMC) applications [6, 7]. Benchmark comparisons in 1D planar geometry demonstrate that the LRP algorithm is generally more accurate than the LP algorithm for these thermal photon transport problems. A significant limitation of these previous benchmark comparisons is the restriction to 1D planar geometry. Brantley and Martos [8] generated three-dimensional (3D) binary stochastic medium benchmark results for a particle transport problem with varying spherical inclusion mean chord lengths and three different spherical radius distributions. That research was an initial effort toward generating 3D benchmark results relevant to particle transport through a stochastic medium composed of inclusions of varying size and shape such as may be encountered in a turbulent medium. An example material realization from that research is shown in Fig. 1(a). (The range of material properties considered in that work was smaller than in the 1D benchmark studies.) Brantley [9] subsequently investigated the accuracy of a 3D Monte Carlo LP algorithm for these 3D benchmark problems (for the exponential spherical radius distribution case) and found that the LP approximation is significantly more accurate than the AM approximation. In addition, the LP approximation qualitatively captures trends in the solution that are not captured by the AM approximation. A limitation of these previous 3D benchmark compar(a) From Ref. [8] (b) From Ref. [10]

Like many other code teams, the developers of the Mercury Monte Carlo Transport code at Lawrence Livermore National Laboratory are being forced by the arrival of GPU-based supercomputers to substantially refactor their application to obtain acceptable performance on new architectures. This paper describes how we have designed, developed, and used Quicksilver, a proxy application for Mercury, to assist in this process. We explain the goals and principles that informed Quicksilver's design and describe the features of Mercury that are and are not represented by the proxy. We recount some experiences testing Quicksilver and changes made to improve its ability to represent Mercury. Finally we explore how we have used Quicksilver and QS lite, a proxy app for Quicksilver, as prototypes to test potential refactorings and provide performance results on GPU and CPU platforms.

High performance computing environments are progressively moving towards many-core computing architectures. The Los Alamos National Laboratory Trinity machine, available in late 2016, will use both Intel Xeon Haswell processors and Intel Xeon Phi Knights Landing many integrated core (MIC) coprocessors. The Lawrence Livermore National Laboratory Sierra machine, available in 2018, will use an IBM PowerPC architecture along with Nvidia graphics processing units (GPUs). Applications that must work in this supercomputing environment must continue to adapt in order to take advantage of the diverse hardware architectures that are coming. A significant consideration is not only the performance of the application on a given platform but also the portability of the application to other platforms. The algorithmic improvements presented in this paper build upon recently-reported work [1] on event-based Monte Carlo transport in the ALPSMC code that models particle transport in one-dimensional binary stochastic media [2]. That paper discussed the lack of available vectorization in the traditional history-based algorithm used for Monte Carlo transport and presented a data parallel event-based algorithm implemented using the Nvidia Thrust library [3] for portability. The performance of the data parallel eventbased algorithm implemented using Thrust was compared to a native CUDA [4] implementation. The conclusions from that work were that the Thrust library abstraction technique caused too significant a loss in performance but that the event-based method was a viable option that should be further investigated. In this paper, we describe algorithmic improvements to the data parallel event-based algorithm previously presented [1]. We made further algorithmic optimizations to the event-based CUDA implementation, most notably: data structure changes, a new conditional particle removal scheme in the event-based process, and the use of multiple GPUs. In addition to improvements to the algorithm, we reimplemented the Thrust version from the now further optimized CUDA version, giving a greater chance for success at a performant abstraction. Finally, we revisited our previous assumptions about the inability of the history-based method to achieve performance on vector style architectures such as the MICs and GPUs, with surprising and promising results.

We describe the challenges posed to the Mercury Monte Carlo particle transport code development team from emerging and future advanced computing architectures. We review recent work to scale Mercury to large numbers of MPI processes as well as to improve compute node parallelism via OpenMP threading and demonstrate these capabilities using a reactor eigenvalue calculation. We then describe initial progress for enabling Mercury for the Intel Xeon Phi-based MIC architecture. We present preliminary results of research investigations into the use of event-based algorithms in a Monte Carlo test code for application to GPU architectures. We then briefly describe work to enable storage of nuclear data in shared memory and to enable the use of the Generalized Nuclear Data format in Mercury via the General Interaction Data Interface.

We numerically investigate the accuracy of two Monte Carlo algorithms originally proposed by Zimmerman [1] and Zimmerman and Adams [2] for particle transport through binary stochastic mixtures. We assess the accuracy of these algorithms using a standard suite of planar geometry incident angular flux benchmark problems and a new suite of interior source benchmark problems. In addition to comparisons of the ensemble-averaged leakage values, we compare the ensemble-averaged material scalar flux distributions. Both Monte Carlo transport algorithms robustly produce physically realistic scalar flux distributions for the benchmark transport problems examined. The base Monte Carlo algorithm reproduces the standard Levermore–Pomraning model [3], [4] results. The improved Monte Carlo algorithm generally produces significantly more accurate leakage values and also significantly more accurate material scalar flux distributions. We also present deterministic atomic mix solutions of the benchmark problems for comparison with the benchmark and the Monte Carlo solutions. Both Monte Carlo algorithms are generally significantly more accurate than the atomic mix approximation for the benchmark suites examined.

A previously described algorithm [T.A. Brunner, T.J. Urbatsch, T.M. Evans, N.A. Gentile, Comparison of four parallel algorithms for domain decomposed implicit Monte Carlo, Journal of Computational Physics 212 (2) (2006) 527–539] for doing domain decomposed particle Monte Carlo calculations in the context of thermal radiation transport has been improved. It has been extended to support cases where the number of particles in a time step are unknown at the beginning of the time step. This situation arises when various physical processes, such as neutron transport, can generate additional particles during the time step, or when particle splitting is used for variance reduction. Additionally, several race conditions that existed in the previous algorithm and could cause code hangs have been fixed. This new algorithm is believed to be robust against all race conditions. The parallel scalability of the new algorithm remains excellent.

The double spherical harmonics angular approximation in the lowest order, i.e. double P0 (DP0), is developed for the solution of time-dependent non-equilibrium grey radiative transfer problems in planar geometry. Although the DP0 diffusion approximation is expected to be less accurate than the P1 diffusion approximation at and near thermodynamic equilibrium, the DP0 angular approximation can more accurately capture the complicated angular dependence near a non-equilibrium radiation wave front. In addition, the DP0 approximation should be more accurate in non-equilibrium optically thin regions where the positive and negative angular domains are largely decoupled. We develop an adaptive angular technique that locally uses either the DP0 or P1 flux-limited diffusion approximation depending on the degree to which the radiation and material fields are in thermodynamic equilibrium. Numerical results are presented for two test problems due to Su and Olson and to Ganapol and Pomraning for which semi-analytic transport solutions exist. These numerical results demonstrate that the adaptive P1–DP0 diffusion approximation can yield improvements in accuracy over the standard P1 diffusion approximation, both without and with flux-limiting, for non-equilibrium grey radiative transfer.

A time-dependent massively-parallel Monte Carlo particle transport calculational module (ParticleMC) for inertial confinement fusion (ICF) applications is described. The ParticleMC package is designed with the long-term goal of transporting neutrons, charged particles, and gamma rays created during the simulation of ICF targets and surrounding materials, although currently the package treats neutrons and gamma rays. Neutrons created during thermonuclear burn provide a source of neutrons to the ParticleMC package. Other user-defined sources of particles are also available. The module is used within the context of a hydrodynamics client code, and the particle tracking is performed on the same computational mesh as used in the broader simulation. The module uses domain-decomposition and the MPI message passing interface to achieve parallel scaling for large numbers of computational cells. The Doppler effects of bulk hydrodynamic motion and the thermal effects due to the high temperatures encountered in ICF plasmas are directly included in the simulation. Numerical results for a three-dimensional benchmark test problem are presented in 3D XYZ geometry as a verification of the basic transport capability. In the full paper, additional numerical results including a prototype ICF simulation will be presented.

We describe an asymptotic analysis of the coupled nonlinear system of equations describing time-dependent three-dimensional monoenergetic neutron transport and isotopic depletion and radioactive decay. The classic asymptotic diffusion scaling of Larsen and Keller [1], along with a consistent small scaling of the terms describing the radioactive decay of isotopes, is applied to this coupled nonlinear system of equations in a medium of specified initial isotopic composition. The analysis demonstrates that to leading order the neutron transport equation limits to the standard time-dependent neutron diffusion equation with macroscopic cross sections whose number densities are determined by the standard system of ordinary differential equations, the so-called Bateman equations, describing the temporal evolution of the nuclide number densities.

A variational analysis is used to derive a mixed P1–DP0 (P1 spherical harmonics–double P0 spherical harmonics) angular approximation to the time-independent monoenergetic neutron transport equation in one-dimensional planar geometry. This mixed angular approximation contains a space-dependent weight factor α(x) that controls the local angular approximation used at a spatial point x: α(x) = 1 yields the standard P1 (diffusion) approximation, α(x) = 0 gives the standard DP0 approximation, and 0 < α(x) < 1 produces a mixed P1–DP0 angular approximation. The diffusion equation obtained differs from the standard P1 diffusion equation only in the definition of the diffusion coefficient. Standard Marshak incident angular flux boundary conditions are also obtained via the variational analysis. We examine the use of this mixed angular approximation coupled with the standard P1 approximation to more accurately treat material interfaces and vacuum boundaries. We propose a simple but effective functional form for the weight factor α(x) that removes the need for the user to specify the value. Numerical results from several test problems are presented to demonstrate that significant improvements in accuracy can be obtained using this method with essentially no computational penalty.