Nonrecursive residual Monte Carlo method for S_N transport discretization error estimation

Nonrecursive residual Monte Carlo method for S N transport discretization error estimation 1 corresp first-author yunhuangz@zju.edu.cn yunhuangz@zju.edu.cn yunhuangz@zju.edu.cn https://orcid.org/0000-0002-3391-6424 Department of Physics, Zhejiang University, Hangzhou 310027, China Department of Physics, Zhejiang University, Hangzhou 310027, China In this paper, we present a nonrecursive residual Monte Carlo method for estimating discretization errors associated with the S N transport solution to radiation transport problems. Although this technique is general, we applied it to the mono-energetic 1-D S N equation with linear-discontinuous finite element method spatial discretization as a demonstration of the theory for the purpose of this study. Two angular flux representations: conforming and simplified representations, were considered in this analysis, and the results were compared. It is shown that the simplified representation dramatically reduces the memory footprint and computational complexity of residual source generation and sampling while accurately capturing the error associated with certain types of responses. residual Monte Carlo discretization error angular flux representation Introduction 1 Efficient and accurate solutions to the radiation transport equation are of fundamental importance in the simulation of fission reactor cores, nuclear fusion, and radiation shielding. The Monte Carlo method[ 1 1 ] and deterministic methods (such as S N [ 2 2 , 3 3 ], P N [ 4 4 ], and SP N [ 5 5 - 9 9 5-9 ]) are two different approaches to the solution of the transport equation. The Monte Carlo method is capable of producing fairly accurate results at the cost of an enormous amount of CPU hours, whereas deterministic methods can be performed at a fraction of the cost but often suffer from model errors and discretization errors. Several efforts have been made to combine these two methods to achieve better efficiency and accuracy. The first method falls in the residual Monte Carlo category[ 10 10 , 11 11 ], where the deterministic method is invoked recursively to accelerate the convergence of the Monte Carlo calculation[ 12 12 - 14 14 12-14 ]. Another method can start with a deterministic solution and use Monte Carlo to perform error correction. This new idea is a nonrecursive approach and is expected to exhibit superior performance in problems for which coarse-mesh deterministic calculations are capable of capturing the fundamental structure of the solution. Consider the mono-energetic 1-D neutron transport equation with the Legendre scattering cross-sectional expansion truncated at order L . For simplicity, we assume an isotropic distributed source: <math display="block" id="M1"><mrow><mi>μ</mi><mfrac><mrow><mtext>d</mtext><mi>ψ</mi><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>x</mi></mrow></mfrac><mo>+</mo><msub><mi>Σ</mi><mi>t</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mi>ψ</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi mathvariant="script">l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover></mstyle><mfrac><mrow><mn>2</mn><mi mathvariant="script">l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msubsup><mi>Σ</mi><mi>s</mi><mi mathvariant="script">l</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>ϕ</mi><mi mathvariant="script">l</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>P</mi><mi mathvariant="script">l</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>+</mo><mi>q</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> where μ is the neutron direction cosine, ψ is the neutron angular flux, ∑ t is the macroscopic total cross-section, <math id="M2"><mrow><mstyle displaystyle="true"><mo>∑</mo><mrow><msubsup><mrow/><mi>s</mi><mi mathvariant="script">l</mi></msubsup></mrow></mstyle></mrow></math> is the <math id="M3"><mi mathvariant="script">l</mi></math> -th Legendre moment of the macroscopic scattering cross-section, <math id="M4"><mrow><msub><mi>ϕ</mi><mi mathvariant="script">l</mi></msub></mrow></math> is the angular flux moment of order <math id="M5"><mi mathvariant="script">l</mi></math> , <math id="M6"><mrow><msub><mi>P</mi><mi mathvariant="script">l</mi></msub></mrow></math> is the Legendre polynomial of order <math id="M7"><mi mathvariant="script">l</mi></math> , and q is the distributed source. Eq. (1) can be solved with the 1-D S N method, which is a deterministic method that collocates the transport equation along a set of discrete ordinates: <math display="block" id="M8"><mtable><mtr><mtd><msub><mi>μ</mi><mi>m</mi></msub><mfrac><mrow><mi>d</mi><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msub></mrow><mrow><mi>d</mi><mi>x</mi></mrow></mfrac><mo>+</mo><msub><mi>Σ</mi><mi>t</mi></msub><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msub><mo>=</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi mathvariant="script">l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover></mstyle><mfrac><mrow><mn>2</mn><mi mathvariant="script">l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msubsup><mi>Σ</mi><mi>s</mi><mi mathvariant="script">l</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mi mathvariant="script">l</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msub><mi>P</mi><mi mathvariant="script">l</mi></msub><mo stretchy="false">(</mo><msub><mi>μ</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>q</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mn>,</mn></mtd></mtr><mtr><mtd><mi>m</mi><mo>=</mo><mn>1</mn><mo>⋯</mo><mi>M</mi><mtext> </mtext><mn>,</mn></mtd></mtr></mtable></math> where m is the discrete direction index running from 1 to M and M is the total number of directions. Along each direction, a first-order ordinary differential equation is solved using a spatial discretization technique such as discontinuous finite element method[ 15 15 , 16 16 ]. The equations along different directions are coupled through the scattering term. The ~ in Eq. (2) indicates that the S N solution is an approximation of the true transport solution, in the sense that the S N solution is subject to angular and spatial discretization errors. In the pursuit of a high-fidelity simulation, it is often desirable to estimate the discretization error, given an approximate solution. This process can be computationally costly because it generally entails mesh refinements in a high-dimensional phase space (in this case, space and angle). The required computation can become prohibitively expensive when the dimensionality is high and/or the solution does not exhibit adequate regularity, which can lead to convergence degradation. We propose a novel approach to tackle this challenge using the Monte Carlo method to evaluate the deterministic numerical error. Let us first rewrite the transport equation in operator form as <math display="block" id="M9"><mrow><mstyle mathvariant="bold" mathsize="normal"><mtext>L</mtext></mstyle><mi>ψ</mi><mo>=</mo><mi>q</mi><mtext> </mtext><mn>,</mn></mrow></math> where L denotes the neutron transport operator, ψ denotes the angular flux, and q denotes the distributed source. Given an approximate solution <math id="M10"><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover></math> , the associated residual is given by <math display="block" id="M11"><mrow><mi mathvariant="script">R</mi><mo>=</mo><mi>q</mi><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>L</mi></mstyle><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mtext> </mtext><mn>.</mn></mrow></math> Subsequently, the error <math id="M12"><mrow><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>≡</mo><mi>ψ</mi><mo>−</mo><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover></mrow></math> satisfies the following transport equation: <math display="block" id="M13"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>L</mi></mstyle><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>=</mo><mi mathvariant="script">R</mi><mtext> </mtext><mn>.</mn></mrow></math> We emphasize that L needs to be the true transport operator for Eq. (5) to hold. To numerically solve for <math id="M14"><mrow><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover></mrow></math> , it might be tempting to substitute L with a high-order deterministic transport operator, such as S N , projected onto a highly refined space-angle mesh. However, as stated earlier in this section, this is computationally expensive. The key to this approach is to use the Monte Carlo method to invert L, which is inherently exempt from discretization errors. The residual Monte Carlo method presented here is significantly more efficient than the standard Monte Carlo method, where the full solution is first calculated, and then, the error is computed by subtracting the deterministic solution from the Monte Carlo solution. This is because the residual Monte Carlo method directly computes the error itself and concentrates computational resources towards phase-space regions where the residual is large by sampling source particles from a probability distribution function that is proportional to the absolute value of the residual. We refer to this procedure as a nonrecursive residual Monte Carlo method to differentiate it from the standard version of the residual Monte Carlo method. The standard residual Monte Carlo method involves a recursive set of calculations combined with phase-space adaptive mesh refinement to achieve exponential convergence with the number of particle histories [ 17 17 - 19 19 17-19 ]. Because the residual function can change sign, a source particle sampled from its absolute value ( <math id="M15"><mrow><mo>|</mo><mi mathvariant="script">R</mi><mo>|</mo></mrow></math> ) carries a weight bearing the sign of the original residual function at its sampled phase-space coordinate. Another difficulty caused by changing the sign is associated with the integration of the probability distribution function over the phase space, which is an indispensable step in source sampling. More specifically, nonsmoothness across the sign-changing boundary (see Fig. 1 Fig. 1 ) is problematic for numerical integration. In the case of 1-D transport calculation with linear finite element spatial representation, one can determine the sign boundary analytically and compute the integral using a divide-and-conquer approach. In multi-D, however, computing the boundary is difficult . To this end, we propose a simplified angular flux representation that is isotopic in angle and piecewise constant in space. As will be shown later in this paper, as long as the simplified representation preserves the original S N scalar flux, the resulting error estimation is exact for angle-space integrated responses such as cell-averaged scalar flux and reaction rate. Note that, although the discussion presented in this paper revolves around mono-energetic 1-D transport calculation, the basic principle can be equally applied to multigroup multi-D transport problems. A flow chart for the entire scheme of the residual Monte Carlo method is presented in Fig. 2 Fig. 2 The remainder of this paper is organized as follows. First, we consider the S N equation with standard linear-discontinuous (LD) finite element spatial discretization and generate the associated residual. We then introduce a simplified spatial representation for the angular flux and explain the residual generation using this simplified representation. Next, we discuss some technical details for sampling source particles from the residual function as well as the boundary conditions and current cell interface preservation techniques. Next, we present numerical test results to demonstrate the efficacy of the method . Finally, we summarize our findings and provide recommendations for future work. Standard Linear Discontinuous Finite Element and Residual Generation 1 The prerequisite for the residual Monte Carlo method is a functional representation of the S N angular flux, which is defined at every space-angle point. This construction is not unique. The representation presented below is relatively straightforward and conformal to the discretization scheme. Piecewise Constant Angular Representation 2 We assume that the angular flux is constant within the angle interval associated with each discrete S N direction. That is, <math display="block" id="M16"><mrow><msup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><msup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msup><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mo>−</mo><mn>1</mn><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mi>μ</mi><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mn>1</mn><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M17"><mrow><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><mo>−</mo><mn>1</mn><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>m</mi><mo>−</mo><mn>1</mn></mrow></munderover></mstyle><msub><mi>w</mi><mi>n</mi></msub><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M18"><mrow><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><mo>−</mo><mn>1</mn><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover></mstyle><msub><mi>w</mi><mi>n</mi></msub><mtext> </mtext><mn>,</mn></mrow></math> and w m is the S N quadrature weight associated with discrete direction μ m . Eq. (7) ensure that the scalar flux computed based on the piecewise constant representation is consistent with the S N quadrature rule: <math display="block" id="M19"><mrow><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msup><mi>Δ</mi><msub><mi>μ</mi><mi>m</mi></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>m</mi><mo>=</mo><mn>1</mn></mrow><mi>M</mi></munderover></mstyle><msup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>m</mi></msup><msub><mi>w</mi><mi>m</mi></msub><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M20"><mrow><mi>Δ</mi><msub><mi>μ</mi><mi>m</mi></msub><mo>=</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mtext> </mtext><mn>.</mn></mrow></math> Linear-Discontinuous Finite Element Method Spatial Discretization 2 Linear-discontinuous finite element method (LD-FEM) is a standard spatial discretization scheme for the first-order S N transport equation. In this method, the spatial solution is projected onto a set of basis functions that are linearly continuous within each cell but discontinuous across cell interfaces. In particular, a pair of basis functions ( B L , B R ) for cell i is defined as <math display="block" id="M21"><mrow><msub><mi>B</mi><mrow><mi>i</mi><mn>,</mn><mi>L</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><mi>x</mi></mrow><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M22"><mrow><msub><mi>B</mi><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M23"><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub><mo>=</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mtext> </mtext><mtext> </mtext><mn>.</mn></mrow></math> Fig. 3 Fig. 3 shows the basis functions for cell i . The angular flux solution is then given by a linear combination of these two basis functions: <math display="block" id="M24"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><msub><mi>B</mi><mrow><mi>i</mi><mn>,</mn><mi>L</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>+</mo><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup><msub><mi>B</mi><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mtext> </mtext><mn>,</mn></mrow></math> where i is the cell-center index, <math id="M25"><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>L</mi></mrow><mi>m</mi></msubsup></mrow></math> is the solution degree of freedom (DoF) on the left in cell i , and <math id="M26"><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup></mrow></math> is the solution DoF on the right in cell i . The solution at the cell interface is determined by the standard upwinding. For example, <math display="block" id="M27"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>m</mi></msubsup><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><msub><mi>μ</mi><mi>m</mi></msub><mo>≥</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><msub><mi>μ</mi><mi>m</mi></msub><mo><</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow><mtext> </mtext><mn>,</mn></mrow></math> where <math id="M28"><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></math> is the cell surface index. An illustration of the LD spatial representation for μ > 0 is shown in Fig. 4 Fig. 4 As illustrated above, one of the major advantages of the finite element method is that the solution is defined everywhere in space. Residual Generation with LD-FEM 2 In the case of LD-FEM spatial representation, the residual generation can be conveniently divided into two parts: cell-interior residual and cell-interface residual. To simplify the discussion, let us assume isotropic scattering. By substituting Eq. (14) and Eq. (9) into Eq. (4), we obtain the cell-interior residual. In particular, for cell i and angle bin m , <math display="block" id="M29"><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mn>,</mn><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow><mn>2</mn></mfrac><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mi>i</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>q</mi><mi>i</mi></msub><mo>−</mo><mi>μ</mi><mfrac><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup><mo>−</mo><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>L</mi></mrow><mi>m</mi></msubsup></mrow><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>−</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mn>,</mn><mi>i</mi></mrow></msub><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mn>,</mn><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mtable><mtr><mtd><mrow><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo><</mo><mi>x</mi><mo><</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mn>1</mn><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mi>μ</mi><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mn>1</mn></mrow></mtd></mtr></mtable><mtext> </mtext><mn>.</mn></mrow></math> The cell-interface residual measures discontinuities in the angular flux across the interfaces. The derivative operation on a discontinuity yields a Dirac delta function. A discontinuity at an interface is approached from either side of the interface depending on the upwind direction, which is determined by μ m . In particular, on cell-interface <math id="M30"><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></math> , <math display="block" id="M31"><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mo>−</mo><mi>μ</mi><msup><mi>δ</mi><mo>+</mo></msup><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><mo>−</mo><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup></mrow><mo>)</mo></mrow><mn>,</mn></mrow></mtd><mtd><mrow><mtext>for </mtext><mn>0</mn><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mi>μ</mi><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mo>+</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mi>μ</mi><msup><mi>δ</mi><mo>−</mo></msup><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><mo>−</mo><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup></mrow><mo>)</mo></mrow><mn>,</mn></mrow></mtd><mtd><mrow><mtext>for </mtext><mo>−</mo><mn>1</mn><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>≤</mo><mi>μ</mi><mo>≤</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo><</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow><mn>,</mn></mrow></math> where δ + is the positive one-sided delta that approaches the discontinuity from the right ( <math id="M32"><mrow><mi>x</mi><mo>></mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow></math> ), and δ - is the negative one-sided delta that approaches the discontinuity from the left ( <math id="M33"><mrow><mi>x</mi><mo><</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow></math> ). An illustration of δ + is provided in Fig. 5 Fig. 5 Simplified Angular Flux Representation and Residual Generation 1 Eq. (15) indicates that <math id="M34"><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover></math> can change sign within a phase-space cell. As will be shown in Section 4, the integration of the residual over each phase-space cell needs to be computed as part of the source sampling procedure. As explained in Section 1, the sign change within the phase-space cell poses a challenge to numerical integration schemes. Moreover, the standard conformal angular flux representation requires a solution at each spatial DoF and for each discrete direction. This requires additional memory costs, as typical S N codes only store moment values. Isotropic and Piecewise Constant Angular Flux Representation 2 To simplify the structure of the residual and thereby simplify the source sampling procedure, we propose a simplified angular flux representation that is isotropic in angle and piece(cell)-wise constant in space. Despite of the dramatic reduction in complexity, it can be shown that, as long as the simplified representation preserves the original cell-integrated scalar flux, the error estimation for the cell-integrated responses, such as cell-averaged scalar flux and reaction rate, is rigorous. The simplified angular flux representation is as follows: <math display="block" id="M35"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo><mo>≡</mo><mfrac><mrow><msub><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow><mi>i</mi></msub></mrow><mrow><mn>2</mn><mo>⋅</mo><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M36"><mrow><msub><mrow><mrow><mo>〈</mo><mrow><mtext> </mtext><mo>⋅</mo><mtext> </mtext></mrow><mo>〉</mo></mrow></mrow><mi>i</mi></msub><mo>≡</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></msub></mrow></mstyle><mo stretchy="false">(</mo><mo>⋅</mo><mo stretchy="false">)</mo><mtext> </mtext><mtext> </mtext><mtext>d</mtext><mi>x</mi><mtext> </mtext><mn>.</mn></mrow></math> As a sanity check on the preservation of the cell-integrated scalar flux, <math display="block" id="M37"><mrow><msub><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mo>〉</mo></mrow><mi>i</mi></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></msub></mrow></mstyle><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup></mrow></mstyle><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub><mtext>d</mtext><mi>μ</mi><mtext> </mtext><mtext>d</mtext><mi>x</mi><mo>=</mo><msub><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow><mi>i</mi></msub><mtext> </mtext><mn>.</mn></mrow></math> To show that <math id="M38"><mrow><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>=</mo><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover></mrow></math> , we present the following proof. We begin with the following error equations: <math display="block" id="M39"><mrow><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>+</mo><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>=</mo><mi>ψ</mi><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M40"><mrow><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>+</mo><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>=</mo><mi>ψ</mi><mtext> </mtext><mn>.</mn></mrow></math> By subtracting Eq. (22) from Eq. (21), <math display="block" id="M41"><mrow><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mo>−</mo><mi>δ</mi><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>=</mo><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mo>−</mo><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mtext> </mtext><mn>.</mn></mrow></math> Integrating Eq. (23) over angle and space, we obtain <math display="block" id="M42"><mrow><mi>δ</mi><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo><mo>−</mo><mi>δ</mi><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mo>〉</mo><mo>=</mo><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mo>〉</mo><mo>−</mo><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo><mo>=</mo><mn>0</mn><mtext> </mtext><mn>.</mn></mrow></math> As in the LD case, the angular flux at the cell interface is also determined by upwinding: <math display="block" id="M43"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo>≥</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo><</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow><mtext> </mtext><mn>.</mn></mrow></math> An illustration of the simplified angular flux representation for μ > 0 is provided in Fig. 6 Fig. 6 . Residual Generation with the Simplified Representation 2 The residual generation is significantly simplified owing to the reduced-order angular flux representation. In the interior of the cell, the residual is piecewise constant and isotropic: <math display="block" id="M44"><mrow><msub><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>Σ</mi><mrow><mi>s</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow><mn>2</mn></mfrac><msub><mover accent="true"><mi>ϕ</mi><mo>¯</mo></mover><mi>i</mi></msub><mo>+</mo><msub><mi>q</mi><mi>i</mi></msub><mo>−</mo><msub><mi>Σ</mi><mrow><mi>t</mi><mn>,</mn><mi>i</mi></mrow></msub><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mtable><mtr><mtd><mrow><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo><</mo><mi>x</mi><mo><</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mn>1</mn><mo>≤</mo><mi>μ</mi><mo>≤</mo><mo>+</mo><mn>1</mn></mrow></mtd></mtr></mtable><mtext> </mtext><mn>,</mn></mrow></math> On the cell interfaces, we obtain <math display="block" id="M45"><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mo>+</mo></msubsup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mi>μ</mi><msup><mi>δ</mi><mo>+</mo></msup><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub></mrow><mo>)</mo></mrow><mn>,</mn><mtext> for </mtext><mn>0</mn><mo>≤</mo><mi>μ</mi><mo>≤</mo><mo>+</mo><mn>1</mn></mrow></math> <math display="block" id="M46"><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mo>−</mo></msubsup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mi>μ</mi><msup><mi>δ</mi><mo>−</mo></msup><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mi>i</mi></msub></mrow><mo>)</mo></mrow><mn>,</mn><mtext> for </mtext><mo>−</mo><mn>1</mn><mo>≤</mo><mi>μ</mi><mo><</mo><mn>0</mn></mrow></math> In this case, the cell interior residual does not change sign within a phase-space cell; therefore, the residual integration is straightforward. From the memory saving perspective, in 1-D, the simplified representation requires storing only one DoF per cell for all directions, whereas LD representation stores two DoFs per cell and for each M discrete direction. Hence, 2× M memory savings are achieved. In the case of 2-D and 3-D, the saving factor is 4 × M ( M +2) and 8 × M ( M +2), respectively. Source Sampling for Residual Monte Carlo 1 From the perspective of implementation, the only major difference between the residual Monte Carlo and standard Monte Carlo methods resides in the source sampling, after which the standard Monte Carlo algorithm can be equally applied to the residual Monte Carlo method for tracking the history of the sampled particles. Unlike the standard Monte Carlo method that samples source particles from a given source probability distribution, the residual Monte Carlo method samples particles from a distribution that is proportional to the absolute value of the residual, <math id="M47"><mrow><mo>|</mo><mi mathvariant="script">R</mi><mo>|</mo></mrow></math> . More specifically, the residual source probability distribution is <math id="M48"><mrow><mo>|</mo><mi mathvariant="script">R</mi><mo>|</mo></mrow></math> re-normalized as follows: <math display="block" id="M49"><mrow><mi>P</mi><mi>r</mi><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mo>|</mo><mi mathvariant="script">R</mi><mo>|</mo></mrow><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><mi>x</mi></mrow></msub></mrow></mstyle><mo>|</mo><mi mathvariant="script">R</mi><mo>|</mo><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>μ</mi></mrow></mfrac><mtext> </mtext><mn>.,</mn></mrow></math> In practice, a hierarchical approach is adopted. First, the angle-space cell index ( m, i ) or cell index ( i ) is sampled depending on the angular flux representation. Subsequently, the actual position and direction of the source particle are sampled from within that cell. Discrete Sampling for Angle-Space Cell Index 2 The angle-space cell index is sampled from a probability mass function, which is the discrete version of Eq. (27) . In the case of the standard LD representation, <math display="block" id="M50"><mrow><mi>P</mi><msub><mi>r</mi><mrow><mi>m</mi><mn>,</mn><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>R</mi><mrow><mi>m</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mrow><mi>m</mi><mn>,</mn><mi>i</mi></mrow></msub><mrow><msub><mi>R</mi><mrow><mi>m</mi><mn>,</mn><mi>i</mi></mrow></msub></mrow></mstyle></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M51"><mrow><msub><mi>R</mi><mrow><mi>m</mi><mn>,</mn><mi>i</mi></mrow></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msub><mi>μ</mi><mi>m</mi></msub></mrow></msub></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></msub></mrow></mstyle><mo>|</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo>|</mo><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>μ</mi><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msubsup><mi>μ</mi><mi>m</mi><mo>+</mo></msubsup></mrow></msub></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msubsup><mi>x</mi><mi>i</mi><mo>+</mo></msubsup></mrow></msub></mrow></mstyle><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>μ</mi><mo>+</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msubsup><mi>μ</mi><mi>m</mi><mo>−</mo></msubsup></mrow></msub></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msubsup><mi>x</mi><mi>i</mi><mo>−</mo></msubsup></mrow></msub></mrow></mstyle><mo stretchy="false">(</mo><mo>−</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">)</mo><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>μ</mi><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M52"><mrow><mi>m</mi><mo>=</mo><mn>1,2,</mn><mo>⋯</mo><mn>,</mn><mi>M</mi><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M53"><mrow><mi>i</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mn>,</mn><mtext> </mtext><mn>1,</mn><mtext> </mtext><mfrac><mn>3</mn><mn>2</mn></mfrac><mn>,</mn><mtext> </mtext><mn>2,</mn><mtext> </mtext><mo>⋯</mo><mn>,</mn><mtext> </mtext><mi>I</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mtext> </mtext><mn>,</mn></mrow></math> and <math id="M54"><mrow><mo stretchy="false">(</mo><mi>Δ</mi><msubsup><mi>μ</mi><mi>m</mi><mo>+</mo></msubsup><mn>,</mn><mi>Δ</mi><msubsup><mi>x</mi><mi>i</mi><mo>+</mo></msubsup><mo stretchy="false">)</mo></mrow></math> denotes the portion of the angle-space domain in which <math id="M55"><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover></math> is positive, <math id="M56"><mrow><mo stretchy="false">(</mo><mi>Δ</mi><msubsup><mi>μ</mi><mi>m</mi><mo>−</mo></msubsup><mn>,</mn><mi>Δ</mi><msubsup><mi>x</mi><mi>i</mi><mo>−</mo></msubsup><mo stretchy="false">)</mo></mrow></math> denotes the portion of the angle-space domain in which <math id="M57"><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover></math> is negative, I is the total number of cells, and half-integer indices (e.g., <math id="M58"><mrow><mfrac><mn>3</mn><mn>2</mn></mfrac></mrow></math> ) are the cell interfaces the widths of which are taken to be <math id="M59"><mrow><mi>ϵ</mi><mo>→</mo><mn>0</mn></mrow></math> . In the case of the simplified representation, angle index m is no longer needed, as no angular dependency is present: <math display="block" id="M60"><mrow><mi>P</mi><msub><mi>r</mi><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>R</mi><mi>i</mi></msub></mrow><mrow><mstyle displaystyle="true"><msub><mo>∑</mo><mi>i</mi></msub><mrow><msub><mi>R</mi><mi>i</mi></msub></mrow></mstyle></mrow></mfrac><mtext> </mtext><mn>,</mn></mrow></math> where <math display="block" id="M61"><mrow><msub><mi>R</mi><mi>i</mi></msub><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup></mrow></mstyle><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></msub></mrow></mstyle><mo>|</mo><msub><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mi>i</mi></msub><mo>|</mo><mtext>d</mtext><mi>x</mi><mtext>d</mtext><mi>μ</mi><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M62"><mrow><mi>i</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mn>,</mn><mtext> </mtext><mn>1,</mn><mtext> </mtext><msup><mrow><mfrac><mn>3</mn><mn>2</mn></mfrac></mrow><mo>−</mo></msup><mn>,</mn><mtext> </mtext><msup><mrow><mfrac><mn>3</mn><mn>2</mn></mfrac></mrow><mo>+</mo></msup><mn>,</mn><mtext> </mtext><mn>2,</mn><mtext> </mtext><mo>⋯</mo><mn>,</mn><mtext> </mtext><mi>I</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mtext> </mtext><mn>,</mn></mrow></math> and the +/- sign on the interface indices indicates the right/left side of the interface. Sampling Angle-Space Coordinates from within a Cell 2 Once the angle-space cell index is determined, the actual spatial position and flight direction of the source particle can be sampled. The sampling procedure is divided between the cell interiors and the cell interfaces. Sampling from Cell Interiors 3 In the case of the standard LD representation with isotropic scattering, <math id="M63"><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup></mrow></math> is linear in both μ and x according to Eq. (16) . As shown in Fig. 7 Fig. 7 , the maximum value of <math id="M64"><mrow><mo>|</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo>|</mo></mrow></math> occurs on one of the four phase-space corners. <math display="block" id="M65"><mrow><mi>max</mi><mo stretchy="false">(</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">)</mo><mo>=</mo><mi>max</mi><mrow><mo>(</mo><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">)</mo><mn>,</mn><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">)</mo><mn>,</mn><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">)</mo><mn>,</mn><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow></math> By using <math id="M66"><mrow><mi>max</mi><mo stretchy="false">(</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">)</mo></mrow></math> as the ceiling, we can sample μ and x from <math id="M67"><mrow><mo>|</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo>|</mo></mrow></math> using the rejection sampling technique. Note that higher-order finite elements or higher-order scattering results in a higher-order <math id="M68"><mrow><mo>|</mo><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mo>|</mo></mrow></math> surface, which is significantly more complex to sample from. Fortunately, none of these matters when we use the simplified representation, of which sampling is simply uniform from the angle and space intervals. This is only possible because <math id="M69"><mrow><mo>|</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo>|</mo></mrow></math> is constant in angle and space, regardless of the actual shape of <math id="M70"><mrow><mo>|</mo><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mo>|</mo></mrow></math> within a given cell. Sampling from Cell Interfaces 3 On interfaces, only μ needs to be sampled because x is dictated by the interface index (in 1-D). For the LD representation, μ is sampled from a linear PDF defined within <math id="M71"><mrow><mo stretchy="false">(</mo><msub><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn><msub><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">)</mo></mrow></math> . For the simplified representation, μ is sampled from a similar distribution, either from the right or left side of the interface. Determining Particle Weights 3 Once the coordinates are sampled, the weight of the sampled particle is determined by <math display="block" id="M72"><mrow><mi>w</mi><mi>t</mi><mo>=</mo><mfrac><mrow><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><mo>|</mo><msubsup><mover accent="true"><mi mathvariant="script">R</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo><mo>|</mo></mrow></mfrac><mn>,</mn></mrow></math> for a standard LD representation or <math display="block" id="M73"><mrow><mi>w</mi><mi>t</mi><mo>=</mo><mfrac><mrow><msub><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mi>i</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo></mrow><mrow><mo>|</mo><msub><mover accent="true"><mi mathvariant="script">R</mi><mo>¯</mo></mover><mi>i</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mn>,</mn><mi>x</mi><mo stretchy="false">)</mo><mo>|</mo></mrow></mfrac><mn>,</mn></mrow></math> for the simplified representation. Preserving the Boundary Flux and Interface Partial Current 1 The residual Monte Carlo algorithm should respect the true transport boundary condition to yield accurate error estimates. In our implementation, this is ensured by setting the upstream angular flux as the incoming boundary flux on the left and right ends, denoted by g L ( μ ) and g R ( μ ), respectively. In particular, we specify <math display="block" id="M74"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>m</mi></msubsup><mo>=</mo><msub><mi>g</mi><mtext>L</mtext></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext> </mtext><mtext>for</mtext><mtext> </mtext><msub><mi>μ</mi><mi>m</mi></msub><mo>≥</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M75"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mrow><mi>I</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mi>m</mi></msubsup><mo>=</mo><msub><mi>g</mi><mtext>R</mtext></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext> </mtext><mtext>for</mtext><mtext> </mtext><msub><mi>μ</mi><mi>m</mi></msub><mo><</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mrow></math> or <math display="block" id="M76"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msub><mi>g</mi><mtext>L</mtext></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo>≥</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M77"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>I</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><msub><mi>g</mi><mtext>R</mtext></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo><</mo><mn>0</mn><mtext> </mtext><mn>,</mn></mrow></math> depending on the used angular flux representation. The partial current across cell interfaces is another aspect that may concern some nuclear engineering applications. Similarly to the aforementioned discussion on cell-averaged scalar flux, in order for the simplified representation to yield an accurate error estimate for the S N partial current, it must preserve the original S N partial current across cell interfaces. This can be achieved by introducing a double discontinuity at the interface. Specifically, instead of upwinding, one restricts the domain of the cell interior angular flux to be strictly within the cell and defines the single-point half-range isotropic angular flux at the interface by another DoF that preserves the S N partial current. That is, we enforce <math display="block" id="M78"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mn>1</mn></msubsup></mrow></mstyle><mi>μ</mi><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mtext>d</mtext><mi>μ</mi><mo>=</mo><mfrac><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mn>2</mn></mfrac><mo>≡</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mn>1</mn></msubsup></mrow></mstyle><mi>μ</mi><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>μ</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>m</mi><mo>|</mo><msub><mi>μ</mi><mi>m</mi></msub><mo>≥</mo><mn>0</mn></mrow></munder></mstyle><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup><mrow><mo>(</mo><mrow><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup><mo>−</mo><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mtext> </mtext><mn>,</mn></mrow></math> <math display="block" id="M79"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>0</mn></msubsup></mrow></mstyle><mi>μ</mi><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mtext>d</mtext><mi>μ</mi><mo>=</mo><mfrac><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mn>2</mn></mfrac><mo>≡</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>0</mn></msubsup></mrow></mstyle><mi>μ</mi><msub><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>μ</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>m</mi><mo>|</mo><msub><mi>μ</mi><mi>m</mi></msub><mo><</mo><mn>0</mn></mrow></munder></mstyle><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><mrow><mo>(</mo><mrow><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup><mo>−</mo><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mtext> </mtext><mn>.</mn></mrow></math> Therefore: <math display="block" id="M80"><mrow><msub><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>m</mi><mo>|</mo><msub><mi>μ</mi><mi>m</mi></msub><mo>≥</mo><mn>0</mn></mrow></munder></mstyle><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mn>,</mn><mi>R</mi></mrow><mi>m</mi></msubsup><mrow><mo>(</mo><mrow><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup><mo>−</mo><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo>≥</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mstyle displaystyle="true"><munder><mo>∑</mo><mrow><mi>m</mi><mo>|</mo><msub><mi>μ</mi><mi>m</mi></msub><mo><</mo><mn>0</mn></mrow></munder></mstyle><msubsup><mover accent="true"><mi>Ψ</mi><mo>˜</mo></mover><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>L</mi></mrow><mi>m</mi></msubsup><mrow><mo>(</mo><mrow><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup><mo>−</mo><msubsup><mi>μ</mi><mrow><mi>m</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mn>,</mn><mtext> </mtext></mrow></mtd><mtd><mrow><mtext>for</mtext><mtext> </mtext><mi>μ</mi><mo><</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow><mtext> </mtext><mn>.</mn></mrow></math> A double discontinuity for the case of μ ≥ 0 is illustrated in Fig. 8 Fig. 8 The double-discontinuity scheme is only presented here for completeness and was not tested in our numerical experiments because we were only concerned with the error in the cell-averaged scalar flux. Numerical Experiment 1 Manufactured Solution Test 2 Problem Setting 3 We verified the efficacy of the proposed residual Monte Carlo method in a test problem defined over a domain of <math id="M81"><mrow><mi>x</mi><mo>∈</mo></mrow></math> [0, 10]cm. A uniform mesh size was set to 0.5cm. The material properties were chosen as ∑ t = 1.0cm -1 and ∑ s = 0.5cm -1 with isotropic scattering. To construct a distributed source together with accompanying boundary conditions from a target isotropic angular flux solution, the manufactured solution technique was used. In the used target solution, a constant value of 3.0n/cm 2 s in the interval [0, 5]cm and 1.0n/cm 2 s in the interval [6,10]cm, with a linear transition region in-betwe en, were taken. We then assumed an S 4 solution that is linear-discontinuous and contains discontinuities at x =0.5, 2.0, and 2.5cm. The angular distribution was assumed to be linearly anisotropic: <math display="block" id="M82"><mrow><msubsup><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover><mi>i</mi><mi>m</mi></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mi>i</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>x</mi><mo>−</mo><msub><mover accent="true"><mi>x</mi><mo>¯</mo></mover><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mi>Δ</mi><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext>for</mtext><mtext> </mtext><mtext> </mtext><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo><</mo><mi>x</mi><mo><</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mn>,</mn></mrow></math> where <math display="block" id="M83"><mrow><msub><mover accent="true"><mi>x</mi><mo>¯</mo></mover><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>x</mi><mrow><mi>i</mi><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub><mo>+</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>+</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msub></mrow><mn>2</mn></mfrac><mn>.</mn></mrow></math> Eq. (45) shows that anisotropy is modulated by a spatial dependency such that the angular flux gets more inward-pointed as it approaches the boundaries. In particular, leakage through the boundaries is zero. Numerical Results and Remarks 3 The objective was to compute the error in the cell-integrated scalar flux <math id="M84"><mrow><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow></mrow></math> . We first computed the error by using the full LD S N angular flux. The results are shown in Fig. 9 Fig. 9 . This is a dual-y-axis plot. On the left axis, ϕ MMS is the manufactured solution, <math id="M85"><mrow><msub><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><msub><mi>S</mi><mi>N</mi></msub></mrow></msub></mrow></math> is the LD S N solution, <math id="M86"><mrow><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow></mrow></math> is the cell-integrated S N scalar flux, <math id="M87"><mrow><mi>δ</mi><msub><mrow><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow></mrow><mrow><mtext>MC</mtext></mrow></msub></mrow></math> is the error computed by the residual Monte Carlo method, and <math id="M88"><mrow><mi>σ</mi><mo>:</mo><mi>δ</mi><msub><mrow><mrow><mo>〈</mo><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mo>〉</mo></mrow></mrow><mrow><mtext>MC</mtext></mrow></msub></mrow></math> is the standard deviation associated with the error estimation. On the right axis, <math id="M89"><mi mathvariant="script">R</mi></math> cell is the cell-interior residual, <math id="M90"><mi mathvariant="script">R</mi></math> interface is the cell-interface residual, and Q is the forward source. For this particular test, 10 6 particles were run in ten batches, with 10 5 particles per batch. It can be seen that the residual Monte Carlo method is capable of quantifying the cell-integrated scalar flux error fairly accurately. With 10 6 particles, the standard deviation is so small that it is difficult to distinguish between the upper and lower bounds of the error bar. The same error was then computed using the isotropic cell-wise constant representation for the angular flux <math id="M91"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover></math> . The results are shown in Fig. 10 Fig. 10 , where the quantities plotted are analogous to those shown in Fig. 9 Fig. 9 . The errors shown in Fig. 10 Fig. 10 are statistically identical to those found in Fig. 9 Fig. 9 in the sense that the difference is within 1- σ . This implies that the error in the cell-integrated scalar flux can be computed precisely when the simplified angular flux representation is substituted for the full S N angular flux. In both Fig. 9 Fig. 9 and Fig. 10 Fig. 10 , the cell-interior residual ( <math id="M92"><mi mathvariant="script">R</mi></math> cell ) is large when the error is large, and the cell-interface residual ( <math id="M93"><mi mathvariant="script">R</mi></math> interface ) is nonzero only if a discontinuity is present in the S N solution. It should be pointed out that for simplicity of illustration, we only show the interface residuals on the x + side, whereas in reality, interface residuals exist on both x + and x - sides. The forward source (Q) may look unphysical between 4cm and 8cm, as it takes on negative values, which is not unexpected considering that it is a hypothetical source generated from the manufactured solution. To illustrate the efficiency of the residual Monte Carlo method compared to the standard Monte Carlo method, the standard deviations are analyzed. Fig. 11 Fig. 11 plots the standard deviation associated with the residual Monte Carlo calculation against those associated with the standard Monte Carlo calculations. From Fig. 11 Fig. 11 , we can see that when running with 1×10 6 particles, σ MC is generally significantly higher than σ RMC . If we focus on a particular cell, for example, Cell 5 ([2, 2.5]cm), the difference is about a factor of 5. To match the performance of the residual Monte Carlo calculation in Cell 5, the standard Monte Carlo calculation should be performed with 25 times the number of particles, which translates to 25 times the CPU time. Boundary Source Test 2 For a more realistic test, we examined the incident flux problem. The problem domain was still [0, 10]cm, whereas the background material was set up to mimic graphite, with neutron cross-sections ∑ t = 4.9435cm -1 and ∑ s = 4.9400cm -1 . The isotropic incident neutron angular flux entered through the left boundary, ψ ( μ < 0, x =0)=20n/cm 2 s. Vacuum boundary conditions were assumed for the right boundary. There was no distributed source throughout the entire problem domain. The results are presented in Fig. 12 Fig. 12 . A zoomed-in plot on [0,2]cm is shown in Fig. 13 Fig. 13 . In these plots, ϕ S N is the piecewise linear scalar flux solution given by an S 2 calculation, ⟨ ϕ ⟩ S N is the cell-integrated scalar flux, ⟨ ϕ ⟩ MC is the same scalar flux computed by the standard Monte Carlo method, and δ ⟨ ϕ ⟩ RMC is the error computed by the residual Monte Carlo method. It can be seen in these figures, especially in the zoomed-in figure, that the residual Monte Carlo method is capable of estimating the additive error associated with the S 2 solution. Conclusion and Future Work 1 In this paper, we present a nonrecursive residual Monte Carlo method for estimating discretization error associated with the S N transport calculation. For cases in which the quantity of interest is an angle-space integrated response, we propose a simplified angular flux representation that greatly simplifies the sampling process while still yielding the exact S N error. We demonstrated the efficacy of our method with a 1-D S 4 test problem. The results show that our methods are effective at estimating the S N errors, and for the case of the angle-space integrated response, the simplified method yields exactly the same results as the full S N method. Moreover, we propose a modified simplified angular flux representation with double discontinuity, which can be employed to preserve the S N error in partial currents across cell interfaces. It should be noted that, although the discussion presented herein is based on LD-FEM spatial discretization, the general methodology can be equally applied to other spatial discretization schemes. In the case of continuous finite element method, which is generally used for the solution of second-order transport equations[ 20 20 ], discontinuities of the angular flux at the cell interfaces disappear, which is the same for the interface residuals. However, the simplified representation is always recommended because of its simplicity when handling residual integral and source sampling. The angular flux discontinuity is inherent to the simplified representation, and we believe that the hassle is worth the benefit. In future work, we plan to extend the 1-D demonstration to multi-D. The advantage of the residual Monte Carlo method can become more profound as the dimensionality increases. This is partly because the fine-mesh SN calculation in multi-D is extremely expensive. Full convergence to an SN solution is unrealistic. Therefore, a viable way of estimating the discretization or iteration error at a fraction of the SN computational cost can be valuable. In addition, we would like to integrate the residual Monte Carlo method into a parallel S N solver framework to determine whether effective parallelization within the existing infrastructure can be achieved. References J. Spanier , E.M. Gelbard , Monte Carlo Principles and Neutron Transport Problems , ( Dover Books , Mineola, New York , 2008 ) Monte Carlo Principles and Neutron Transport Problems E.E. Lewis , W.F. Miller , Computational Methods of Neutron Transport , ( American Nuclear Society , La Grange Park, Illinois , 1993 ) Computational Methods of Neutron Transport J.E. Morel , B.T. Adams , T. Noh , et al. , Spatial discretizations for self-adjoint forms of the radiative transfer equations . J. Comput. Phys. 214 , 12 - 40 ( 2006 ). doi: 10.1016/j.jcp.2005.09.017 http://doi.org/10.1016/j.jcp.2005.09.017 Spatial discretizations for self-adjoint forms of the radiative transfer equations E.M. Gelbard , Applications of spherical harmonics method to reactor problems. Tech. Rep. WAPD-BT-20 , Bettis Atomic Power Laboratory ( 1960 ) Applications of spherical harmonics method to reactor problems. Tech. Rep. WAPD-BT-20 E.M. Gelbard , Simplified spherical harmonics equations and their use in shielding problems. Tech. Rep. WAPD-T-1182 , Bettis Atomic Power Laboratory ( 1961 ) Simplified spherical harmonics equations and their use in shielding problems. Tech. Rep. WAPD-T-1182 E.M. Gelbard , Applications of simplified spherical harmonics equations in spherical geometry. Tech. Rep. WAPD-TM-294 , Bettis Atomic Power Laboratory ( 1962 ) Applications of simplified spherical harmonics equations in spherical geometry. Tech. Rep. WAPD-TM-294 E.W. Larsen , J.E. Morel , J.M. McGhee , Asymptotic derivation of the multigroup p1 and simplified pn equations with anisotropic scattering . Nucl. Sci. Eng. 123 , 328 - 342 ( 1996 ). Asymptotic derivation of the multigroup p1 and simplified pn equations with anisotropic scattering G.C. Pomraning , Asymptotic and variational derivations of the simplified pn equations . Ann. Nucl. Energy 20 , 623 - 637 ( 1993 ). Asymptotic and variational derivations of the simplified pn equations Y. Zhang , J.C. Ragusa , J.E. Morel , Iterative performance of various formulations of the spn equations . J. Comput. Phys. 252 , 558 - 572 ( 2013 ). doi: 10.1016/j.jcp.2013.06.009 http://doi.org/10.1016/j.jcp.2013.06.009 Iterative performance of various formulations of the spn equations T. Evans , T. Urbatsch , H. Lichtenstein , et al. , A residual monte carlo method for discrete thermal radiative diffusion . J. Comput. Phys. 189 , 539 - 556 ( 2003 ). doi: 10.1016/S0021-9991(03)00233-X http://doi.org/10.1016/S0021-9991(03)00233-X A residual monte carlo method for discrete thermal radiative diffusion J. Willert , H. Park , Residual monte carlo high-order solver for moment-based accelerated thermal radiative transfer equations . J. Comput. Phys. 276 , 405 - 421 ( 2014 ). doi: 10.1016/j.jcp.2014.07.039 http://doi.org/10.1016/j.jcp.2014.07.039 Residual monte carlo high-order solver for moment-based accelerated thermal radiative transfer equations J.H. Halton , Sequential monte carlo . Math. Proc. Cambridge. 58 , 57 - 78 ( 1962 ). doi: 10.1017/S0305004100036227 http://doi.org/10.1017/S0305004100036227 Sequential monte carlo J. Halton , Sequential monte carlo techniques for the solution of linear systems . J. Sci. Comput. 9 ,. Sequential monte carlo techniques for the solution of linear systems L. Chacón , G. Chen , D. Knoll , et al. , Multiscale high-order/low-order (holo) algorithms and applications . J. Comput. Phys. 330 , 21 - 45 ( 2017 ). doi: 10.1016/j.jcp.2016.10.069 http://doi.org/10.1016/j.jcp.2016.10.069 Multiscale high-order/low-order (holo) algorithms and applications W.H. Reed , T.R. Hill , F.W. Brinkley , et al. , Triplet: A two-dimensional, multigroup, triangular mesh, planar geometry, explicit transport code. Tech. Rep. LA-5428-MS , Los Alamos Scientific Laboratory ( 1973 ) Triplet: A two-dimensional, multigroup, triangular mesh, planar geometry, explicit transport code. Tech. Rep. LA-5428-MS J. Morel , T. Wareing , K. Smith , A linear-discontinuous spatial difference scheme for sn radiative transfer calculations . J. Sci. Comput. 128 , 445 - 462 ( 1996 ). A linear-discontinuous spatial difference scheme for sn radiative transfer calculations J.A. Favorite , H. Lichtenstein , Exponential monte carlo convergence of a three-dimensional discrete ordinates solution. Tech. Rep. LA-UR-99-3198 , Los Alamos Scientific Laboratory ( 1999 ) Exponential monte carlo convergence of a three-dimensional discrete ordinates solution. Tech. Rep. LA-UR-99-3198 J.R. Peterson , J.E. Morel , J.C. Ragusa , Residual monte carlo for the one-dimensional neutron transport equation . SIAM J. Sci. Comput. 38 , B941 - B96 ( 2016 ). doi: 10.1137/16M1057619 http://doi.org/10.1137/16M1057619 Residual monte carlo for the one-dimensional neutron transport equation S. Bolding , M. Cleveland , J. Morel , A high-order low-order algorithm with exponentially-convergent monte carlo for thermal radiative transfer . Nuclear Science & Engineering: M&C 2015 Special Issue . A high-order low-order algorithm with exponentially-convergent monte carlo for thermal radiative transfer Y. Zhang , J.E. Morel , J.C. Ragusa , Convergence behavior of second-order transport equations in near-void problems . J. Quant. Spectrosc. Ra. 244 , 106843 ( 2020 ). doi: 10.1016/j.jqsrt.2020.106843 http://doi.org/10.1016/j.jqsrt.2020.106843 Convergence behavior of second-order transport equations in near-void problems 10.1007/s41365-022-01042-w 2021-12-19 2022-03-17 2022-04-05 2022-06-16 2022-07-26 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01042-w http://dx.doi.org/10.1007/s41365-022-01042-w . 2022 List of symbols ψ Angular flux ϕ Scalar flux μ Directional cosine x x -coordinate ∑ t Total cross-section ∑ s Scattering cross-section <math id="M94"><mi mathvariant="script">l</mi></math> Legendre expansion/moment order <math id="M95"><mrow><msub><mi>P</mi><mi mathvariant="script">l</mi></msub></mrow></math> l -th order Legendre moment m Discrete direction index M Total number of directions q Distributed source <math id="M96"><mover accent="true"><mi>ψ</mi><mo>˜</mo></mover></math> S N angular flux solution L Transport operator <math id="M97"><mi mathvariant="script">R</mi></math> Residual w Angular quadrature weight i Cell index B Finite element basis function Ψ Finite element coefficient for angular flux δ +/- ( x ) Positive/negative one-sided Dirac delta function <math id="M98"><mover accent="true"><mi>ψ</mi><mo>¯</mo></mover></math> Simplified representation of angular flux ⟨⋅⟩ Cell-integrated quantity δψ Angular flux error δ ⟨ϕ⟩ Cell-integrated scalar flux error Pr Residual source probability distribution R Cell-angle integrated absolute residual I Total number of cells wt Weight of sampled particle g Angular flux boundary condition σ Standard deviation Np Number of particles LD-FEM Linear-discontinuous finite element C-FEM Continuous finite element Nuclear Science and Techniques 05 33