Monte Carlo Chord Length Sampling for d-dimensional Markov binary mixtures

# Monte Carlo Chord Length Sampling for d-dimensional Markov binary mixtures

Coline Larmier Den-Service d’Etudes des Réacteurs et de Mathématiques Appliquées (SERMA), CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, FRANCE. Adam Lam School of Nuclear Science & Engineering, Oregon State University, Corvallis, Oregon, USA. Patrick Brantley Lawrence Livermore National Laboratory, P.O. Box 808, Livermore, CA 94551. Fausto Malvagi Todd Palmer Andrea Zoia
###### Abstract

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

###### keywords:
Chord Length Sampling, Markov geometries, benchmark, Monte Carlo, Tripoli-4 ®, Mercury
journal: Journal of Quantitative Spectroscopy & Radiative Transfer

## 1 Introduction

Several applications in nuclear science and engineering involve linear particle transport theory in stochastic media. Examples include neutron diffusion in pebble-bed reactors or randomly mixed water-vapor phases in boiling water reactors Pomraning (1991); Larsen and Vasques (2011); Levermore et al. (1986); Sanchez (1986); Levermore et al. (1988), and inertial confinement fusion Zimmerman (1990); Zimmerman and Adams (1991); Haran et al. (2000). Particle propagation in random media emerges more broadly in material and life sciences and in radiative transport Torquato (2002); Barthelemy et al. (2009); Davis and Marshak (2004); Kostinski and Shaw (2001); Malvagi et al. (1992); Tuchin (2007); Brantley et al. (2017b). Assuming that particles undergo single-speed transport with isotropic scattering, the angular particle flux for each physical realization of the system obeys the linear Boltzmann equation

 ω⋅∇φ+Σ(r)φ=Σs(r)Ωd∫φ(r,ω′)dω′+S. (1)

Here and denote the position and direction variables, respectively, being the total cross section and the source term. The quantity is the surface area of the unit sphere in dimension , being the Gamma function. The quantities , and are in principle random variables, since the materials composing the traversed medium are assumed to be possibly distributed according to some statistical law. The physical observable of interest is typically the ensemble-averaged angular particle flux , or more generally some ensemble-averaged functional of the particle flux, namely,

 ⟨φ(r,ω)⟩=∫P(q)φ(q)(r,ω)dq, (2)

where is the solution of the Boltzmann equation (1) corresponding to a single realization , and is the stationary probability of observing the state for the functions , and  Pomraning (1991); Zuchuat et al. (1994).

Exact solutions for can be in principle obtained in the following way: first, a realization of the medium is sampled from the underlying mixing statistics; then, the linear transport equation (1) corresponding to this realization is solved by either deterministic or Monte Carlo methods, and the physical observables of interest are determined; a sufficiently large collection of realizations is produced; and ensemble averages are finally taken for the physical observables.

Reference solutions are very demanding in terms of computational resources, especially if transport is to be solved by Monte Carlo methods in order to preserve the highest possible accuracy in solving the Boltzmann equation. In principle, it would be thus desirable to directly derive a single equation for the ensemble-averaged flux . A widely adopted model of random media is the so-called binary stochastic mixing, where only two immiscible materials (say and ) are present Pomraning (1991). Then, by averaging Eq. (1) over realizations having material at , we obtain the following equation for

 [ω⋅∇+Σα]pα⟨φα⟩=pαΣs,αΩd∫⟨φα(r,ω′)⟩dω′ +pβ,α⟨φβ,α⟩−pα,β⟨φα,β⟩+pαSα (3)

where is the probability of finding the material of index at position . Here represents the probability per unit length of crossing the interface from material to material for a particle located at and travelling in direction . The quantity denotes the angular flux averaged over those realizations where there is a transition from material to material for a particle located at and travelling in direction . The cross sections and are those of material . The equation for is immediately obtained from Eq. (3) by permuting the indices and . Excluding the special case of particle transport in the absence of scattering, we are thus led to an infinite hierarchy for in Eqs. (3).

In order to explicitly derive the ensemble-averaged flux , it is therefore necessary to introduce a closure formula, which will in general depend on the underlying mixing statistics Pomraning (1991); Zuchuat et al. (1994); Su and Pomraning (1995). The celebrated Levermore-Pomraning model assumes for instance for homogeneous Markov mixing statistics, with

 pi,j(r,ω)=piΛi(ω), (4)

where is the mean chord length for trajectories crossing material in direction  Pomraning (1991). Several generalisations of this model have been later proposed, including higher-order closure schemes Pomraning (1991); Su and Pomraning (1995).

