© P.A. Vaquer et al., Published by EDP Sciences, 2025

1. Introduction

Accuracy is paramount in the field of nuclear engineering and Monte Carlo radiation transport software, such as the MCNP^®1 code, are often considered the “gold standard” for accurate nuclear simulations [1]. The key reason is that Monte Carlo methods rely on very few assumptions about particle transport and geometry. Monte Carlo codes replicate the random nature of nuclear interactions by tracing individual particle tracks through a geometry, and they directly sample from fundamental cross-section probability distributions, which are often represented at high resolution (i.e. continuous-energy cross sections). They provide the flexibility to model even the most complex nuclear systems, and can be used in multi-physics frameworks by inputting/outputting field quantities onto a common geometric mesh that is shared with other physics software.

One of the drawbacks of current Monte Carlo simulations is that they often use the same nuclear cross-section data over large regions of space, neglecting cross-section differences caused by spatial variations in material density and temperature. This is analogous to simulate the movement of a cue ball on a worn-out billiard table without considering any local imperfections of the table.

In some cases, it is acceptable to ignore spatial variations in cross section because there may not be much variation in material densities or temperatures. However, for applications such as fusion reactor simulations, this simplification is inadequate. Plasma densities and temperatures in a fusion reactor can vary by orders of magnitude over short distances, requiring a much better representation of continuously varying cross sections.

Fig. 1. Comparison of element-wise-constant to linear finite-element cross sections, where color is used to represent cross-section magnitude for a homogeneous material. (a) Element-wise-constant cross sections. (b) Linear finite-element cross sections.

This paper presents a first-order-in-space method for modeling continuously varying cross sections in 3-D mesh-based Monte Carlo neutral-particle transport simulations. This method is based on Brown’s direct-sampling approach, which was previously applied to 1-D problems [2]. However, this method defines cross sections as linear nodal finite elements, such that cross-section values are specified on the vertices of a tetrahedral mesh. Monte Carlo particle tracks are then sampled directly from these finite-element cross sections. This method is primarily intended for multi-physics frameworks, such that linear finite-element material properties can be directly imported into Monte Carlo codes.

In 1972, Carter, Cashwell, and Taylor were likely the first to model continuously varying cross sections in Monte Carlo codes [3]. Specifically, they implemented a delta-tracking method that does not require knowledge of the majorant in a modified version of the MCN code (a predecessor to the MCNP^® code) [1, 4, 5]. Carter, Cashwell, and Taylor considered direct sampling to be too computationally expensive because it requires numerical integration, ∫₀^sΣ(s′)ds′, to determine the number of mean free paths to the collision site. Belanger also critiqued the direct sampling method because generally it requires iterative methods to determine a particle’s distance to collision [6]. Fortunately, these computational burdens can be avoided when using linear simplicial finite elements to model cross section. The method derived in this paper only requires solving a quadratic equation to determine a particle’sdistance to collision, circumventing the need for numerical integration or iterative methods.

The main advantage of this direct sampling approach over delta tracking is that the track length in each mesh element is known, thus enabling the use of track-length tallies. Track-length tallies are particularly useful in regions of space with few collisions, where collision-based tallies would contain much larger statistical uncertainties. For example, using direct sampling along with track-length tallies could be particularly advantageous for estimating particle flux at various altitudes in Earth’s atmosphere, because at higher altitudes, density decreases dramatically, resulting in a minimal number of collisions.

The advantage of using linear finite elements instead of element-wise-constant cross sections is illustrated by Figure 1. For linear nodal finite elements, the local cross section at any point within an element can be linearly interpolated based on the cross-section values defined on the element’s vertices. In contrast, element-wise-constant cross sections require averaging the cross section over an element, which does a poorer job of preserving cross section on a point-wise basis. Although the plots in Figure 1 are 2-D, the methods in this paper are applied to 3-D tetrahedral meshes.

Note, this paper is an extended version of a conference paper published for SNA + MC 2024 [7]. Unlike the original conference paper, this paper provides more details about the derivation of the direct sampling method, discusses linearization strategies, includes a wider variety of test problems, and compares the results to three different element-wise-constant cross-section models.

2. Theory for direct-sampling algorithm

This section demonstrates how to directly sample Monte Carlo flight paths in 3-D tetrahedral meshes with linear nodal finite-element cross sections.

2.1. Linear nodal tetrahedral finite element

As previously mentioned, the cross section at any point within a linear nodal tetrahedral finite element can be determined by linear interpolation of the element’s vertex values. This means the equation for the linear nodal tetrahedral finite-element cross section is

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\overrightarrow{r})=\sum _{j=1}^4{\mathrm{\Sigma }}_j{\lambda }_j(\overrightarrow{r}),\end{array}$$(1)

where

2.2. Tetrahedral barycentric coordinates

Barycentric coordinates are a convenient way to describe the location of a point relative to other points. Each vertex in a tetrahedron has a corresponding barycentric coordinate, so there are a total of four barycentric coordinates,

