LEACS: lagrange error-adaptive collision sampling for multiphysics coupled Monte carlo transport

1 Introduction

Monte Carlo particle transport has long been regarded as the reference solution method in reactor physics, owing to its ability to faithfully represent continuous-energy nuclear data, arbitrarily complex geometry, and the full phase-space distribution of the neutron population without the condensation and discretization errors inherent to deterministic methods. In its conventional form, particles are tracked through a constructive solid geometry (CSG) or CAD representation in which each cell is assigned a single material composition at a single temperature; cross sections are pre-broadened to that temperature and assumed piecewise-constant along the particle’s flight path, enabling efficient collision distance sampling.

Over the past decade, however, high-performance computing advances have driven neutronics toward tightly coupled multiphysics simulation, in which Monte Carlo transport is iterated with thermal-hydraulics, fuel performance, or structural-mechanics solvers. Monte Carlo codes such as OpenMC, Serpent2, MCNP, RMC (Romano et al., 2015; Leppänen et al., 2025; Werner et al., 2018; Wang, Kan et al., 2024) have all demonstrated coupled calculations with sub-channel and CFD codes. Recent work has pushed further still by integrating Monte Carlo transport with finite-element solvers in the MOOSE framework (Novak et al., 2023).

What these coupled solves make clear is that the temperature and density fields inside a reactor are not piecewise-constant. They vary continuously and often steeply in space, with fuel pellet centerline-to-surface gradients of several hundred Kelvin and coolant density changes of tens of percent across an axial span. One can no longer neglect spatial variations in these fields because the feedback they impart to neutron transport and all resulting tallies is precisely what a high-fidelity coupled calculation is meant to capture. Several methodologies have been developed to address this complication, though no method is without its flaws.

The work described in this paper introduces a new method of solving the problem of coupling Monte Carlo transport with spatially varying material properties, named Lagrange Error-Adaptive Collision Sampling (LEACS). The narrative begins with a discussion of the theory behind collision distance sampling, a synopsis of the main contemporary methods, and a general discussion of the logic behind LEACS and how it is different from existing methods. The following are more detailed descriptions of the methodology, including algorithmic LEACS processes. Finally, two sections devoted to implementation results are included. The first details analytical error estimation, and the second presents the results and computational cost of implementing LEACS in OpenMC.

2 Background

This section first details the theory behind collision sampling in Monte Carlo transport. Next, the complications that arise from multiphysics coupling are discussed, followed by the current state-of-the-art solutions to this challenge. Finally, the LEACS methodology, the thesis of this work, is explained along with the logic behind its mechanisms.

2.1 Collision distance estimation

During Monte Carlo particle transport, a single particle’s history is simulated as detailed in Algorithm 1.

Algorithm 1 Single Particle History in Monte Carlo Transport.

Step 4 of Algorithm 1, the estimation of distance to collision, is the focus of this work. Estimating this collision distance begins with the probability density function (PDF) of the collision along a line, defined in Equation 1 and with sub-definitions in Equations 2–4.where

f s d s = P NC δ s − s b d s + d τ s d s e − τ s d s (1)

P NC = e − τ s b (2)

τ s = ∫ 0 s Σ t T s ′ , N s ′ d s ′ (3)

d τ s d s = Σ t s (4)

Here P NC is the probability of not colliding, s and s b are the distance traveled along the particle’s flight path and the distance to the cell boundary, respectively. The optical depth is denoted by τ, and Σ t is the macroscopic total cross-section. Note that the macroscopic cross-section is a function of position because it depends on temperature T and number density N that are both position-dependent in a realistic system.

By integrating this function, one obtains the cumulative distribution function (CDF), seen in Equation 5.

F s = ∫ 0 s f s d s = P NC s b H s − s b + ∫ 0 s d τ s d s e − τ s d s (5)

Here H is the Heaviside step function. The utility of this formulation is that one can now sample a uniform random variable between 0 and 1, ξ, set it equal to the CDF, and invert out the distance to collision. In traditional Monte Carlo transport, material definitions are constant and distances to boundaries are neglected since there will be a later comparison against the boundary distance (see Step 4 of the particle history process above). This simplifies the logic greatly and reduces collision distance sampling to Equation 6.