In parallel, a family of Monte Carlo algorithms have been conceived in order to approximate the ensemble-averaged solutions to various degrees of accuracy Zimmerman and Adams (1991); Donovan et al. (2003); Donovan and Danon (2003). Their common feature is that they allow a simpler treatment of transport in stochastic mixtures (typically by neglecting the correlations on particle trajectories induced by the spatial disorder), which might be convenient in practical applications. In this context, a prominent role is played by the so-called Chord Length Sampling (CLS) algorithm, which is supposed to solve the Levermore-Pomraning model for Markovian binary mixing Zimmerman and Adams (1991); Sahni (1989a, b). The basic idea behind CLS is that the interfaces between the constituents of the stochastic medium are sampled on-the-fly during the particle displacements by drawing the distances to the following material boundaries from a distribution depending on the mixing statistics. The free parameters of the CLS model are the average chord length through each material and the volume fraction . Since the spatial configuration seen by each particle is regenerated at each particle flight, the CLS corresponds to an annealed disorder model, as opposed to the quenched disorder of the reference solutions, where the spatial configuration is frozen for all the traversing particles. Generalization of these Monte Carlo algorithms including partial memory effects due to correlations for particles crossing back and forth the same materials have been also proposed Zimmerman and Adams (1991).

In order to quantify the accuracy of the various approximate models, comparisons with respect to reference solutions are mandatory. For instance, although originally formulated for Markov statistics, CLS has been extensively applied also to randomly dispersed spherical inclusions into background matrices, with application to pebble-bed and very high temperature gas-cooled reactors Donovan et al. (2003); Donovan and Danon (2003), and several benchmark problems have been examined in two and three dimensions Donovan et al. (2003); Donovan and Danon (2003); Brantley and Martos (2011); Brantley (2014). Some methods to mitigate the errors between CLS and the reference solutions have been presented in the context of eigenvalue calculations, e.g., in Ji and Brown (2013). For Markov mixing, a number of benchmark problems comparing CLS and reference solutions have been proposed in the literature so far Adams et al. (1989); Brantley (2011); Zuchuat et al. (1994); Brantley and Palmer (2009); Brantley (2009), with focus exclusively on geometries, either of the rod or slab type. Flat two-dimensional configurations have received less attention Haran et al. (2000).

In a series of recent papers, some of the authors have provided reference solutions for particle transport in extruded two-dimensional and full three-dimensional random media with Markov statistics Larmier et al. (2017a, b), where the spatial disorder has been generated by means of homogeneous and isotropic -dimensional Poisson tessellations Larmier et al. (2016). In this work, we will compare the CLS simulation results to the reference solutions for the classical benchmark problem proposed by Adams, Larsen and Pomraning for transport in stochastic media Adams et al. (1989) and revisited by Brantley Brantley (2011). The case of slab disorder has been considered previously in the literature Adams et al. (1989); Brantley (2011); Zuchuat et al. (1994); Brantley and Palmer (2009); Brantley (2009) and will be reported here for the sake of completeness. In addition, we will also consider extruded and full Markov mixing configurations. The physical observables of interest will be the particle flux , the transmission coefficient and the reflection coefficient : we will examine the discrepancies between reference and CLS simulation results as a function of the benchmark configurations and of the system dimensionality . In order to verify the consistency of the proposed results, the CLS calculations will be performed by using two independent Monte Carlo implementations of the CLS algorithm, in the Tripoli-4 ® code Brun et al. (2015) and in the Mercury code Brantley et al. (2017a, c), respectively.

This paper is organized as follows: in Sec. 2 we will recall the benchmark specifications that will be used in this work. In Secs. 3 and 4 we will detail the methods and the algorithms that we have adopted in order to produce reference and CLS results, respectively. Simulation findings will be illustrated and discussed in Sec. 5. Conclusions will be finally drawn in Sec. 6.

## 2 Benchmark specifications

In order for the paper to be self-contained, we start by recalling the benchmark specifications that have been selected for this work, which are essentially taken from those originally proposed in Adams et al. (1989) and Zuchuat et al. (1994), and later extended in Brantley (2011); Brantley and Palmer (2009); Brantley (2009).

We consider single-speed linear particle transport through a stochastic binary medium with homogeneous Markov mixing. The medium is non-multiplying, with isotropic scattering. The geometry consists of a cubic box of side (in arbitrary units), with reflective boundary conditions on all sides of the box except two opposite faces (say those perpendicular to the axis), where leakage boundary conditions are imposed 111In Adams et al. (1989) and Zuchuat et al. (1994), system sizes and were also considered, but in this work we will focus on the case , which leads to more physically relevant configurations.. Two kinds of non-stochastic sources will be considered: either an imposed normalized incident angular flux on the leakage surface at (with zero interior sources), or a distributed homogeneous and isotropic normalized interior source (with zero incident angular flux on the leakage surfaces). Following the notation in Brantley (2011), the benchmark configurations pertaining to the former kind of source will be called suite I, whereas those pertaining to the latter will be called suite II. The material properties for the Markov mixing are entirely defined by assigning the average chord length for each material , namely , which in turn allows deriving the homogeneous probability of finding material at an arbitrary location within the box, namely

 pi=ΛiΛi+Λj. (5)

