Ray tracing combined CSG/CAD computational geometries

2023 journal article
2022 conference paper

One-sentence summary

The Shift Monte Carlo radiation transport code was modified to support computational geometries specified using an arbitrary combination of Constructive Solid Geometry (CSG) and Computer-Aided Design (CAD) models.

Computational Geometries

Monte Carlo (MC) radiation transport simulations are used to assess the safety and performance of nuclear reactors and other nuclear systems. This method involves simulating the paths of particles (i.e., neutrons and/or photons) as they move through a nuclear system. By simulating a large number of particles, statistical methods can be used to estimate radiological quantites such as particle flux, biological dose rate, and the neutron multiplication factor (which determines if a reactor is subcritical, critical, or supercritical).

MC codes track the movement of particles within a computational geometry using ray tracing, the same algorithm used for rendering 3D graphics. The implementation of the ray tracing algorithm depends on how the computational geometry is specificed. Production-level MC codes typically have ray tracing routines that operate on models specified in the following formats:

1. Combinatorial Solid Geometry (CSG): geometry cells are specified by combining geometric primitives (e.g., planes, cuboids, and spheres) using Boolean logic operations, in a text-based format.
2. Computer-Aided Design (CAD): the boundaries of geometry cells are specified using surface meshes with triangular elements.

There are significant human effort and performance trade-offs between these formats, and as a result, both are widely used. However, it can be extremely tedious and time-consuming to combine models in different formats.

Layered Geometry

Shift is a production-level MC code developed at Oak Ridge National Laboratory designed to scale from laptops to supercomputers. A new "layered" geometry type has been implemented in the Shift MC radiation transport code. With layered geometry, CAD and CSG models can be arbitrarily translated, rotated, clipped, and placed in overlapping layers, as seen in the example below. New ray tracing routines have been added to Shift in order to track particles within layered geometries. This drastically reduces the human effort required to combine models in multiple formats.

Demonstration

The layered geometry feature within the Shift MC code was demonstrated with a Transformational Challenge Reactor (TCR) problem. As seen in the diagram below, a layered geometry was created by overlaying complex CAD models of fuel and reflector elements over a CSG core model to form a transport-ready geometry. A k-eigenvalue calculation was carried out on this layered geometry. The converged neutron source and neutron flux distributions appear below.

GPU-friendly algorithm for free gas elastic scattering

2021 conference paper
2021 conference presentation
YouTube video of presentation

One-sentence summary

A new algorithm has been developed for simulating the process of neutrons scaterring off of nuclei, in a way that is conducive to GPU execution.

Singular Value Decomposition (SVD) of Monte Carlo variance reduction parameters

One-sentence summary

The Singular Value Decomposition (SVD) matrix factorization technique is shown to reduce RAM requirements for radiation transport simulation parameters by an order of magnitude, as demonstrated on the Summit supercomputer with a realistic nuclear reactor problem.

SVD basics

Singular Value Decomposition (SVD) is a matrix factorization technique that has applications in data compression and categorization problems. For an $m \times n$ matrix $\mathbf{A}$ the SVD is $$ \mathbf{A}_{m \times n} = \mathbf{U}_{m \times m} \mathbf{\Sigma}_{m \times n} \mathbf{V}_{n \times n}^T, $$ where $\Sigma$ is a diagonal matrix. In this diagonal matrix, the entries decrease monotonically. As a result, an estimate of $A$ can be formed by truncating the $\mathbf{U}$, $\mathbf{\Sigma}$, and $\mathbf{V}$ matrices. Specifically, an estimate of the $a_{ij}$ entry in the $\mathbf{A}$ matrix can be obtained via: $$ a_{ij} \approx \sum_{k=1}^K u_{ik} \sigma_{kk} v_{jk}, $$ with $1 < K < n$. In other words, an estimate of the $\mathbf{A}$ matrix — referred to as a "reconstruction" — can be obtained by storing only the first $K$ rows in $\mathbf{U}$, the first $K$ entries along the diagonal of $\mathbf{\Sigma}$, and the first $K$ columns in $\mathbf{V}$, which results in a factor of $\Gamma$ memory savings, where $\Gamma$ is given by $$ \Gamma = \frac{mn}{mK + K + nK}. $$ For example, consider the following $1333 \times 1000$ pixel image of a cave waterfall in Jackson County, AL: This image can be represented by a $1333 \times 1000$ matrix, where each entry represents a pixel with a gray-scale value between 0 and 1. Three reconstructions of this image are shown below. The $K=100$ reconstruction provides a fairly accurate representation of the image and requires $\Gamma=5.71$ times less memory.