ξ = e − Σ t s → s = − l n ξ / Σ t (6)

2.2 The multiphysics problem

When coupling particle transport with a multiphysics calculation, the assumption of constant cross-sections becomes a problem. With changing temperatures and densities, cross-sections become spatially dependent within the original material bounds. This adds mathematical complication, though the collision sampling can be reframed by using the method described by Introduced by Brown (Brown and Martin, 2003) and described by Ellis (Ellis, 2017). After first evaluating P NC and sampling a random number to determine if the collision will happen before the boundary, the CDF can be renormalized to Equation 7.

F s = ∫ 0 s f s d s = ξ = 1 − e − τ s 1 − e − τ s b (7)

While this is closer to a tractable form, it does not resolve the main issue: optical depth now depends on the multiphysics solution, the specific path of the particle through that solution, and the manner in which multiphysics information is transferred.

It is this last component that is most important. For example, typical thermal-hydraulics solvers use finite element methods, which represent a temperature distribution on a mesh of elements made up of nodes with associated shape functions that reconstruct the intra-element temperature variation. In order to apply this information to evaluating Equation 7, some transformation of the data is required. Addressing this problem has led to the creation of many solutions, though they primarily fall into one of three major categories.

2.3 Current options

2.3.1 Delta tracking

This method involves calculating a majorant cross-section Σ maj, the largest cross-section in the geometry. Distance to collision is then sampled using Equation 6, followed by a rejection check against Σ t ( r ) / Σ maj. If the random number is below this value the collision is accepted, otherwise the particle flight continues.

This rejection sampling method was originally created by Woodcock (Woodcock et al., 1965) and has a contemporary implementation with improvements in Serpent 2 (Leppänen, 2013). This method is particularly useful for complex geometries such as those containing TRISO particles as the particle does not need to update cross-sections at material boundaries meaning that a large number of surfaces does not impact computational speed.

The drawbacks of this method are twofold.

• If there is a region with extremely high cross-sections, such as a gadolinium absorber in the system, there can be significant computational overhead from oversampling collisions in moderator regions. A hybrid method would typically help with this, in which too many rejections cause the code to revert to default transport, but in this case there is no viable fallback, because standard transport cannot accommodate spatially varying cross-sections.
• Calculating the true majorant cross-section is complicated when coupling to a simulation that provides continually varying temperature and density ranges. Doppler broadening’s effect on cross-section magnitude is highly energy-dependent, particularly near resonance peaks. Additionally, determining whether to reject the collision requires evaluating the cross-section at the particle’s position. This means there must be a representation of the multiphysics solution that the transport solver can access dynamically.

2.3.2 Adaptive binning

The most recently developed method of handling spatially varying cross-sections during multiphysics comes from a recent addition to the CARDINAL code (Novak et al., 2022). The MoabSkinner user object automatically regenerates the DAGMC (Wilson et al., 2010) geometry by binning the geometry into cells of equal temperature and density (Novak et al., 2024). This results in a more complex geometry but permits use of the simple Equation 6. This avoids any changes to the underlying sampling algorithm, and CARDINAL is already a code meant for neutronic-multiphysics coupling, making the addition particularly well-suited to this application.

There are two cons to this solution, one obvious and one less obvious. The obvious con is that a domain with large gradients and a fine binning tolerance will result in substantially more geometry cells, increasing computational cost during transport.

The less obvious con is that the resolution of the multiphysics representation in the neutronic solve is tied to the mesh that the multiphysics was solved on. The neutronic geometry cannot be subdivided into regions smaller than the single element size. This is unlikely to be an explicit issue, as only tetrahedral meshes can currently take advantage of the adaptive binning process. That presents another problem, however, as tetrahedral meshes are not particularly efficient for solving a geometry like a pin cell that is meters long but only a centimeter wide. A tradeoff emerges where one must sacrifice efficiency on the multiphysics side for the ability to accurately represent the material gradients on the neutronics side. The common choice of using a second-order solution with lower mesh refinement becomes less attractive in this case.