By definition, the material probability yields the volume fraction for material . The cross sections for each material will be denoted as customary for the total cross section and for the scattering cross section. The average number of particles surviving a collision in material will be denoted by . The physical parameters for the benchmark configurations are recalled in Tabs. 1 and 2: the benchmark specifications include three cases (numbered , and , corresponding to different materials), and three sub-cases (noted , and , corresponding to different for a given material) for each case Adams et al. (1989). The so-called atomic mix limit Pomraning (1991), where one assumes that the statistical disorder can be approximated by simply taking a full homogenization of the physical properties based on the ensemble-averaged cross sections, has been examined, e.g., in Brantley (2011) for and in Larmier et al. (2017a) for and and will not be considered here.

The physical observables of interest for the proposed benchmark will be the ensemble-averaged outgoing particle currents on the two surfaces with leakage boundary conditions, the ensemble-averaged scalar particle flux along , and the total scalar flux . For the suite I configurations, the outgoing particle current on the side opposite to the imposed current source will represent the ensemble-averaged transmission coefficient, namely, , whereas the outgoing particle current on the side of the current source will represent the ensemble-averaged reflection coefficient, namely, . For the suite II configurations, the outgoing currents on opposite faces are expected to be equal (within statistical fluctuations), for symmetry reasons. In this case, we also introduce the average leakage current .

## 3 Reference solutions

For particle transport in the presence of quenched disorder with -dimensional Markov mixing, reference solutions for the ensemble-averaged scalar particle flux and the currents and have been first obtained and thoroughly described in Larmier et al. (2017a). For this work, CEA has produced a new set of reference solutions with lower statistical uncertainty 222In order to reduce computer time for highly fragmented geometries including hundreds of thousands of polyhedra, in Larmier et al. (2017a) we used thin empty boxes adjacent to the two free surfaces in order to sample the source particles for suite I configurations and to record the weights of the particles escaping from the viable domain. In the new simulations, we have also decreased the thickness of these boxes in order to improve the overall accuracy.. Here we will briefly detail the methods and the simulation parameters that have been used and recall the changes introduced with respect to Larmier et al. (2017a).

### 3.1 Poisson tessellations

Random tessellations are stochastic aggregates of disjoint and space-filling cells obeying a given distribution Santalo (1976). Poisson tessellations are obtained by partitioning a domain of a -dimensional space by sampling -dimensional hyper-planes from an auxiliary Poisson process Santalo (1976); Miles (1964, 1972). An explicit construction amenable to Monte Carlo realizations for two-dimensional homogeneous and isotropic Poisson geometries of finite size has been established in Switzer (1964). A generalization of this algorithm to -dimensional domains has been recently proposed Ambos and Mikhailov (2011). The construction of Poisson stochastic geometries depends on a single free parameter , which takes the name of tessellation density, and is such that an arbitrary segment of length will have on average intersections with the random hyper-planes.

The algorithm for the slab tessellations is recalled in Adams et al. (1989), based on the Poisson process on the line. For the extruded tessellations, we begin by creating an isotropic Poisson tessellation of a square of side , according to the algorithm detailed in Lepage at al. (2011). The full geometrical description for the cube is simply achieved by extruding the random polygons of the plane along the orthogonal (say ) axis. The algorithm for tessellations of a cube of side by drawing random planes has been detailed in Larmier et al. (2016).

Isotropic Poisson geometries satisfy a Markov property: for domains of infinite size, arbitrary drawn lines will be cut by the -surfaces of the -polyhedra into segments whose lengths are exponentially distributed, with average chord length  Santalo (1976). The quantity intuitively defines the correlation length of the Poisson geometry, i.e, the typical linear size of a volume composing the random tessellation.

### 3.2 Colored stochastic geometries

Binary Markov mixtures required for the benchmark specifications are obtained as follows: first, a -dimensional Poisson tessellation is constructed as described above. Then, each polyhedron of the geometry is assigned a material composition by formally attributing a distinct ‘color’, say or , with associated complementary probabilities and  Pomraning (1991). This gives rise to (generally) non-convex and clusters, each composed of a random number of convex polyhedra. It can be shown that the average chord length through clusters with composition is related to the correlation length of the geometry via , and for we similarly have . This yields , and we recover

 pα=ΛΛβ=ΛαΛα+Λβ. (6)

Based on the formulas above, and using , the parameters of the colored Poisson geometries corresponding to the benchmark specifications provided in Tab. 1 are easily derived.

### 3.3 Particle transport and ensemble averages

For each benchmark case and sub-case, a large number of geometries has been generated, and the material properties have been attributed to each volume as described above. Then, for each realization of the ensemble, linear particle transport has been simulated by using the production Monte Carlo code Tripoli-4 ®, developed at CEA Brun et al. (2015). Tripoli-4 ® is a general-purpose stochastic transport code capable of simulating the propagation of neutral and charged particles with continuous-energy cross sections in arbitrary geometries. In order to comply with the benchmark specifications, constant cross sections adapted to mono-energetic transport and isotropic angular distributions have been prepared. The number of simulated particle histories per configuration is . For a given physical observable , the benchmark solution is obtained as the ensemble average

 ⟨O⟩=1MM∑k=1Ok, (7)

where is the Monte Carlo estimate for the observable obtained for the -th realization. Specifically, currents and at a given surface are estimated by summing the statistical weights of the particles crossing that surface. Scalar fluxes have been tallied using the standard track length estimator over a pre-defined spatial grid containing uniformly spaced meshes along the axis.