$$\begin{array}{c}\hfill \lambda =({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4).\end{array}$$(2)

The barycentric coordinates corresponding to the 1st vertex are λ = (1, 0, 0, 0). Similarly, the barycentric coordinates corresponding to the 2nd, 3rd, and 4th vertex are λ = (0, 1, 0, 0), λ = (0, 0, 1, 0) and λ = (0, 0, 0, 1), respectively. Furthermore, the sum of the barycentric coordinates must always equal one. Figure 2 demonstrates how tetrahedral barycentric coordinates describe the location of a point relative to the tetrahedron’s vertices.

Fig. 2. The barycentric coordinate λj at a point is determined by first constructing a sub-tetrahedron (the translucent gray region) composed of the point of interest and all the vertices except the jth vertex. The volume of this sub-tetrahedron is then divided by the total tetrahedral volume formed by the four vertices. Notice that as the point of interest moves closer to the jth vertex, then λj approaches one, and as the point moves further away λj approaches zero.

Tetrahedral barycentric coordinates can be described as volume ratios,

$$\begin{array}{c}\hfill {\lambda }_j=\frac{V_j}{V_{}\text{tot}}\text{for}j\in \}1,2,3,4\{.\end{array}$$(3)

where

The total volume of the tetrahedral element can be calculated using the determinant,

$$\begin{array}{c}\hfill V_{}\text{tot}=\frac{1}{6}det\left(\begin{array}{cccc}1& 1& 1& 1\\ x_1& x_2& x_3& x_4\\ y_1& y_2& y_3& y_4\\ z_1& z_2& z_3& z_4\end{array}\right).\end{array}$$(4)

where (x_j, y_j, z_j) are the Cartesian coordinates of the jth vertex. The total volume can be expanded to

$$\begin{array}{c}\hfill V_{}\text{tot}=\frac{1}{6}({\delta }_1+{\delta }_2+{\delta }_3+{\delta }_4)\end{array}$$(5)

where

$$\begin{array}{cc}\hfill {\delta }_1& =x_2{y}_3{z}_4-x_2{y}_4{z}_3-x_3{y}_2{z}_4\hfill \\ \hfill & +x_3{y}_4{z}_2+x_4{y}_2{z}_3-x_4{y}_3{z}_2,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\delta }_2& =-x_1{y}_3{z}_4+x_1{y}_4{z}_3+x_3{y}_1{z}_4\hfill \\ \hfill & -x_3{y}_4{z}_1-x_4{y}_1{z}_3+x_4{y}_3{z}_1,\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\delta }_3& =x_1{y}_2{z}_4-x_1{y}_4{z}_2-x_2{y}_1{z}_4\hfill \\ \hfill & +x_2{y}_4{z}_1+x_4{y}_1{z}_2-x_4{y}_2{z}_1,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill {\delta }_4& =-x_1{y}_2{z}_3+x_1{y}_3{z}_2+x_2{y}_1{z}_3\hfill \\ \hfill & -x_2{y}_3{z}_1-x_3{y}_1{z}_2+x_3{y}_2{z}_1.\hfill \end{array}$$(9)

The volumes of the four sub-tetrahedra, V₁, V₂, V₃, and V₄, can be expressed as