2.3.3 Numerical estimation

The final method is based on numerically estimating the optical depth along the particle flight path by sampling the cross-section at different points and fitting a quadratic function to the shape. This approach began with Brown, Martin, and Greisheimer (Brown and Martin, 2003; Griesheimer, 2005), where nuclide densities were prescribed with function expansions in the radial and axial dimension. These functional expansions are evaluated along the particle’s flight path, then cross-sections are computed and the integrals estimated via Simpson’s rule. This sampling and integration is highly efficient on modern CPUs, and accuracy is only limited by the fit of the basis set to the original data.

The present author has explored the applicability of this method extensively (Kendrick, 2024) and found that the primary limitation of this method’s applicability to multiphysics coupling lies in the representation of functional expansions on unstructured meshes. Functional expansions work well in a single dimension, but ensuring orthogonality for three-dimensional domains is non-trivial. In an unstructured mesh, this orthogonality becomes nearly impossible to guarantee. The domains that the basis sets represent do not perfectly match polyhedral mesh elements, particularly if there has been any deformation (Specht and Novak, 2025).

2.4 Proposed solution: LEACS

Both the adaptive binning and numerical estimation solutions take the original multiphysics solution and convert it into a lower-order form that is accessible to the transport solver. The original finite element solution is converted into constant average cells or subcells, or the solution is fitted to polynomial approximation of the original spatial distribution. Delta tracking is the only method that does not do this, and that’s only because there is no concrete definition for how to retrieve true cross-sections when evaluating pseudo-collisions.

The purpose of this work is to propose an alternate solution that has the ability to preserve nearly all of the original information from the multiphysics solution. Computational cost is minimized by spacing sample points to efficiently recreate the collision probability distributions. This solution is the LEACS method, the Lagrange error-adaptive collision sampling method.

2.4.1 Method basics

The LEACS method is essentially an extension of the numerical estimation method developed by Greisheimer and colleagues (Brown and Martin, 2003; Griesheimer, 2005; Ellis, 2017), where collision distance is estimated by fitting quadratic functions to cross-sections that are sampled along the flight path at runtime. The difference for LEACS is that instead of sampling cross-sections by relying on a precalculated polynomial fit to a material quantity (e.g., density), these quantities are sampled directly from the finite element solution during particle transport. This means accessing the native shape functions from the multiphysics solution, ensuring that the accuracy of the information transfer is limited only by the number of sample points, and is therefore user-controllable.

For the work presented, only temperature data is transferred from the multiphysics solution to the Monte Carlo transport, although density would follow the same logic, if that information is also stored in the same finite element solution. These temperatures are used to generate Doppler broadened cross-sections, which are then used to estimate optical depths and therefore collision distance.

Notably, accessing the finite element solution is not without computational cost. The general process (explicitly detailed in Section 3.2) consists of identifying the element containing the sample point, and then evaluating the element’s shape functions to reconstruct the local solution. Because this process is expensive, each sample point must justify its computational cost, a consideration that is discussed further in the following subsection.

2.4.2 Optimal sample spacing

The largest difference between previous numerical estimation methods and LEACS concerns how sample positions are determined. The default spacing for any piecewise representation of the true cross-section profile would be equidistant spacing, as it is the simplest method to implement. The downside to this is that for cases with large optical depth (high cross-section, large distance, or both), the most likely collision distances will be closest to the initial position. Equidistant refinement adds as many points in the unimportant region (near the boundary) as in the region where the cross-section profile is most critical (near the initial position).

Because sampling the finite element solution is anticipated to be the most computationally expensive component of the transport calculation, each sample location should be optimally chosen. LEACS is premised on the principle that accuracy requires reconstructing the correct shape of the distance to collision CDF, as this is the equation inverted to obtain the collision distance. By making an assumption about how the CDF will likely look, one can select sampling locations that are best suited to rebuilding the true shape.