The error affecting the average observable results from two separate contributions, the dispersion

 σ2G=1MM∑k=1Ok2−⟨O⟩2 (8)

of the observables exclusively due to the stochastic nature of the geometries and of the material compositions, and

 σ2O=1MM∑k=1σ2Ok, (9)

which is an estimate of the variance due to the stochastic nature of the Monte Carlo method for particle transport, being the dispersion of a single calculation Donovan and Danon (2003); Donovan et al. (2003). The statistical error on is then estimated as

 σ[⟨O⟩]=√σ2GM+σ2O. (10)

In order to reduce the dispersion of the observables due to the statistical nature of the geometries, for the new set of reference solutions computed for this work we have increased the number of realizations for the benchmark configurations displaying larger correlation lengths (i.e., larger material chunks).

For slab tessellations, we have taken for sub-cases ; for sub-cases ; and for sub-cases . Otherwise, we have used the same number of realizations as in Larmier et al. (2017a), namely, for sub-case of the suite II, and for the remaining cases and sub-cases.

For the extruded tessellations, we have taken for the sub-cases and for the sub-cases . Otherwise, we have used the same number of realizations as in Larmier et al. (2017a), namely, .

Finally, for the tessellations we have taken for the sub-case of the suite II; for all the other sub-cases of case ; for the sub-case of the suite II; and for all the other sub-cases of case . For all remaining cases and sub-cases, we have used the same number of realizations as in Larmier et al. (2017a), namely, .

Actually, increasing the dimension implies a better statistical mixing (in other words, a single realization is more representative of the average behaviour), at the expense of increasing the computational burden (each realization takes longer both for generation and for Monte Carlo transport).

Transport calculations have been run on a cluster based at CEA, with Intel Xeon E5-2680 V2 2.8 GHz processors. The average computer time globally increases as a function of dimension, but depends also on the correlation lengths, volume fractions, and material properties such as cross sections and scattering probabilities. For the simulations discussed here we have largely benefited from a feature implemented in the code Tripoli-4 ®, namely the possibility of reading pre-computed connectivity maps for the volumes composing the geometry. During the generation of the Poisson tessellations, care has been taken so as to store the indices of the neighbouring volumes for each realization, which means that during the geometrical tracking a particle will have to find the following crossed volume in a list that might be considerably smaller than the total number of random volumes composing the box (depending on the features of the random geometry).

## 4 The Chord Length Sampling approach

Reference solutions based on the quenched disorder approach are computationally expensive, so that intensive research efforts have been devoted to the development of Monte Carlo-based annealed disorder models capable of approximating the ensemble observables on-the-fly, i.e., with a single particle transport simulation. The pioneering work by Zimmerman and Adams Zimmerman (1990); Zimmerman and Adams (1991) has led to a family of algorithms that go now under the name of Chord Length Sampling methods. In particular, it has been shown that the standard form of the CLS (Algorithm A in Zimmerman and Adams (1991)) formally solves the Levermore-Pomraning equations, i.e., Eq. (3) with the closure formula (4), corresponding to Markov mixing with the approximation that memory of the crossed material interfaces is lost at each particle flight Sahni (1989a, b).

Algorithm A proceeds as follows Zimmerman and Adams (1991): each particle history begins by sampling position, angle and velocity from the specified source, as customary. Moreover, the particle is assigned a supplementary attribute, the material label, which is sampled from the probability . Then we need to compute three distances, denoted respectively , , and . The quantity is the distance to the next physical boundary, along the current direction of the particle. The quantity is the distance to the next collision, which is determined by using the material cross section that has been chosen at the previous step: if the particle has a material , e.g., then will be drawn from an exponential distribution of parameter . Finally, the quantity is the distance to the next material interface, which is sampled from an exponential distribution with parameter , i.e., the average chord length of material , if the particle has a material label (whence the name of CLS).

Then, the minimum distance among and must be selected: if the minimum is , the particle is moved along a straight line until it hits the external boundary; if the minimum is , the particle is moved to the collision point, and the outgoing particle features are selected according to the collision kernel pertaining to the current material label. If the minimum is , the particle is moved to the interface between the two materials, and the material label is switched. If the particle is not absorbed, a new set of distances and are determined. During the time spent within the random medium, the particle will be thus either colliding within a random chunk, or crossing the interface between two chunks; the particle will ultimately get absorbed in one of the chunks, or escape out of the boundaries of the random medium. The Monte Carlo estimators for the scalar flux and the currents are the same as those for the reference solutions described above.

As observed above, Algorithm A assumes that the particle has no memory of its past history, and in particular the crossed interfaces are immediately forgotten (which is coherent with the closure formula of the Levermore-Pomraning model). In this respect, CLS is an approximation of the exact treatment of disorder-induced spatial correlations (actually, it can be shown that CLS is exact only for pure absorbers). As a result, Algorithm A is expected to be less accurate in the presence of strong scatterers with optically thick mean material chunk length. A thorough discussion of the shortcomings of the CLS approach for can be found, e.g., in Liang and Ji (2011).