Application of SVD to radiation transport

Monte Carlo (MC) radiation transport calculations are used to assess the safety and performance of nuclear reactors and other nuclear systems. The MC method uses a random walk technique to simulate the paths particles (i.e., neutrons and/or photons) take as they move through a nuclear system, from birth (via fission or an external source) to termination (through absorption reactions or by leaving the spatial boundaries of the system). Particle paths are recorded as they move through regions of interest, known as tallies. By simulating a large number of particles, statistical methods can be used to estimate radiological quantities such as particle flux, biological dose rate, and reactions rates within these tallies.

One significant downside of MC radiation transport is that when tallies are located in regions of low particle flux (i.e., regions that few particles pass through), the convergence rates of tallies may be prohibitively slow. To address this problem most MC codes employ a class of biased random sampling techniques known as MC variance reduction (VR). VR method allow the code to drastically increase the fraction of compute time spent simulating particles that are likely to reach the tally, without introducing systematic bias into the tally results.

Industry-standard VR methods require a low-order estimate of the adjoint flux distribution, which describes the "importance" of position/energy phase space regions to the final result. This information is supplied via "importance maps": 3D spatial meshes with each mesh cell tagged with energy-wise adjoint fluxes. These important maps require a significant amount of RAM, which can limit the fidelity of the simulation that can be achieved.

This work explores the novel application of SVD to compress importance maps for reduced memory requirements. In the example above, the $\mathbf{A}$ matrix represents the two spatial dimensions of the waterfall image. Here, the 3 spatial dimensions of the importance maps are collapsed such that the 2 dimensions of the $\mathbf{A}$ matrix represent the spatial mesh cells and discrete energy groups of the adjoint flux. This allows an $A$ matrix to be formulated as: $$ \mathbf{A} = \begin{bmatrix} \phi^{\dagger \, g = 0}_{c = 0} & \phi^{\dagger \, g = 1}_{c = 0} & \dots & \phi^{\dagger \, g = N_g - 1}_{c = 0} \\[0.3cm] \phi^{\dagger \, g = 0}_{c = 1} & \phi^{\dagger \, g = 1}_{c = 1} & \dots & \phi^{\dagger \, g = N_g - 1}_{c = 1} \\[0.3cm] \vdots & \vdots & \ddots & \vdots \\[0.3cm] \phi^{\dagger \, g = 0}_{c = N_{c} - 1} & \phi^{\dagger \, g = 1}_{c = N_{c} - 1} & \dots & \phi^{\dagger \, g = N_g - 1}_{c = N_{c} - 1} \\[0.3cm] \end{bmatrix}, $$ where $\phi^{\dagger}$ denotes the adjoint flux, $g$ denotes the energy group index, and $c$ denotes the spatial mesh cell index. Since adjoint flux distributions may vary by 10 orders of magnitude or more, and SVD in log-space was also considered, i.e., $$ \mathbf{A} = \begin{bmatrix} \log\left(\phi^{\dagger \, g = 0}_{c = 0}\right) & \log\left(\phi^{\dagger \, g = 1}_{c = 0}\right) & \dots & \log\left(\phi^{\dagger \, g = N_g - 1}_{c = 0}\right) \\[0.3cm] \log\left(\phi^{\dagger \, g = 0}_{c = 1}\right) & \log\left(\phi^{\dagger \, g = 1}_{c = 1}\right) & \dots & \log\left(\phi^{\dagger \, g = N_g - 1}_{c = 1}\right) \\[0.3cm] \vdots & \vdots & \ddots & \vdots \\[0.3cm] \log\left(\phi^{\dagger \, g = 0}_{c = N_{c} - 1}\right) & \log\left(\phi^{\dagger \, g = 1}_{c = N_{c} - 1}\right) & \dots & \log\left(\phi^{\dagger \, g = N_g - 1}_{c = N_{c} - 1}\right) \\[0.3cm] \end{bmatrix}. $$ By storing the truncated SVDs of these $\mathbf{A}$ matricies, rather than the full matrices themselves, a signicant amount of RAM can be saved. Values of the adjoint flux can then be reconstructed on-the-fly (i.e., during a simulation). Both formulations were implemented and performance tested in the Shift Monte Carlo radiation transport code.