In the initial phase of the particle’s movement through a material, no temperatures (and therefore cross-sections) have been sampled, but one can obtain an estimation of the CDF shape by using the cell-averaged temperature, consistent with the default OpenMC treatment. The estimated CDF would then be (normalized) Equation 8.

F s = 1 − e − τ s 1 − e − τ s b = 1 − e − Σ t s 1 − e − Σ t s b (8)

Two spacing methods follow from this formulation.

• Equal probability sampling.
• 2nd-order curvature-based sampling.

The most straightforward way to select sample spacing based on the CDF is to simply space each point equally in terms of the CDF value. For example, if the user specifies five samples per crossing, this method would select the beginning point, boundary point, and then the s values that return F ( s ) = 0.25 , 0.5 , 0.75.

In functional form, this is expressed as Equation 9:

s i = − 1 Σ t l n 1 − i N 1 − e − Σ t s b (9)

For second-order curvature-based sampling, one can rely on a theorem from Lagrange and Taylor, hence the name LEACS. Taylor’s theorem has some well-defined behaviors when representing curves as piecewise k-th order polynomials. Associated with this is the Lagrange form of the remainder, derived from the mean value theorem, also known as the Lagrange Error Bound or Taylor’s Inequality, seen in Equation 10.

R k x = f k + 1 ξ L k + 1 ! x − a k + 1 (10)

for some real number ξ L between a and x. For the purposes of this work, ( x − a ) is the spacing between points Δ s, ξ L is some point within the interval [ a , x ], k is the order of the approximation being used, and f here is the CDF F ( s ). LEACS uses a quadratic fit, so k = 2.

The key result from Equation 10 is seen in Equation 11:

error ∝ | F k + 1 s | Δ s k + 1 (11)

Because optimal spacing requires the error to be uniform across all intervals, one can rearrange to Equation 12.

Δ s s ∝ | F k + 1 s | − 1 / k + 1 (12)

A cumulative sample count function takes the form of Equation 13.

G s = ∫ 0 s 1 Δ s s ′ d s ′ = ∫ 0 s | F k + 1 s ′ | 1 / k + 1 d s ′ (13)

The physical meaning of this function is the total number of samples from 0 to s. Setting G ( s i ) = i and normalizing by the total sum gives Equation 14.

G s i = i N G s b (14)

which allows one to finally solve for s i. For the quadratic fit employed by LEACS when estimating optical depth, the optimal spacing function is given by Equation 15.

s i = − 3 Σ t l n 1 − i N 1 − e − Σ t s b / 3 (15)

This is the biggest advancement of the LEACS method; every sample location is analytically optimal for accurately predicting the true collision distance. This is how LEACS aims to limit the magnitude of the slowdown that sampling the finite element solution will invariably cause.

2.9 Windowed multipole formalism

In previous sections and subsections, the notion of sampling a temperature and then returning an accurately Doppler broadened cross-section has been introduced without detailed discussion of whether this is feasible. The use of windowed multipole formalism allows for broadening of cross-sections on the fly, with a much higher accuracy than linear interpolation between single-temperature cross-section libraries could provide (Ducru et al., 2021).

The windowed multipole equation for a cross-section is Equation 16.

σ T z ≈ 1 z 2 R π ∑ j ≥ 1 r j i β ⋅ w z − p j β + ∑ n ≥ − 2 a n D β n z (16)

where z = E, β = k b T 4 A, r j is the jth residue, and p j is the jth pole. The second term on the right-hand side of Equation 16 is the polynomial fit describing broadening outside of the window.

In addition to being more accurate than interpolation, recent advances have made the broadening process extremely efficient by way of evaluating w, the Faddeeva function, as a rational fraction (Forget et al., 2022).

3 Methods

This section contains the detailed descriptions of the two defining algorithms behind the LEACS method. First is the collision sampling algorithm that occurs when sampling a distance to collision, and second is the finite element sampling algorithm invoked when retrieving data from the finite element solution.

3.1 Collision sampling algorithm

Algorithm 2 is called every time a particle intends to move inside a material flagged for the LEACS method. As mentioned previously, the current implementation only considers temperature coupling, although the inclusion of density would follow largely the same logic.