### 4.1 Slab geometries

For mono-energetic particle transport in slab geometries with isotropic scattering, the Boltzmann equation (1) yields

 μ∂∂xφ+Σ(x)φ=Σs(x)2∫1−1dμ′φ(x,μ′), (11)

where is the angular particle flux for particles at position with a direction cosine with respect to the axis. The source and the boundary conditions depend on the benchmark specifications.

Correspondingly, the CLS algorithm that formally solves the Levermore-Pomraning model as applied to Eq. 11 is the following. For suite I, the source particle position is set to , and the direction cosine is sampled from a cosine distribution, namely,

 μ=√ξ, (12)

where is a uniform random number in , in order to ensure the isotropic incident flux condition. For suite II, the starting position is sampled uniformly in , and the direction cosine is sampled uniformly in in order to ensure the uniform and isotropic source condition. According to the Levermore-Pomraning prescription, the distance to material interfaces for a particle in material is sampled from an exponential distribution as follows:

 di=−Λα|μ|ln(1−ξ), (13)

where the factor accounts for the projection of the distance along the axis. The distance to the next collision is sampled from the exponential distribution of parameter , and the distance to the boundary is computed as customary. For isotropic scattering, the cosine direction after collision is sampled uniformly in .

### 4.2 Two-dimensional extruded geometries

Assuming again mono-energetic particle transport with isotropic scattering, the Boltzmann equation for two-dimensional geometries extruded in the axis direction yields

 √1−μ2cos(ϕ)∂∂xφ+√1−μ2sin(ϕ)∂∂yφ= Σ(x,y)φ+Σs(x,y)4π∫1−1dμ′∫2π0dϕ′φ(x,y,μ′,ϕ′), (14)

where is the angular particle flux for particles being at position with a direction cosine with respect to the axis and a polar angle with respect to the axis.

The CLS algorithm that formally corresponds to solving the Levermore-Pomraning model as applied to Eq. 14 is the following. For suite I, the source particle positions are set to and taken uniformly in . Then we sample a direction cosine (with respect to the axis) from

 μ′=√ξ (15)

where is taken in , and a polar angle (with respect to the axis) uniform in . The initial particle direction is

 ω0={μ′Q,√1−μ′2cos(ϕ′)Q}, (16)

with

 Q=√μ′2+(1−μ′2)cos2(ϕ′), (17)

in order to ensure the isotropic incident flux condition, and the initial direction cosine is defined by

 μ0=√1−μ′2sin(ϕ′). (18)

For suite II, the starting positions are sampled uniformly in , the direction cosine is sampled uniformly in and the polar angle is sampled uniformly in in order to ensure the uniform and isotropic source condition, which yields the initial particle direction

 ω0={cos(ϕ),sin(ϕ)}. (19)

According to the Levermore-Pomraning prescription, the distance to material interfaces for a particle in material is sampled from an exponential distribution as follows:

 di=−Λα√1−μ2ln(1−ξ), (20)

where the factor again accounts for the projection of the distance on the plane. The distance to the next collision is sampled from the exponential distribution of parameter , and the distance to the boundary is computed as customary. For isotropic scattering, the cosine direction after collision is sampled uniformly in , and the polar angle is sampled uniformly in ; the particle direction is then given by

 ω={cos(ϕ),sin(ϕ)}. (21)

### 4.3 Three-dimensional geometries

The Boltzmann equation for mono-energetic transport with isotropic scattering in three-dimensional geometries yields

 (√1−μ2cos(ϕ)∂∂x+√1−μ2sin(ϕ)∂∂y+μ∂∂z)φ= Σ(x,y,z)φ+Σs(x,y,z)4π∫1−1dμ′∫2π0dϕ′φ(x,y,z,μ′,ϕ′), (22)

where is the angular particle flux for particles being at position with a direction cosine with respect to the axis and a polar angle with respect to the axis.

The CLS algorithm that formally corresponds to solving the Levermore-Pomraning model as applied to Eq. 22 is the following. For suite I, the source particle positions are set to and taken uniformly in . Then we sample a direction cosine (with respect to the axis) from

 μ′=√ξ (23)

where is taken in , and a polar angle (with repect to the axis) uniform in . The initial particle direction is

 ω0={μ′,√1−μ′2cos(ϕ′),√1−μ′2sin(ϕ′)} (24)

in order to ensure the isotropic incident flux condition. For suite II, the starting positions are sampled uniformly in , the direction cosine is sampled uniformly in and the polar angle is sampled uniformly in in order to ensure the uniform and isotropic source condition, which yields the initial particle direction

 ω0={√1−μ2cos(ϕ),√1−μ2sin(ϕ),μ}. (25)

According to the Levermore-Pomraning prescription, the distance to material interfaces for a particle in material is sampled from an exponential distribution as follows:

 di=−Λαln(1−ξ). (26)

The distance to the next collision is sampled from the exponential distribution of parameter , and the distance to the boundary is computed as customary. For isotropic scattering, the cosine direction after collision is sampled uniformly in , and the polar angle is sampled uniformly in ; the particle direction is then given by

 ω={√1−μ2cos(ϕ),√1−μ2sin(ϕ),μ}. (27)

