Motivation
I’ve been tasked with implementing Delta Tacking in OpenMC. There are several motivations for this:
- TRISO particle / pebble bed reactors
- Profiling has shown that in OpenMC the surface crossings between layers of TRISO particles or in and out of pebbles in these reactor designs take up the majority of our simulation time.
- GPUs
- Delta-tracking is potentially more performant on GPUs than the ray tracing method as its decision tree is more streamlined.
- Multi-physics coupling
- Coupling volumetric meshes to OpenMC will be significantly easier using Delta tracking as point containment queries on these meshes is far less complex than attempting a ray tracing tracking algorithm through the elements.
Delta Tracking
Most Monte Carlo Radiation Transport (MCRT) codes currently rely on ray tracing or “surface tracking” during particle transport. This tracking algorithm chooses a next-event distance based on the total mean free path, $\Sigma_{t}$, for the particle and its expected exponential attenuation to describe the probability, $p$, that a particle will undergo an interaction with the medium at a distance $x$:
\[p = e^{-\Sigma_{t} x}\]For a sampled probability between zero and one, we can rearrange this function to describe the distance, $x$, a particle has traveled.
\[x_{coll} = x = -\frac{\ln(p)}{\Sigma_{t}}\]This is typically referred to as the distance to collision, $x_{coll}$ . To ensure this collision is truly occurring in the media associated with $\Sigma_{t}$ , $x_{coll}$ is compared to the distance to the nearest surface crossing along the particle’s trajectory, $x_{surf}$ . This then results in two options for the particle’s next event:
- if $x_{coll} \lt x_{surf} \rightarrow$ particle undergoes a collision
- if $x_{coll} \ge x_{surf} \rightarrow$ particle crosses the surface
After either of these occurs, $\Sigma_{t}$ is updated and the process continues.
While ray tracing is still the most commonly implemented form of particle tracking in Monte Carlo codes, a shift toward Delta tracking is rapidly occuring due to it’s need for only one robust geometric operation (point containment) in a code’s geometry kernel. The remainder of this section will summarize the Delta tracking algorithm.
Delta tracking (or Woodcock Delta tracking) was proposed in 1965 as a way to advance particles along their trajectory without halting them at surface crossings as is done in traditional ray tracing tracking (Viitanen & Leppanen, 2015). This is accomplished by playing a statistical game based on what is referred to as the majorant cross section, $\Sigma_{m}$. The majorant cross section is defined as a cross section greater than or equal to any cross section for the media in a given problem that a particle may be traveling through. This majorant is then used to sample incremental particle steps forward rather than the actual total cross section of the medium the particle is traveling in. Because the majorant is at least as large as the largest total cross section in the problem, the particle moves forward a distance less than or equal to the distance to the next collision/event in its history.
At each substep, a true collision will occur with the probability:
\[p = \frac{\Sigma_{t}}{\Sigma_{m}}\]where $\Sigma_{t}$ is the total cross section for the particle’s current location. If no collision occurs, then the particle will continue to the next sub-step. The computational trade here is a point containment search at each sub-step to lose the calculation the nearest surface crossing along the particle’s trajectory and a subsequent point containment search when crossing into the next cell/volume. In terms of the geometric complexity this is generally a win. There are cases in which the particle will take many sub-steps inside the same cell, but caching the last cell/cross section is an easy way to circumvent a full domain search. The last cross section can be cached here too as well if a particle undergoes many virtual collisions.
Outside of geometric considerations, there are other challenges with this tracking method to consider as well.
- Delta tracking cannot provide track-length tally estimation, reducing the FOM per particle in the simulation.
- Generation of reliable but optimal majorant cross sections can be difficult
Item one is fairly clear, there is no way to recover this information because there is no way to determine how far a particle has traveled in a given region without determining surface crossings. While cells and cross sections can be cached, particles taking exceedingly small sub-steps due to an overly conservative majorant will slow the simulation. At the other end of this spectrum, an incorrect majorant cross section will cause bias in the simulation for regions of phase space where it is below the total cross section of the problem’s media. Additionally, media with different (or changing) temperatures take extra consideration when forming a majorant cross section due to thermal effects (Doppler Broadening, density changes, etc.). These effects lead to the majority of complications in forming the majorant cross section and improving on existing methods will be a major focus of our work during Delta tracking implementation in OpenMC.
Developing Temperature Dependent Majorant Cross Sections
SIGMA1
The SIGMA1 method involves preprocessing point-wise data for each unique temperature in the problem. Given the point-wise data, it is relatively easy to form a majorant cross section using the maximal value of all data sets for each bin in the energy segmentation used by the code. This method has fallen out of favor due to the amount of time needed to preprocess the data and the large amount of memory it consumes.
Target Motion Sampling
As we look to define majorant cross sections in OpenMC, we’ll begin with what has been done previously. The Serpent Monte Carlo code has relied heavily on Delta tracking for some time now and has significant experience in defining temperature dependent majorant cross sections. In Serpent, Doppler broadening is modeled using the target motion sampling (TMS) approach in which the relative velocity between the incident neutron and target is sampled from a the Maxwell-Boltzmann distribution of thermal motion on an event-to-event basis.
Maxwell-Boltzmann Distribution:
\[f(v) = \Big(\frac{m}{2 \pi k T}\Big)^{\frac{3}{2}} 4 \pi v^{2} e ^{-\frac{m v^{2}}{2kT}}\]This distribution is defined any velocity from zero to infinity, so restrictions are placed on the minimum and maximum velocity sampled to increase sampling efficiency without biasing the result. The end result of this method is that the behavior of the system will approach what we would expect if we were able to calculate the cross section at the specific temperature of the current medium. In other words, it approaches what we’d expect if we were able to calculate the cross section at a given temperature (i.e. the average behavior of the system). When determining a majorant cross section value however, Serpent then has to account for any value within these energy bounds.
Windowed Multipole
The Multipole format is an analytic representation of cross section data using a set of redisudes and poles in the complex plane (Forget et al., 2013). Values of the cross section are evaluated using the following expression:
\[\sigma(u,T) = \frac{1}{2 u^{2} \sqrt{\xi}} \sum_{j} \Re \Big[ i r_{j,x} \sqrt{\pi}W(z) - \frac{r_{j,x}}{\sqrt{\pi}}C \Big(\frac{p_{j}}{\sqrt{\pi}},\frac{u}{2\sqrt{\xi}}\Big) \Big]\]where
\[u = \sqrt{E}\] \[\xi = \frac{k_{b}T}{4A}\] \[z = \frac{u - p_{j}}{2 \sqrt{\xi}}\] \[C \Big(\frac{p_{j}}{\sqrt{\pi}},\frac{u}{2\sqrt{\xi}}\Big) = 2p_{j} \int_{0}^{\infty} du' \frac{ exp \big[- \frac{(u + u')^{2}}{4\xi} \big]}{ p_{j}^{2} - u'^{2}}\]$W$ is the Faddeeva function and $J$ is the number of pole-residue pairs in the cross section representation. Noting the temperature dependence, this form can be evaluated at any temperature to provide a Doppler-broadened value of the cross section on-the-fly. Several thousand pole-residue pairs may be needed to represent a cross section, however, and the evaluation of this function for each one is prohibitively expensive during simulation. Out of this need came the Windowed Multipole format (Josey, 2015).
The windowed multipole format segments the incident neutron energy spectrum into windows where each window represents local pole-residue pairs exactly using the form above, while pairs outside the window are represented using a curve-fit. This advancement greatly reduces the computational cost of evaluating the cross section at a marginal increase in memory and a negligible loss in accuracy.
Based on the analytic form of the windowed multipole format, it is possible to determine a priori the maximum value of a nuclide’s cross section across the set of temperatures present in the problem. Here is a demonstration of a majorant calculation for U238’s 6.7 eV resonance:
These majorants were generated using OpenMC’s Python API and the SciPy toolkit to find the maximum across the temperature ranges shown. Note that the transitions between different temperatures are not L1 discontinuous, but transition smoothly from curve to curve. This is because the maximum is being found on a continuous temperature spectrum.
For this resonance, the result is fairly predictable and majorant values could be determined based on a relatively simple maximum location scheme, but other regions of the spectrum are not as clear. Interplay between adjacent or overlapping resonances can have unpredictable effects on the majorant cross section value - potentially creating multiple maxima and minima in the temperature range of the problem. In this case, a simple global maximum search may fail and an invalid or inefficient majorant will be returned. Such an example can be seen here in the 19-20 eV range for U238.
Here, we can see some cross-talk between the resonances causing a combined rise in the cross section between the two.
It was originally hoped that majorants could be pre-processed for the windowed multipole format and the rumtime processing of majorants avoided altogether, but this would still result in significant inefficiencies near the peak of resonances if a 0K base temperature were always used. However, trends in the resulting majorant temperature profile using various base temperatures indicate that it may be possible to customize the majorant for a problem’s unique temperature range.
References
- Viitanen, T., & Leppanen, J. (2015). Temperature majorant cross sections in Monte Carlo neutron tracking. Nuclear Science and Engineering, 180(2), 209–223. https://doi.org/10.13182/NSE14-46
- Forget, B., Xu, S., & Smith, K. (2013). Direct Doppler broadening in Monte Carlo simulations using the multipole representation. https://doi.org/10.1016/j.anucene.2013.09.043
- Josey, C. (2015). Windowed Multipole-An Efficient Doppler Broadening Technique for Monte Carlo by Signature redacted.