Algorithm 2 LEACS Collision distance sampling.

Note that temperatures are also fit with piecewise quadratic functions. This is done so the particle sqrtkT value can be updated at the end of the process. This is not strictly necessary for collision distance estimation, but other functions in the solver rely on this value. The computational cost of fitting Σ t and T as opposed to just Σ t is negligible.

Additionally, note that because cross-sections now vary spatially in the material, track-length tallies of all types require an analogous quadratic fit of the cross-section relevant to the requested score type. For this work, all tallies, eigenvalue estimation, and normalization are based on collision estimation, which is not impacted by this process.

3.2 Finite element sampling algorithm

Algorithm 3 contains the logic behind the finite element solution sampling. This process operates on a vector of N requested sample positions.

Algorithm 3 Finite-element sampling at a point.

Note that step 1 of the algorithm indicates that element locating is accelerated by a kd-tree. This tree is generated using Nanoflann Blanco and Rai, (2014). Not included in this algorithm description are some of the more mundane but necessary inclusions, like initial instantiation of the finite element solution into memory, generation of the kd-tree, and safety checks to handle edge cases.

4 Analytical error estimation

For an initial assessment of the accuracy of the LEACS method in a controlled environment, two pseudo-pin test cases were generated. Each case is meant to test the upper and lower limits of potential neutron flight distance in the average fuel pin. One case uses a distance of 0.95 cm to emulate a purely radial crossing, while the other uses a distance of 364.85 cm to emulate a purely axial crossing.

Accurate temperature representation is most consequential for predictions of resonance impacts on particle transport. The most significant resonance effect in a typical light water reactor concerns resonance absorption in U-238. The three lowest energy resonances thermally broaden, absorbing more neutrons before they can reach thermal energies. This behavior is one of the strongest temperature-based reactivity feedback effects in the core. To magnify this effect, the material in both cases is set to uniform-density U238.

One-dimensional thermal simulations were conducted in MOOSE Harbour et al. (2025); Shemon et al. (2023) for each case in order to generate a temperature profile to sample from. The temperature and macroscopic cross-section profile versus position are shown in Figure 1. These cross-section plots show the different temperature representations tested in this section: constant average temperature, equidistant sampling, and Lagrange error-adaptive sampling (labeled as LEACS). Note that the LEACS spacing concentrates sample points near the particle’s initial position, as this is where collisions are most likely when the cross-section is high.

The most informative measure of the quality of the quadratic fit is obtained by inspecting the CDF directly. One can generate a CDF and compare it against the control (the same 1,001 equidistant sampled points as before) to identify where collisions are over- or underpredicted. This comparison of CDF errors can be seen in Figure 2 for a constant temperature, 5-sample equidistant spacing, and 5-sample LEACS spacing. Positive errors indicate an overprediction in probability of collision, negative errors indicate an underprediction in probability of collision.

As a final quantification of error from the sampling method, for each description of the spatial variation of the cross-section along the path s, one can calculate the expected distance to collision based on the PDF f ( s ), seen in Equation 17.

s ̄ = ∫ 0 s b s f s d s = s b e − τ s b + ∫ 0 s b s Σ t s e − τ s d s (17)

This value can be used to evaluate the accuracy of each method and parameter combination (number of points sampled, quadrature choice for quadratic fitting) by assigning one finely sampled case as the control. For this work, the control an equidistantly sampled, piecewise quadratic fit with 1,001 sampled points. This error is defined in Equation 18 as:

ϵ ̄ i = s ̄ i − s ̄ 0 (18)

where s ̄ 0 is the expected distance of the control method and s ̄ i is the expected distance of method i. This error has a physical interpretation: it represents the mean discrepancy in distance traveled by a neutron, in centimeters. By evaluating this quantity across a range of neutron energies, the average collision distance error can be compared across methods. This is done in Figure 3.