## 5 Simulation results

The simulation results for the total scalar flux , the transmission coefficient and the reflection coefficient are provided in Tabs. 3 to 5 for the benchmark cases corresponding to suite I, and in Tabs. 6 to 8 for the benchmark cases corresponding to suite II, respectively. The reference solutions have been computed by following the procedure detailed in Sec. 3, based on Larmier et al. (2017a).

The CLS results have been obtained with both Tripoli-4 ® and Mercury Monte Carlo codes by following the procedure described in Sec. 4. We will denote by the resulting statistical uncertainty associated to each physical observable . For the Tripoli-4 ® CLS simulations of the -dimensional benchmark configurations we have used 10 particle histories. Mercury is a Monte Carlo particle transport code being developed at Lawrence Livermore National Laboratory Brantley et al. (2017a, c). The Monte Carlo Levermore-Pomraning CLS algorithm was previously implemented in Mercury Brantley (2014) in a manner consistent with the algorithmic descriptions in Zimmerman and Adams (1991); Brantley (2011) and Sec. 4. The Mercury Levermore-Pomraning implementation has been demonstrated Brantley (2014) to accurately reproduce the independent one-dimensional slab geometry Monte Carlo Levermore-Pomraning results in Brantley (2011). We modelled the three-dimensional benchmark suites I and II using the Mercury Levermore-Pomraning CLS implementation with 10 particle histories. We obtained results that are generally statistically equivalent to the Tripoli-4 ® CLS results to typically within three standard deviations for the reflection and transmission coefficients and the scalar flux distributions (agreement to typically four to five digits). For this paper, we will present only the Tripoli-4 ® simulation results. Computer times for the reference and CLS solutions are also provided in the same tables: not surprisingly, the CLS approach is much faster than the reference method, since a single transport simulation is needed.

As a general remark, the accuracy of CLS with respect to reference solutions increases with increasing system dimensionality . This is expected on physical grounds, since the higher and the smaller is the impact of the spatial correlations: a particle undergoing back-scattering is less likely to cross exactly the same material interface as the one crossed during the previous flight. In other words, the approximations introduced in the CLS algorithm by neglecting spatial correlations will have a weaker effect on particle transport. Nonetheless, simulation results show a few exceptions among the examined configurations. Moreover, the accuracy of CLS also generally improves when increasing the tessellation density, i.e., decreasing the average chord length: configurations pertaining to case globally show a better agreement than those of case , and those of case show a better agreement than those of case .

The effects of system dimensionality on the discrepancies between CLS and exact solutions are stronger for configurations with smaller average chord lengths. This behaviour is again consistent with the fact that increasing the chord length induces larger chunks of materials, and for chunks that span a large fraction of the entire geometry the impact of dimensionality must be rather weak: in this regime, particle transport is mostly influenced by the material volume fractions (i.e., the coloring probability).

The behaviour of suite II configurations is quite similar to that of suite I configurations, and no specific trend due to the source and/or initial conditions can be easily detected.

The spatial scalar flux within the box is illustrated in Figs. 1 to 3 for case to case , respectively. The discrepancies between CLS and reference solutions for this observable have the same behaviour as for the scalar quantities described above. The discrepancy decreases with increasing system dimensionality and with decreasing average chord length. For dense geometries (case ) the effects of dimensionality on the discrepancy are rather strong, and become less appreciable for less dense geometries. The kind of source and/or initial conditions plays again a minor role. This analysis is confirmed by plotting the differences between reference and CLS solutions (see Figs. 4 to 6 for case to case , respectively). Since both reference and CLS solutions are affected by a statistical uncertainty, the error bars on have been computed by taking the combined variance

 σ[Δ[O]]=√σ2[⟨O⟩]+σ2CLS[O] (28)

for each observable .

## 6 Conclusions

The Chord Length Sampling algorithm efficiently provides approximate ensemble-averaged observables corresponding to the Levermore-Pomraning model for Markovian binary mixing. The interfaces between the constituents of the random medium are sampled on-the-fly during the particle displacements by drawing the distances to the following material boundaries from a distribution depending on the mixing statistics: the correlations on particle trajectories induced by the spatial disorder are thus neglected. Comparisons of CLS solutions with respect to reference results are mandatory in order to quantify the degree of approximations introduced in these models. For Markov mixing, a number of benchmark problems have been proposed in the literature for this purpose, but so far analyses have been conducted in one-dimensional media of the rod or slab type.

In this work we have contrasted CLS simulation results to the reference solutions for the classical benchmark problem proposed by Adams, Larsen and Pomraning, and recently revisited by Brantley, for particle propagation in stochastic media with binary Markov mixing. In particular, we have examined the evolution of the particle flux, the transmission coefficient and the reflection coefficient as a function of the benchmark configurations and of the system dimension .