Test problem

The Westinghouse AP1000 is a Generation III+ nuclear reactor currently being constructed at several sites in in the US and China. SVD compression of importance maps was tested for the "excore dosimetry" problem shown below. This problem consists of a quarter-core model with reflecting boundaries, with a spherical neutron detector in the excore region (i.e., the region outside the reactor core), on the $z$-midplane. The image on the left shows a $z$-midplane slice of this geometry, with the detector shown in green. The image on the right shows the neutron source distribution from fisison. The objective of the problem is to transport particles from the fission source, through the sheilding (shown in burgundy), to the detector, which serves as a tally. The Denovo 3D discrete ordinates (S$_N$) code was used to obtain an adjoint flux distribution (importance map) on a 44 $\times$ 44 $\times$ 49 spatial mesh with 28 energy groups. The contour plot below shows the adjoint flux for the 1.4227–1.8268 MeV energy group (i.e., fission neutrons). This plot shows that fission neutrons experience more than 6 orders of magnitude of attenuation between the source and detector, which indicates that VR will be necessary. Because 28 energy groups were used, the $\mathbf{A}$ matrix has an $n$ of 28. Shift was run using $K = [1, 28]$ with the SVD performed in both linear and log space. Ten independent trials were run for each $K$. Each trial was run on the Summit supercomputer on 264 nodes with 7 processes per node (due to memory constraints). The plot below shows the speedups and corresponding standard deviations for these trials, where the speedup is defined as the ratio of the convergence rate with VR, to the convergence rate without VR. This plot shows that the two best options are 1) use a rank $K=1$ approximation with linear SVD, which provides a $\Gamma \approx 28*$ while incurring a $\sim3\times$ performance penalty, or 2) use a $K=2$ approximation with a log SVD, which provides a $\Gamma \approx 14$ with no performance penalty. The results show that taking the SVD of the adjoint flux is an effective method for reducing the memory requirements of importance maps without a significant sacrifice in performance.

* Note for the nuclear engineering audience: this work using the CADIS method and these $\,\Gamma$ values refer only to weight windows. For the $K=1$ SVD in linear space, the biased source can also be compressed by a factor of 28, but for the $K=2$ SVD in log space, the biased source cannot be compressed. These detailed were omitted for brevity. A full explanation is found in the journal article .

Optimizing gamma dose rate calculations for nuclear fusion shielding

2017 journal article
2016 dissertation

One-sentence summary

Stocastic radiation transport simulations are optimized for an important nuclear fusion analysis scenario in which neutrons generate radioisotopes in structural materials, resulting in a prolonged source of potentially harmful gamma radiation.

Theory

One of the challenges of designing and operating nuclear fusion reactors (such as ITER) is the creation of radioactive materials. For most viable nuclear fusion fuels, the process of fusing atomic nuclei releases high-energy neutrons. These neutrons can penetrate through the structural materials of the reactor and convert normally harmless materials such as the iron in stainless steel or the hydrogen in water into radioactive isotopes. Even after the reactor is switched off, these isotopes will persist and emit high-energy photons (gamma rays) as they decay. Some of these photons travel to areas where people may be taking measurements or doing maintenance, so it is necessary to estimate the potential biological dose rate. This is generally done using Monte Carlo radiation transport, which involves simulating the paths that particles take as they move through the reactor and interact with the materials.

From the standpoint of computer processing power, these calculations are not possible without employing a class of techniques known as Monte Carlo variance reduction. Variance reduction techniques are used to optimize Monte Carlo transport by first finding an approximate solution to an adjoint formulation of the problem (though not all techniques are explicitly formulated this way). This is done with respect to a detector response function. In terms of a standard optimization problem, the variance of the detector response function can be thought of as an objective function that should be minimized. For decay photons, defining this detector response function is straightforward, as I demonstrated in earlier work [2017 journal article]. For the neutrons that produce the decay photons, this process is much more complicated.