Note that for radial, or shorter, flight distances, the error between methods is less pronounced. It is at longer flight distances where the error improvement from LEACS is most apparent. It is also unclear as to whether this error comparison reflects what one would see in a real-life pin cell case, as the true neutron spectrum in a reactor is not confined to the vicinity of a single resonance.

To assess the average collision distance errors under a more realistic neutron energy spectrum, an OpenMC 2D pincell model was simulated with 3.2% enriched uranium dioxide fuel and the detailed flux spectrum in the geometry was tallied. This flux spectrum was then normalized and converted to a probability distribution function, and subesquently integrated to obtain a CDF. This permitted sampling of neutron energies from a spectrum more representative of a PWR pincell, enabling estimation of average collision distance error for each method across varying numbers of sample counts, seen in Figure 4. Note that the stochastic uncertainty is included in the error bars.

This shows that the LEACS method can significantly reduce errors in neutron flight distance within materials exhibiting large temperature gradients. The gains from the Lagrange spacing method compared to equidistant spacing increase as the number of points sampled increases, though these gains are most pronounced for large total flight distances.

5 OpenMC implementation

The previous section successfully shows that the LEACS method improves prediction of collision distance in thermally-varying materials, but there are two quantifications that cannot be gained from this analysis. First, error in average distance to collision is useful in that it has a physical meaning for a single neutron being sampled, but most reactor physics simulations are performed to calculate the system eigenvalue or to obtain specific tally quantities. The impact of LEACS on spatial neutron capture profiles constitutes a more practically relevant metric. Secondly, the computational cost of integrating LEACS into the Monte Carlo transport process has yet to be quantified, and is of high importance for evaluating the method’s utility.

By implementing the LEACS method into OpenMC, both of these outstanding questions can be addressed.

5.1 Problem descriptions

Similarly to the reasoning behind the choice of U238 as a material in Section 4, the geometry and materials of the test cases for this section are chosen to both reflect typical PWR conditions and to emphasize the impact of Doppler broadening. There are four test cases studied, two being pin cell models and two being assembly models.

The pin cell geometry is based on AP1000 specifications (Westinghouse, 2012), with generic water moderator, Zircaloy-4 cladding (from Detwiler et al., 2021), and 1% mixed-oxide fuel. This fuel definition originates from Mosteller (Mosteller, 2007) who recorded the largest Doppler impact for this fuel composition, hence why it was chosen for this test case. This is the main test case for determining spatial impact to reaction rates and convergence per number of samples.

The assembly geometry and materials are based on the C5G7 benchmark problem (Cathalau et al., 1996; Lewis et al., 2001), featuring a mix of UO 2 and mixed-oxide fuel pins, along with water reflector regions on the peripheries. This assembly test case is included to aid in quantifying computational cost of the LEACS method in a larger, non-infinite lattice problem. The computational penalty for LEACS will be relatively smaller when solving a problem that already has significant computational penalties, like a large DAGMC assembly. Additionally, because neutron energy determines whether the LEACS method is invoked and only fuel regions use the LEACS method, the differing neutron spectrum in the reflector regions will affect the overall computational cost.

Each of the cases, the pin cell and assembly, has a 2D and 3D version. The 2D case is used when analyzing spatial variation of reaction rates in the radial direction. The 3D case is used when analyzing that same spatial variation in the axial direction. It also provides additional data for assessing computational cost, similar to the inclusion of the C5G7 assembly model in this work.

The thermal finite element solution for each test case is generated in MOOSE using the same material conductivities, and a volumetric heating rate derived from the peak AP1000 linear power rating (Westinghouse, 2012). For the 3D test cases, a cosine axial shape is applied to the volumetric heating rate. The intent is to create a case such that the resulting temperature profile remains consistent with the radial pin temperature predictions from Yadigaroglu (2005), and is guaranteed to have a strong effect on neutron transport so the work can analyze how the LEACS method performs in capturing these effects. This temperature profile is the same for every pin cell in the assembly case, and all 2D cases use the radial temperature profile at the mid-plane, where axial power is greatest. The geometry and temperature profiles of the pin cell can be seen in Figure 5.

5.2 Spatial tally results