Two main trends have been detected: the accuracy of CLS algorithm with respect to reference solutions generally increases with increasing system dimensionality. Moreover, the accuracy of the CLS algorithm increases for decreasing average chord length, i.e., for denser stochastic tessellations. The impact of dimensionality is particularly relevant for case configurations (which have smaller chord lengths), and progressively diminishes for configurations having larger material chunks. The considerations presented in this paper, although derived strictly speaking for the Adams, Larsen and Pomraning benchmark considered here, seem to be quite general.

This work represents a first step towards extensive comparisons between CLS and reference solutions for Markov mixing statistics in higher dimensions. Furthermore, extension of these comparisons to reference solutions for other types of -dimensional mixing statistics based on spatial tessellations (such as the Poisson-Voronoi model presented in Larmier et al. (2017b)) would be interesting topics for future research.

## Acknowledgements

TRIPOLI-4 ® is a registered trademark of CEA. C. Larmier, A. Zoia and F. Malvagi wish to thank Électricité de France (EDF) for partial financial support and A. Mazzolo and E. Dumonteil for fruitful discussions. Work of P. Brantley performed under the auspices of the U. S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

## References

• Adams et al. (1989) Adams ML, Larsen EW, Pomraning GC. Benchmark results for particle transport in a binary Markov statistical medium. J Quant Spectrosc Radiat Transfer 1989;42:253-66.
• Brantley (2011) Brantley PS. A benchmark comparison of Monte Carlo particle transport algorithms for binary stochastic mixtures. J Quant Spectrosc Radiat Transfer 2011;112:599-618.
• Pomraning (1991) Pomraning GC. Linear kinetic theory and particle transport in stochastic mixtures. River Edge, NJ, USA: World Scientific Publishing; 1991.
• Larsen and Vasques (2011) Larsen EW, Vasques R. A generalized linear Boltzmann equation for non-classical particle transport. J Quant Spectrosc Radiat Transfer 2011:112;619-31.
• Levermore et al. (1986) Levermore CD, Pomraning GC, Sanzo DL, Wong J. Linear transport theory in a random medium. J Math Phys 1986;27:2526-36.
• Sanchez (1986) Sanchez R. Linear kinetic theory in stochastic media. J Math Phys 1988;30:2498-2511.
• Levermore et al. (1988) Levermore CD, Pomraning GC, Wong J. Renewal theory for transport processes in binary statistical mixtures. J Math Phys 1988;29:995-1004.
• Zimmerman (1990) Zimmerman GB. Recent developments in Monte Carlo techniques. Lawrence Livermore National Laboratory Report UCRL-JC-105616; 1990.
• Zimmerman and Adams (1991) Zimmerman GB, Adams ML. Algorithms for Monte Carlo particle transport in binary statistical mixtures. Trans Am Nucl Soc 1991:66;287.
• Haran et al. (2000) Haran O, Shvarts D, Thieberger R. Transport in 2D scattering stochastic media: simulations and models. Phys Rev E 2000:61;6183-89.
• Torquato (2002) Torquato S. Random heterogeneous materials: microstructure and macroscopic properties. New York, USA: Springer-Verlag; 2002.
• Barthelemy et al. (2009) Barthelemy P, Bertolotti J, Wiersma DS. A Lévy flight for life. Nature 2009:453,495-98.
• Davis and Marshak (2004) Davis AB, Marshak A. Photon propagation in heterogeneous optical media with spatial correlations. J Quant Spectrosc Radiat Transfer 2004:84;3-34.
• Kostinski and Shaw (2001) Kostinski AB, Shaw RA. Scale-dependent droplet clustering in turbulent clouds. J Fluid Mech 2001:434;389-98.
• Malvagi et al. (1992) Malvagi F, Byrne RN, Pomraning GC, Somerville RCJ. Stochastic radiative transfer in partially cloudy atmosphere. J Atm Sci 1992:50;2146-58.
• Tuchin (2007) Tuchin V. Tissue optics: light scattering methods and instruments for medical diagnosis. Cardiff, UK: SPIE Press; 2007.
• Brantley et al. (2017b) Brantley PS, Gentile NA, Zimmerman GB. Beyond Levermore-Pomraning for implicit Monte Carlo radiative transfer in binary stochastic media. In: Proceedings of M&C 2017 - International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering, Jeju, Korea. April 16-20, 2017 [on USB].
• Zuchuat et al. (1994) Zuchuat O, Sanchez R, Zmijarevic I, Malvagi F. Transport in renewal statistical media: benchmarking and comparison with models. J Quant Spectrosc Radiat Transfer 1994;51:689-722.
• Su and Pomraning (1995) Su B, Pomraning GC. Modification to a previous higher order model for particle transport in binary stochastic media. J Quant Spectrosc Radiat Transfer 1995;54:779-801.
• Donovan et al. (2003) Donovan TJ, Sutton TM, Danon Y. Implementation of Chord Length Sampling for transport through a binary stochastic mixture. In: Proceedings of the nuclear mathematical and computational sciences: a century in review, a century anew, Gatlinburg, TN. La Grange Park, IL: American Nuclear Society; April 6-11, 2003 [on CD-ROM].
• Donovan and Danon (2003) Donovan TJ, Danon Y. Application of Monte Carlo chord-length sampling algorithms to transport through a two-dimensional binary stochastic mixture. Nucl Sci Eng 2003;143:226-39.
• Sahni (1989a) Sahni DC. Equivalence of generic equation method and the phenomenological model for linear transport problems in a two-state random scattering medium. J Math Phys 1989:30; 1554-9.
• Sahni (1989b) Sahni DC. An application of reactor noise techniques to neutron transport problems in a random medium. Ann Nucl Energy 1989:16;397-408.
• Brantley and Martos (2011) Brantley PS, Martos JN. Impact of spherical inclusion mean chord length and radius distribution on three-dimensional binary stochastic medium particle transport. In: Proceedings of the international conference on mathematics, computational methods & reactor physics (M&C2011), Rio de Janeiro, RJ, Brazil; May 8-12, 2011 [on CD-ROM].
• Brantley (2014) Brantley PS. Benchmark investigation of a 3D Monte Carlo Levermore-Pomraning algorithm for binary stochastic media. Trans Am Nucl Soc, Anaheim, CA; November 9-13 2014.
• Ji and Brown (2013) Liang C, Ji W, Brown FB. Chord Length Sampling method for analyzing stochastic distribution of fuel particles in continuous energy simulations. Ann Nucl Energy 2013:53;140-6.
• Brantley and Palmer (2009) Brantley PS, Palmer TS. Levermore-Pomraning model results for an interior source binary stochastic medium benchmark problem. In: Proceedings of the international conference on mathematics, computational methods & reactor physics (M&C2009), Saratoga Springs, New York. La Grange Park, IL: American Nuclear Society; May 3-7, 2009 [on CD-ROM].
• Brantley (2009) Brantley PS. A comparison of Monte Carlo particle transport algorithms for binary stochastic mixtures. In: Proceedings of the international conference on mathematics, computational methods & reactor physics (M&C2009), Saratoga Springs, New York. La Grange Park, IL: American Nuclear Society; May 3-7, 2009 [on CD-ROM].
• Larmier et al. (2017a) Larmier C, Hugot FX, Malvagi F, Mazzolo A, Zoia A. Benchmark solutions for transport in -dimensional Markov binary mixtures. J Quant Spectrosc Radiat Transfer 2017:189;133–148.
• Larmier et al. (2017b) Larmier C, Zoia A, Malvagi F, Dumonteil E, Mazzolo A. Monte Carlo particle transport in random media: The effects of mixing statistics. J Quant Spectrosc Radiat Transfer 2017:196;270–86.
• Larmier et al. (2016) Larmier C, Dumonteil E, Malvagi F, Mazzolo A, Zoia A. Finite-size effects and percolation properties of Poisson geometries. Phys Rev E 2016:94;012130.
• Brun et al. (2015) Brun E, Damian F, Diop CM, Dumonteil E, Hugot FX, Jouanne C, Lee YK, Malvagi F, Mazzolo A, Petit O, Trama JC, Visonneau T, Zoia A. TRIPOLI-4, CEA, EDF and AREVA reference Monte Carlo code. Ann Nucl Energy 2015:82;151-60.
• Brantley et al. (2017a) Brantley PS, Bleile RC, Dawson SA, McKinley MS, O’Brien MJ, Pozulp MM, Procassini RJ, Richards DF, Sepke SM, Stevens DE, Walsh JA. Mercury User Guide: Version 5.6. Lawrence Livermore National Laboratory Report LLNL-SM-560687 (Modification #12) 2017.
• Brantley et al. (2017c) Brantley PS, Bleile RC, Dawson SA, Gentile NA, McKinley MS, O’Brien MJ, Pozulp MM, Richards DF, Stevens DE, Walsh JA, Childs H. LLNL Monte Carlo Transport Research Efforts for Advanced Computing Architectures. In: Proceedings of M&C 2017 - International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering, Jeju, Korea. April 16-20, 2017 [on USB].
• Santalo (1976) Santaló LA. Integral geometry and geometric probability. Reading, MA, USA: Addison-Wesley; 1976.
• Miles (1964) Miles RE. Random polygons determined by random lines in a plane. Proc Nat Acad Sci USA 1964: 52; 901-7.
• Miles (1972) Miles RE. The random division of space. Suppl Adv Appl Prob 1972:4;243-66.
• Switzer (1964) Switzer P. Random set process in the plane with Markov property. Ann Math Statist 1965:36;1859-63.
• Ambos and Mikhailov (2011) Ambos AYu, Mikhailov GA. Statistical simulation of an exponentially correlated many-dimensional random field. Russ J Numer Anal Math Modelling 2011:26;263-73.
• Lepage at al. (2011) Lepage T, Delaby L, Malvagi F, Mazzolo A. Monte Carlo simulation of fully Markovian stochastic geometries. Prog Nucl Sci Techn 2011:2;743-48.
• Liang and Ji (2011) Liang C, W Ji. On the Chord Length Sampling in 1-d binary stochastic media. Transp. Theory Stat Phys 2011:40;282-303.