$$\begin{array}{c}\hfill V_j=\frac{1}{6}({\alpha }_{j}x+{\beta }_{j}y+{\gamma }_{j}z+{\delta }_j)\text{for}j\in \}1,2,3,4\{,\end{array}$$(10)

where (x, y, z) are the Cartesian coordinates which are being converted into barycentric coordinates, and

$$\begin{array}{c}\hfill {\alpha }_1=-y_2{z}_3+y_2{z}_4+y_3{z}_2-y_3{z}_4-y_4{z}_2+y_4{z}_3,\end{array}$$(11)

$$\begin{array}{c}\hfill {\alpha }_2=y_1{z}_3-y_1{z}_4-y_3{z}_1+y_3{z}_4+y_4{z}_1-y_4{z}_3,\end{array}$$(12)

$$\begin{array}{c}\hfill {\alpha }_3=-y_1{z}_2+y_1{z}_4+y_2{z}_1-y_2{z}_4-y_4{z}_1+y_4{z}_2,\end{array}$$(13)

$$\begin{array}{c}\hfill {\alpha }_4=y_1{z}_2-y_1{z}_3-y_2{z}_1+y_2{z}_3+y_3{z}_1-y_3{z}_2,\end{array}$$(14)

$$\begin{array}{c}\hfill {\beta }_1=x_2{z}_3-x_2{z}_4-x_3{z}_2+x_3{z}_4+x_4{z}_2-x_4{z}_3,\end{array}$$(15)

$$\begin{array}{c}\hfill {\beta }_2=-x_1{z}_3+x_1{z}_4+x_3{z}_1-x_3{z}_4-x_4{z}_1+x_4{z}_3,\end{array}$$(16)

$$\begin{array}{c}\hfill {\beta }_3=x_1{z}_2-x_1{z}_4-x_2{z}_1+x_2{z}_4+x_4{z}_1-x_4{z}_2,\end{array}$$(17)

$$\begin{array}{c}\hfill {\beta }_4=-x_1{z}_2+x_1{z}_3+x_2{z}_1-x_2{z}_3-x_3{z}_1+x_3{z}_2,\end{array}$$(18)

$$\begin{array}{c}\hfill {\gamma }_1=-x_2{y}_3+x_2{y}_4+x_3{y}_2-x_3{y}_4-x_4{y}_2+x_4{y}_3,\end{array}$$(19)

$$\begin{array}{c}\hfill {\gamma }_2=x_1{y}_3-x_1{y}_4-x_3{y}_1+x_3{y}_4+x_4{y}_1-x_4{y}_3,\end{array}$$(20)

$$\begin{array}{c}\hfill {\gamma }_3=-x_1{y}_2+x_1{y}_4+x_2{y}_1-x_2{y}_4-x_4{y}_1+x_4{y}_2,\end{array}$$(21)

$$\begin{array}{c}\hfill {\gamma }_4=x_1{y}_2-x_1{y}_3-x_2{y}_1+x_2{y}_3+x_3{y}_1-x_3{y}_2.\end{array}$$(22)

Next, by substituting Equation (10) into Equation (3),

$$\begin{array}{c}\hfill {\lambda }_j=\frac{1}{6V_{}\text{tot}}({\alpha }_{j}x+{\beta }_{j}y+{\gamma }_{j}z+{\delta }_j)\text{for}j\in \}1,2,3,4\{.\end{array}$$(23)

Note, Equations (6)–(9) and Equations (11)–(22) are similar to corresponding equations in the SNA + MC 2024 paper [7], but multiplied by a factor of ( − 1), yet they still result in the same values for m_j and ℓ_j in Equations (25) and (26). A Python implementation of Equations (6)–(9) and Equations (11)–(22) is shown in Appendix A.

2.3. Track parametrization of barycentric coordinates

Barycentric coordinates can be parametrized as linear functions corresponding to a particle’s distance from its starting location along a predefined track,

$$\begin{array}{c}\hfill {\lambda }_j({\overrightarrow{r}}_o+s^{\prime }\widehat{\mathrm{\Omega }})=m_j{s}^{\prime }+{\ell }_j\text{for}j\in \}1,2,3,4\{,\end{array}$$(24)

where

By evaluating Equation (23) at (x_o + s′Ω_x,  y_o + s′Ω_y z_o+s′Ω_z) and applying Equation (24) to the left-hand side of Equation (23),

$$\begin{array}{c}\hfill \left(\begin{array}{c}{m}_1{s}^{\prime }+{\ell }_1\\ m_2{s}^{\prime }+{\ell }_2\\ m_3{s}^{\prime }+{\ell }_3\\ m_4{s}^{\prime }+{\ell }_4\end{array}\right)=\frac{1}{6V_{}\text{tot}}\left(\begin{array}{cccc}{\alpha }_1& {\beta }_1& {\gamma }_1& {\delta }_1\\ {\alpha }_2& {\beta }_2& {\gamma }_2& {\delta }_2\\ {\alpha }_3& {\beta }_3& {\gamma }_3& {\delta }_3\\ {\alpha }_4& {\beta }_4& {\gamma }_4& {\delta }_4\end{array}\right)\left(\begin{array}{c}{x}_o+s^{\prime }{\mathrm{\Omega }}_x\\ y_o+s^{\prime }{\mathrm{\Omega }}_y\\ z_o+s^{\prime }{\mathrm{\Omega }}_z\\ 1\end{array}\right).\end{array}$$

Next, the m_j coefficients and ℓ_j constants can be evaluated separately,

$$\begin{array}{c}\hfill m_j=\frac{{\alpha }_j{\mathrm{\Omega }}_x+{\beta }_j{\mathrm{\Omega }}_y+{\gamma }_j{\mathrm{\Omega }}_z}{6V_{}\text{tot}}\text{for}j\in \}1,2,3,4\{,\end{array}$$(25)

$$\begin{array}{c}\hfill {\ell }_j=\frac{{\alpha }_j{x}_o+{\beta }_j{y}_o+{\gamma }_j{z}_o+{\delta }_j}{6V_{}\text{tot}}\text{for}j\in \}1,2,3,4\{.\end{array}$$(26)

2.4. Track parametrization of cross section

Similarly, cross sections can be parametrized as a linear function corresponding to a particle’s distance from its starting location along a predefined track. This is accomplished by substituting Equation (24) into Equation (1),

$$\begin{array}{c}\hfill \mathrm{\Sigma }({\overrightarrow{r}}_o+s^{\prime }\widehat{\mathrm{\Omega }})=\sum _{j=1}^4{\mathrm{\Sigma }}_j(m_j{s}^{\prime }+{\ell }_j).\end{array}$$(27)

Next, by introducing the variables a and b, this can be further simplified to

$$\begin{array}{c}\hfill \mathrm{\Sigma }({\overrightarrow{r}}_o+s^{\prime }\widehat{\mathrm{\Omega }})=2a{s}^{\prime }+b,\end{array}$$(28)

where

$$\begin{array}{c}\hfill a=\frac{1}{2}\sum _{j=1}^4{\mathrm{\Sigma }}_j{m}_j,\end{array}$$(29)

$$\begin{array}{c}\hfill b=\sum _{j=1}^4{\mathrm{\Sigma }}_j{\ell }_j.\end{array}$$(30)

A Python implementation of Equations (25), (26), (29), and (30) is shown in Appendix A.

Until now, Sections 2.1–2.4 have only discussed how cross sections can be expressed as linear nodal finite elements and how cross sections can be parametrized as a linear function of distance along a track. This theory can be applied to a variety of particle-tracking algorithms, such as delta tracking [4]. However, Sections 2.5–2.9 will focus on the direct-sampling algorithm.

2.5. Optical depth

The “optical depth” is the number of mean free paths to a collision site along a predefined particle track. However, in continuously varying media, the concept of mean free path no longer refers to a single, fixed length. Instead, the mean free path depends on the local cross section, so the total optical depth is found by summing (integrating) many small segments of “fractional” mean free paths along a particle track,

$$\begin{array}{c}\hfill \tau ({\overrightarrow{r}}_o\to {\overrightarrow{r}}_o+s\widehat{\mathrm{\Omega }})=\underset{0}{\overset{s}{\int }}\mathrm{\Sigma }({\overrightarrow{r}}_o+s^{\prime }\widehat{\mathrm{\Omega }})\mathrm{d}{s}^{\prime },\end{array}$$(31)

where

Furthermore, by analytically integrating Equation (28), the optical depth between ${\overrightarrow{r}}_o$ and ${\overrightarrow{r}}_o+s\widehat{\mathrm{\Omega }}$ can be simplified to

$$\begin{array}{c}\hfill \tau ({\overrightarrow{r}}_o\to {\overrightarrow{r}}_o+s\widehat{\mathrm{\Omega }})=a{s}^2+bs.\end{array}$$(32)

2.6. Probability of no collision within an element

Before determining the probability of no collision within an element, one must first compute the distance to the next element face, s_f, as shown in Figure 3a. This can be accomplished using a ray-triangle intersection algorithm, such as the Möller-Trumbore algorithm [8].

Fig. 3. An example of computing a particle’s distance to the next event. The variation in color represents the linear variation in cross section within the element.

The probability of no collision, P_NC, between the current location, ${\overrightarrow{r}}_o$, and the next element, ${\overrightarrow{r}}_o+s_f\widehat{\mathrm{\Omega }}$, along a predefined track, is

$$\begin{array}{c}\hfill P_{}\text{NC}=exp(-{\tau }_f),\end{array}$$(33)

where ${\tau }_f\equiv \tau ({\overrightarrow{r}}_o\to {\overrightarrow{r}}_o+s_f\widehat{\mathrm{\Omega }})$. Next, by substituting Equation (32) into Equation (33), the probability of no collision is simply

$$\begin{array}{c}\hfill P_{}\text{NC}=exp(-a{s}_f^2-b{s}_f).\end{array}$$(34)

2.7. Determine the optical depth to the collision site

After determining the probability of no collision, a random number, ξ, is generated from a uniform distribution between 0 and 1 to simulate whether a collision occurs within the current element,

$$\begin{array}{c}\hfill \{\begin{array}{c}\text{if}\xi \le P_{}\text{NC},\text{particle exits element without colliding}\hfill \\ \text{if}\xi >P_{}\text{NC},\text{particle experiences collision before exiting}.\hfill \end{array}\end{array}$$

If ξ >  P_NC, then the distance to collision, s_c, must be determined as shown in Figure 3b. This is accomplished by first determining the optical depth to the collision site.

The optical depth to the collision site, τ_c, is computed by randomly sampling from the cumulative distribution function (CDF). Generally, the region of interest is 0 ≤ τ_c <  ∞ and the CDF is F(τ_c)=1 − exp(−τ_c). However, if the optical depth is confined to 0 ≤ τ_c <  τ_f, then the CDF must be normalized to

$$\begin{array}{c}\hfill F({\tau }_c)=\frac{1-exp(-{\tau }_c)}{1-P_{}\text{NC}}.\end{array}$$(35)

Now that the CDF is known, a second random number is needed to uniformly sample from the CDF. Conveniently, the previous random number ξ can be recycled because at this point it is known to satisfy P_NC <  ξ ≤ 1, and therefore using the number (ξ − P_NC)/(1 − P_NC) is equivalent to sampling from a uniform distribution between 0 and 1. By setting this random number equal to the right-hand side of Equation (35),

$$\begin{array}{c}\hfill \frac{\xi -P_{}\text{NC}}{1-P_{}\text{NC}}=\frac{1-exp(-{\tau }_c)}{1-P_{}\text{NC}},\end{array}$$(36)

and simplifying,

$$\begin{array}{c}\hfill {\tau }_c=-ln[1-(\xi -P_{}\text{NC})].\end{array}$$(37)

2.8. Determine the distance to collision

By evaluating Equation (32) at s_c and setting τ_c = −c, the following quadratic equation can be derived

$$\begin{array}{c}\hfill a{s}_c^2+b{s}_c+c=0.\end{array}$$(38)

The two roots of the quadratic equation are

$$\begin{array}{c}\hfill s_c=\frac{-b\pm \sqrt{b^2-4ac}}{2a},\end{array}$$(39)

where the true distance to collision, s_c, is the smallest positive root.

For the direct-sampling algorithm, Equation (39) will only yield real solutions. The reason is that the distance to the next element face, s_f, and cross section are both real and non-negative (s_f, Σ ∈ ℝ_≥0). Additionally, according to Equations (31) and (32), the optical depth to the next face, τ_f = as_f² + bs_f, must also be real and non-negative (τ_f ∈ ℝ_≥0). This means that the discriminant, b² + 4aτ_f, must satisfy b² + 4aτ_f ≥ 0 (because a negative discriminant would yield complex solutions). Furthermore, because τ_c ≤ τ_f, then b² + 4aτ_c ≥ 0 must also be true, and therefore Equation (39) cannot have complex solutions.

A Python implementation of Equation (39) is shown towards the end of Appendix A.

2.9. Algorithm

This direct-sampling algorithm is performed using element-by-element tracking, such that when a particle enters a new element the material properties for that element are queried to determine a new distance to collision. Each element is considered in isolation, as depicted in Figure 3. Algorithm 1 summarizes the direct-sampling procedure.

3. Results

Seven analytical test problems are used to verify the direct-sampling algorithm. Each problem consists of a purely absorbing medium (Σ_t ≡ Σ_a) and a unique function that specifies the cross-section variation along the x dimension. A brief description of each test case is provided in Table 1.

Fig. 4. Four tetrahedral meshes, from coarsest (a) to finest (d), generated with the Cubit software, and used for all test cases. (a) x-dimension mesh size = 2 (6 elements), (b) x-dimension mesh size = 1 (12 elements), (c) x-dimension mesh size = 0.5 (24 elements), (d) x-dimension mesh size = 0.25 (48 elements).

These test problems originate from Brown’s paper [2]. The original test cases simulated Monte Carlo transport of particles through a one-dimensional slab with a thickness of 2 units. This study converts the 1-D slabs into equivalent 3-D prisms with the following dimensions,

$$\begin{array}{c}\hfill 0\le x\le 2,0\le y\le 1,0\le z\le 1.\end{array}$$

The Cubit^®2 software is used to generate unstructured tetrahedral meshes for the rectangular prism for all test cases [9]. For each test case, the mesh is refined along x to verify that simulation error decreases as the mesh is refined. Figure 4 shows the four meshes used for all of the test cases.

Each test case irradiates the left face (x = 0) with a unidirectional beam of particles traveling in the +x direction. The analytical solution for the flux in the purely absorbing medium is

$$\begin{array}{c}\hfill \varphi (x)=exp(-\underset{0}{\overset{x}{\int }}\mathrm{\Sigma }(x^{\prime })\mathrm{d}{x}^{\prime })\text{for}0\le y,z\le 1.\end{array}$$(40)

Therefore, the transmission, T, at x = 2 is

$$\begin{array}{c}\hfill T=exp(-\underset{0}{\overset{2}{\int }}\mathrm{\Sigma }(x)\mathrm{d}x)\text{for}0\le y,z\le 1,\end{array}$$(41)

and the volume-averaged flux, $\overline{\varphi }$, in the rectangular prism is

$$\begin{array}{c}\hfill \overline{\varphi }=\frac{1}{2}\underset{0}{\overset{2}{\int }}exp(-\underset{0}{\overset{x}{\int }}\mathrm{\Sigma }(x^{\prime })\mathrm{d}{x}^{\prime })\mathrm{d}x.\end{array}$$(42)

Table 2 lists the analytical transmission and volume-averaged flux for the seven test cases.

3.1. Linearization

Since the methods in this paper rely on linear nodal finite elements, non-linearly varying cross sections must be “linearized” prior to simulating particle transport. One way to linearize a non-linear function is to simply evaluate the non-linear function at each vertex in the tetrahedral mesh, and assume the function varies linearly in between the vertices. However, for an exponentially decreasing cross section (e.g. test case 4) or an exponentially increasing cross section (e.g. test case 5), this simple linearization strategy will consistently overestimate cross-section values.

An improved strategy is to treat linearization as a constrained-minimization problem

$$\begin{array}{c}\hfill min\Vert \tilde{\mathrm{\Sigma }}-\mathrm{\Sigma }\Vert \text{subject to}\tilde{\mathrm{\Sigma }}\ge 0,\end{array}$$(43)

where

In this study, a second constraint is added to the minimization problem, so that the linear model preserves the integrated cross section,

$$\begin{array}{c}\hfill min\Vert \tilde{\mathrm{\Sigma }}-\mathrm{\Sigma }\Vert \text{subject to}\{\begin{array}{c}\tilde{\mathrm{\Sigma }}\ge 0,\hfill \\ [6pt]\int \tilde{\mathrm{\Sigma }}(x)dx=\int \mathrm{\Sigma }(x)dx.\hfill \end{array}\end{array}$$(44)

This linearization strategy is performed in 1-D using the |scipy.optimize.minimize| function, and the integrals are computed using the |scipy.integrate.quad| function of the SciPy Python package [10]. The results are then mapped onto the corresponding 3-D tetrahedral mesh.

Figure 5 shows how the continuously varying cross section in test cases 4 through 7 are linearized (note: test cases 1 through 3 do not require linearization).

Fig. 5. Linearized cross sections for test cases 4 through 7, for mesh sizes of 2 and 0.5. (a) Linearizations for test case 4, (b) Linearizations for test case 5, (c) Linearizations for test case 6, (d) Linearizations for test case 7.

Fig. 6. Mesh convergence results for test case 1 (constant: Σ(x)=1). Smaller values are better. The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 7. Mesh convergence results for test case 2 (linearly decreasing: Σ(x)=2 − x). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 8. Mesh convergence results for test case 3 (linearly increasing: Σ(x)=x). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 9. Mesh convergence results for test case 4 (exponentially decreasing: Σ(x)=exp(−3x)). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 10. Mesh convergence results for test case 5 (exponentially increasing: Σ(x)=0.1exp(2x)). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 11. Mesh convergence results for test case 6 (sharp Gaussian: $\mathrm{\Sigma }(x)=\frac{40}{\sqrt{\pi }}exp(-{\left(\frac{x-1}{0.05}\right)}^2)$). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

Fig. 12. Mesh convergence results for test case 7 (broad Gaussian: $\mathrm{\Sigma }(x)=\frac{2}{\sqrt{2\pi }}exp(-{(x-1)}^2)$). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.

3.2. Element-wise-constant cross sections

The accuracy of the linear finite-element model for cross section is compared to three simpler piecewise-constant (element-wise-constant) models for cross section: a “conservative” model, an “averaged” model, and a “centroid” model.

3.2.1. The conservative element-wise-constant model

The “conservative” element-wise-constant model preserves the integral of the continuously varying cross section over a tetrahedral element.

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_i=\frac{1}{V_i}\underset{V_i}{\int }\mathrm{\Sigma }(\overrightarrow{r})\mathrm{d}V,\end{array}$$(45)

where

The integral over each tetrahedron is performed using 10th-order Gaussian quadrature.

3.2.2. The averaged element-wise-constant model

The “averaged” element-wise-constant model simply averages the cross-section values for the four vertices of the already linearized cross-section model.

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_i=\frac{1}{4}\sum _{j=1}^4{\tilde{\mathrm{\Sigma }}}_{i,j},\end{array}$$(46)

where

3.2.3. The centroid element-wise-constant model

The “centroid” element-wise-constant model simply evaluates the continuously varying cross section at the centroid (midpoint) of each tetrahedral element,

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_i=\mathrm{\Sigma }({\overrightarrow{r}}_{i,\text{centroid}}).\end{array}$$(47)

where

3.3. Mesh convergence study

For the mesh convergence study, the meshes are generated by setting the mesh size for x dimension in the Cubit software to 2, 1, 0.5, and 0.25, which results in meshes composed of 6, 12, 24, and 48 tetrahedral elements, respectively (as shown in Fig. 4). Each mesh is generated independently rather than subdividing a coarser one. The Cubit commands for generating the finest mesh are provided in Appendix B.

A Python package is used to import the Cubit meshes and perform direct sampling through the tetrahedral linear finite-element mesh, as described by Algorithm 1. Each simulation consists of 10 000 000 particle histories.

Figures 6–12 present the mesh-convergence results for all seven test cases. As shown, the linear finite-element model yields the exact particle-transmission result (surface-flux tally at x = 2) for every problem within statistical error bounds. The linear model also generally provides the most accurate result for the volume-averaged flux tally when compared to element-wise-constant approximations.

Test cases 2 and 3 are particularly important for verification because they demonstrate that the linear finite-element model is exact for both transmission and volume-averaged flux tallies when cross sections vary linearly (as shown in Figs. 7 and 8).

Test cases 4–7 are intended to demonstrate mesh convergence as well as compare the linear finite-element model to various element-wise-constant approximations. For these test cases, the linear finite-element model is generally the most accurate, followed by the “conservative”, “averaged”, and, lastly, the “centroid" element-wise-constant approximations, with a few deviations from this overall ranking (as shown in Figs. 9–12). Note that test case 6 is particularly difficult for all the methods because the mesh sizes used in the simulations are rather large compared to the width of the Gaussian distribution; by further refining the mesh, the results in test case 6 (a sharp Gaussian) should approach results similar to test case 7 (a broad Gaussian).

4. Conclusions and future work

This paper extends Brown’s 1-D direct-sampling method to 3-D tetrahedral-mesh-based Monte Carlo simulations [2]. The algorithm uses cross-section values at vertices of a tetrahedral mesh to construct a linear finite-element representation of cross section within each element. This enables the cross section to be parametrized as a linear function of distance along a predefined track, and its path integral is quadratic. Therefore, the distance to collision can be computed by finding the smallest positive root of the quadratic equation, rather than using numerical integration and iterative methods as in Brown’s original paper [2]. Because many multiphysics workflows already exchange nodal finite-element quantities, this method integrates seamlessly and eliminates the need for re-meshing, interpolation, or averaging of field quantities.

Seven analytical pure-absorber test problems from Brown’s study serve as verification cases [2]. The original 1-D slabs are converted to equivalent 3-D tetrahedral meshes, and each cross-section profile is linearized to preserve the 1-D integrated value. For the first three cases, which have constant or linearly varying cross-sections, the linear finite-element model yields exact transmission and volume-averaged flux within statistical uncertainties. For the remaining four cases, which feature nonlinear cross-section variations, mesh-refinement results show that the model converges toward the analytical solutions. Across all cases, the linear finite-element approach is generally more accurate than element-wise-constant approximations.

These results demonstrate the benefit of using the linear finite-element direct-sampling method instead of element-wise-constant cross-sections in mesh-based Monte Carlo codes. Although the derivation presented here is detailed, integrating the direct-sampling algorithm into an existing mesh-based Monte Carlo code requires only modest source-code changes and can streamline coupling to multiphysics frameworks.

Future work includes further exploring linearization strategies and comparing the computational cost of this direct-sampling method to other delta-tracking and direct-sampling approaches [2, 4, 6, 11, 12]. It is also important to test the method on problems including scattering and fission; both are expected to be straightforward to implement and are essential for realistic simulations.

Acknowledgments

The authors would like to thank the editors of The European Physical Journal: Nuclear Sciences and Technologies (EPJ-N) for improving the quality and presentation of this work.

Funding

This work is supported by the Department of Energy (DOE) Advanced Simulation and Computing (ASC) program at Los Alamos National Laboratory (LANL) operated by Triad National Security, LLC, for the National Nuclear Security Administration (NNSA) under Contract No. 89233218CNA000001.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

Raw simulation data is not available, but can be reproduced using Algorithm 1 (which is partially implemented in Appendix A).

Author contribution statement

Pablo A. Vaquer and Michael E. Rising conceptualized the direct-sampling idea. Pablo A. Vaquer derived the direct-sampling algorithm and linearization strategy, implemented them in Python, and verified them. Joel A. Kulesza, Michael E. Rising, and Colin A. Weaver provided feedback along the way. Pablo A. Vaquer wrote the original draft of the article. Joel A. Kulesza, Michael E. Rising, and Colin A. Weaver edited the final version.

References

Appendix A A subset of Python functions used for the direct-sampling algorithm

A Python function for computing several tetrahedral-geometry parameters using Equations (6)–(9) and Equations (11)–(22) is shown below:

def get_alpha_beta_gamma_delta(vertices):

    # get x,y,z coordinates for 4 vertices
    x1, y1, z1 = vertices[0]
    x2, y2, z2 = vertices[1]
    x3, y3, z3 = vertices[2]
    x4, y4, z4 = vertices[3]

    alpha=[]; beta=[]; gamma=[]; delta=[]

    # alpha_1, alpha_2, alpha_3, alpha_4
    alpha.append(-y2*z3 + y2*z4 + y3*z2 -
                 y3*z4 - y4*z2 + y4*z3)
    alpha.append(y1*z3 - y1*z4 - y3*z1 +
                 y3*z4 + y4*z1 - y4*z3)
    alpha.append(-y1*z2 + y1*z4 + y2*z1 -
                 y2*z4 - y4*z1 + y4*z2)
    alpha.append(y1*z2 - y1*z3 - y2*z1 +
                 y2*z3 + y3*z1 - y3*z2)

    # beta_1, beta_2, beta_3, beta_4
    beta.append(x2*z3 - x2*z4 - x3*z2 +
                x3*z4 + x4*z2 - x4*z3)
    beta.append(-x1*z3 + x1*z4 + x3*z1 -
                x3*z4 - x4*z1 + x4*z3)
    beta.append(x1*z2 - x1*z4 - x2*z1 +
                x2*z4 + x4*z1 - x4*z2)
    beta.append(-x1*z2 + x1*z3 + x2*z1 -
                x2*z3 - x3*z1 + x3*z2)

    # gamma_1, gamma_2, gamma_3, gamma_4
    gamma.append(-x2*y3 + x2*y4 + x3*y2 -
                 x3*y4 - x4*y2 + x4*y3)
    gamma.append(x1*y3 - x1*y4 - x3*y1 +
                 x3*y4 + x4*y1 - x4*y3)
    gamma.append(-x1*y2 + x1*y4 + x2*y1 -
                 x2*y4 - x4*y1 + x4*y2)
    gamma.append(x1*y2 - x1*y3 - x2*y1 +
                 x2*y3 + x3*y1 - x3*y2)

    # delta_1, delta_2, delta_3, delta_4
    delta.append(x2*y3*z4 - x2*y4*z3 -
                 x3*y2*z4 + x3*y4*z2 +
                 x4*y2*z3 - x4*y3*z2)
    delta.append(-x1*y3*z4 + x1*y4*z3 +
                 x3*y1*z4 - x3*y4*z1 -
                 x4*y1*z3 + x4*y3*z1)
    delta.append(x1*y2*z4 - x1*y4*z2 -
                 x2*y1*z4 + x2*y4*z1 +
                 x4*y1*z2 - x4*y2*z1)
    delta.append(-x1*y2*z3 + x1*y3*z2 +
                 x2*y1*z3 - x2*y3*z1 -
                 x3*y1*z2 + x3*y2*z1)

    return alpha, beta, gamma, delta

A Python function for obtaining parametrization parameters using Equations (25), (26), (29), and (30) is shown below:

import numpy as np

def get_a_and_b(
    xs_at_vertices,
    alpha, beta, gamma, delta,
    position, direction):

    x_o, y_o, z_o = position
    omega_x, omega_y, omega_z = direction
    denom = np.sum(delta)

    # compute "m" and "l" arrays
    m = np.zeros(4)
    l = np.zeros(4)
    for i in range (4):
        m[i] = (alpha[i]*omega_x +
                beta[i]*omega_y +
                gamma[i]*omega_z) / denom
        l[i] = (alpha[i]*x_o +
                beta[i]*y_o +
                gamma[i]*z_o +
                delta[i]) / denom

    # compute "a" and "b" parameters
    a = 0.5 * np.sum(xs_at_vertices*m)
    b = np.sum(xs_at_vertices*l)
    return a, b

A Python function for computing the distance to collision using Equations (37) and (39) is shown below:

import numpy as np

def get_dist_to_collision(p_nc, rn, a, b):

    assert rn > p_nc

    c = np.log(1 - (rn-p_nc))
    if np.abs(a) < 1e-14:
        return -c / b

    d = b**2 - 4*a*c # discriminant
    assert d >= 0

    x1 = (-b - np.sqrt(d)) / (2*a)
    x2 = (-b + np.sqrt(d)) / (2*a)
    return min(x for x in [x1, x2] if x>=0)

Appendix B Cubit commands for generating a tetrahedral mesh

The Cubit commands for generating a tetrahedral mesh for a rectangular prism with dimensions 2 × 1 × 1 with a maximum mesh-element-edge length of 0.25 along x are:

brick x 2 y 1 z 1
curve 2 4 6 8 size 0.25
curve 1 3 5 7 9 10 11 12 size 1
htet vol 1

These commands generate the mesh shown in Figure 4d.

All Tables

All Figures

	Fig. 1. Comparison of element-wise-constant to linear finite-element cross sections, where color is used to represent cross-section magnitude for a homogeneous material. (a) Element-wise-constant cross sections. (b) Linear finite-element cross sections.
In the text

Fig. 2. The barycentric coordinate λj at a point is determined by first constructing a sub-tetrahedron (the translucent gray region) composed of the point of interest and all the vertices except the jth vertex. The volume of this sub-tetrahedron is then divided by the total tetrahedral volume formed by the four vertices. Notice that as the point of interest moves closer to the jth vertex, then λj approaches one, and as the point moves further away λj approaches zero.

In the text

	Fig. 3. An example of computing a particle’s distance to the next event. The variation in color represents the linear variation in cross section within the element.
In the text

	Fig. 4. Four tetrahedral meshes, from coarsest (a) to finest (d), generated with the Cubit software, and used for all test cases. (a) x-dimension mesh size = 2 (6 elements), (b) x-dimension mesh size = 1 (12 elements), (c) x-dimension mesh size = 0.5 (24 elements), (d) x-dimension mesh size = 0.25 (48 elements).
In the text

	Fig. 5. Linearized cross sections for test cases 4 through 7, for mesh sizes of 2 and 0.5. (a) Linearizations for test case 4, (b) Linearizations for test case 5, (c) Linearizations for test case 6, (d) Linearizations for test case 7.
In the text

	Fig. 6. Mesh convergence results for test case 1 (constant: Σ(x)=1). Smaller values are better. The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 7. Mesh convergence results for test case 2 (linearly decreasing: Σ(x)=2 − x). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 8. Mesh convergence results for test case 3 (linearly increasing: Σ(x)=x). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 9. Mesh convergence results for test case 4 (exponentially decreasing: Σ(x)=exp(−3x)). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 10. Mesh convergence results for test case 5 (exponentially increasing: Σ(x)=0.1exp(2x)). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 11. Mesh convergence results for test case 6 (sharp Gaussian: $\mathrm{\Sigma }(x)=\frac{40}{\sqrt{\pi }}exp(-{\left(\frac{x-1}{0.05}\right)}^2)$). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text

	Fig. 12. Mesh convergence results for test case 7 (broad Gaussian: $\mathrm{\Sigma }(x)=\frac{2}{\sqrt{2\pi }}exp(-{(x-1)}^2)$). The vertical axis is logarithmic above 10−3 and linear below. Error bars represent one standard deviation of statistical uncertainty.
In the text