As discussed in Section 5.1, the 2D pin cell test case is meant to show how the LEACS method resolves the true radial reaction rate shape in light of the fuel pin temperature gradient (Figure 6). illustrates this effect, with a clear radial variation in the (n,γ) reaction rate (Figure 7). compares the radial profile between a variety of different number of sampled LEACS runs. The root mean square difference is below 0.35% for all cases when compared to the 21-sample case, and the average z-score of the tally differences is below 1.0, suggesting that even three sample points are likely sufficient to capture the true spatial variation in reaction rate.

For the 3D pin cell model (Figure 8), shows the axial comparison between the constant temperature case and the 5-sample LEACS case. As in the 2D radial case, there is a clear impact on the axial (n,γ) reaction rate profile when using LEACS to account for the spatial temperature profile (Figure 9). shows the comparison between the 3-sample and 5-sample LEACS case. There is no clear difference beyond the stochastic noise, suggesting the increased number of sampled points change little in the resolution of the solution.

Although this is not an exhaustive exploration of this topic, the current results indicate that capturing the spatial neutronic effect of the true temperature profile is potentially achievable with a minimal number of samples.

5.3 Computational cost

The final analysis of the LEACS implementation in this work deals with the computational cost of the method. This is perhaps the most critical component of assessing the viability of a method like this, since much of this work has focused on the Lagrange spacing technique to optimize each finite element sample.

The fundamental breakdown of each component of the LEACS method’s time cost can be seen in Table 1, in which each phase corresponds to one of the algorithms previously. As a reference for the magnitude of each time cost, the default collision distance sampling only takes an average of 0.068 microseconds per call. The LEACS timing values are significantly higher, with the point locating and inverse mapping taking the majority of the time. The next highest cost comes from cross-section calculation.

TABLE 1. Timing breakdown by phase.

FE sampling data is from 3.65M calls, see Algorithm 3; Collision distance sampling is from 730K calls, see Algorithm 2. This case was using N = 5 samples along the flight path.

A more comprehensive display of computational cost can be seen in Table 2. Note that the slowdown is proportionally smaller for the larger assembly cases. This is likely due to the previously mentioned effects. A combination of the base calculation being more expensive which reduces the relative penalty of LECAS, and the neutrons spending a greater fraction of their histories in the regions that do not invoke the LEACS method. A slowdown factor exceeding 16 would likely render the method impractical for production calculations, whereas a factor of 1.3 may be acceptable in many applications.

TABLE 2. Simulation slowdown comparison between default and LEACS methods.

The default method lists the simulation time per particle, while slowdown factors are the multiplicative increase to computation time that comes with using the specific LEACS implementation. All simulations were performed with a single thread on a 3.0 GHz Intel i7-1185G7 processor.

6 Conclusion

This work introduces a novel method of accounting for spatially dependent temperatures during Monte Carlo neutron transport by dynamically evaluating the coupled finite element solution. This LEACS method leverages optimally efficient sample spacing that maximizes the utility of each sample point in reconstructing the true collision probability distribution. The implementation of this method in OpenMC demonstrated a clear capacity to represent the thermal impact on neutron transport with a limited number of samples, but also showed that current computational overhead will likely need to be reduced in future work.

This should be achievable through more efficient element search algorithms and inverse mapping processes during the finite element sampling process, although the method will never be as fast as traditional collision distance sampling. A further planned improvement is the incorporation of density sampling, an important capability for use with molten salt reactors and boiling water reactors.

This author foresees that the best use for the LEACS method would be its integration into the CARDINAL code suite, given that CARDINAL already provides a framework for coupling OpenMC with finite element solvers. A critical future analysis would be a direct comparison of LEACS against the adaptive binning approach in terms of both neutronic accuracy and computational cost.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI was used for the initial code development during OpenMC implementation, for grammar and spelling checks of the manuscript, and for suggested edits for brevity.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Mark D. DeHart, Abilene Christian University, United States

Reviewed by: Wen Yin, Harbin Engineering University, China

Hicham Satti, Abdelmalek Essaâdi University, Morocco