As a component of my dissertation work I developed and implemented a technique for automatically generating the neutron detector response function in order to optimize Monte Carlo radiation transport for the aforementioned neutron problem. The most difficult part of forming this function was approximating the process of "neutron activation" in which material exposed to neutrons creates a vast network of radioactive products. For example when $^{62}\mathrm{Ni}$ is exposed to neutrons, the following reactions, and many more, occur: $$ {^{62}\mathrm{Ni}} \xrightarrow{\textrm{(n, $\alpha$)}} {^{59}\mathrm{Fe}} \xrightarrow{\beta^-} {^{59}\mathrm{Co}} \begin{array}{ll} \xrightarrow{\textrm{(n, $\gamma$)}} {^{60}\mathrm{Co}} \\ \xrightarrow{\textrm{(n, $\gamma$)}} {^{60\mathrm{m}}\mathrm{Co}} \xrightarrow{\mathrm{IT}} {^{60}\mathrm{Co}} \end{array} $$ In general, composition of materials before and after nuclear irradiation can be related via the matrix exponential equation: $$ \vec{N}(t) = e^{\mathbf{A}t} \vec{N}_0,$$ where $\vec{N}_0$ and $\vec{N}(t)$ are vectors describing the composition of a material before and after neutron irradiation for time $t$. The $\mathbf{A}$ is a matrix of reaction rates that depend on the rate in which neutrons pass through a given area (the neutron flux). I developed a set of approximations that can be used to simplify this matrix exponential such that it becomes a linear operator with respect to neutron flux. The details of this process can be found in my dissertation/journal article linked at the top of this page, but just for fun I will show the final answer for the optimal neutron detector response function ($q^+_n$) that results from the approximation of $\mathbf{A}$: $$ \begin{align} q^+_n(E_n, E_p) = \int_{E_p} \text{d}{E_p} \phi^+_p(E_p) \sum\limits_{c} \lambda_{c, i_c} b_{c, i_c}(E_p) N_{c,1}(0) \Bigg\lbrack \Big\lbrack &\sigma_{c, 1 \rightarrow 2}(E_n) \frac{t_{\text{irr}}^{i_c-1}}{(i_c-1)!} \prod\limits_{j=3}^{i_c} P_{c, j} \Big\rbrack e^{-d_{c, i_c} t_{\text{decay}}}\\ & + \sum_{j=2}^{i_c-1} \Big\lbrack \sigma_{c, 1 \rightarrow 2}(E_n) \frac{t_{\text{irr}}^{j-1}}{(j-1)!} \prod\limits_{k=3}^{j} P_{c, k} \Big\rbrack B_{c,i_c,j}(t_{\text{decay}}) \Bigg\rbrack. \end{align} $$

Demonstration

Here I'll show that my method, which is called Groupwise Transmutation Consistent Adjoint Driven Importance Sampling (GT-CADIS) improves the convergence rate of Monte Carlo neutron transport by a factor of $O(10^2)$ relative to the industry-standard method known as Forward-Weighted (FW)-CADIS. Consider a simple problem shown below consisting of a block of stainless steel (gray) with a U-shaped air-duct, a neutron source (blue) and a region (orange) where we'd like to estimate the biological dose rate from photons after neutron irradiation. Using FW-CADIS, neutron detector response function is shown below. Recall that the variance of this function can be thought of as the guess for the objective function. With GT-CADIS, the neutron detector response function has a much different shape: the intensity is much higher in the regions of stainless steel that border the air duct. Note that only the shape is important here, not the magnitude (or even units!). I generated variance reduction parameters using both detector response functions and tracked the rate of convergence of Monte Carlo simulations using each set. The image below shows the relative error in the neutron flux as a function of neutron transport processor time using GT-CADIS and FW-CADIS variance reduction parameters, as well as no variance reduction. Without variance reduction, the neutron flux converges fastest near the source. With FW-CADIS the neutron flux converges at a roughly uniform rate throughout the problem. With GT-CADIS, the neutron flux converges fastest in the area immediately surrounding the area of interest, as expected. In order to quantify how much faster the simulation converges with GT-CADIS, the Monte Carlo figure of merit was used, as defined by $$ {\text{figure of merit}} = \frac{1}{R^2t_{\text{proc}}}, $$ where R is the relative error of the photon flux in the region of interest (e.g. the orange sphere in the geometry image) and $t_{\text{proc}}$ is the neutron transport processor time. The figure of merit converges to a constant value during a simulation, as seen in the plot below. The ratio of figure of merits is the speedup. GT-CADIS was found to provide speedups of $200 \pm 100$ relative to FW-CADIS and $9 \pm 5 \times 10^5$ relative to simulations without variance reduction. In other words, GT-CADIS drastically reduces the required computational resources necessary for calculating these photon dose rates.