Introduction

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] and deterministic methods (such as SN[2, 3], PN[4], and SPN[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, 11], where the deterministic method is invoked recursively to accelerate the convergence of the Monte Carlo calculation[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: $\mu \frac{\text{d}\psi (x\mathrm{,}\mu )}{\text{d}x}+{\Sigma }_t(x)\psi ={\displaystyle \sum _{\mathcal{l}=0}^L}\frac{2\mathcal{l}+1}{2}{\Sigma }_s^{\mathcal{l}}(x){\varphi }_{\mathcal{l}}(x)P_{\mathcal{l}}(\mu )+q(x)\mathrm{,}$ (1) where μ is the neutron direction cosine, ψ is the neutron angular flux, ∑t is the macroscopic total cross-section, ${\displaystyle \sum {}_s^{\mathcal{l}}}$ is the $\mathcal{l}$-th Legendre moment of the macroscopic scattering cross-section, ${\varphi }_{\mathcal{l}}$ is the angular flux moment of order $\mathcal{l}$, $P_{\mathcal{l}}$ is the Legendre polynomial of order $\mathcal{l}$, and q is the distributed source. Eq. (1) can be solved with the 1-D SN method, which is a deterministic method that collocates the transport equation along a set of discrete ordinates: $\begin{array}{c}{\mu }_m\frac{d{\tilde{\psi }}_m}{dx}+{\Sigma }_t{\tilde{\psi }}_m=={\displaystyle \sum _{\mathcal{l}=0}^L}\frac{2\mathcal{l}+1}{2}{\Sigma }_s^{\mathcal{l}}(x){\tilde{\varphi }}_{\mathcal{l}}(x)P_{\mathcal{l}}({\mu }_m)+q(x)\mathrm{,}\\ m=1\cdots M\mathrm{,}\end{array}$ (2) 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, 16]. The equations along different directions are coupled through the scattering term. The ~ in Eq. (2) indicates that the SN solution is an approximation of the true transport solution, in the sense that the SN 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 $\text{L}\psi =q\mathrm{,}$ (3) where L denotes the neutron transport operator, ψ denotes the angular flux, and q denotes the distributed source. Given an approximate solution $\tilde{\psi }$, the associated residual is given by $\mathcal{R}=q-L\tilde{\psi }\mathrm{.}$ (4) Subsequently, the error $\delta \tilde{\psi }\equiv \psi -\tilde{\psi }$ satisfies the following transport equation: $L\delta \tilde{\psi }=\mathcal{R}\mathrm{.}$ (5) We emphasize that L needs to be the true transport operator for Eq. (5) to hold. To numerically solve for $\delta \tilde{\psi }$, it might be tempting to substitute L with a high-order deterministic transport operator, such as SN, 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-19].

Because the residual function can change sign, a source particle sampled from its absolute value ($\left|\mathcal{R}\right|$) 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) 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 SN scalar flux, the resulting error estimation is exact for angle-space integrated responses such as cell-averaged scalar flux and reaction rate.

Fig. 1Full-screen View

Nonsmoothness in the absolute value of 𝒭

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. 2Full-screen View

Residual Monte Carlo flow chart

The remainder of this paper is organized as follows. First, we consider the SN 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

The prerequisite for the residual Monte Carlo method is a functional representation of the SN 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.

2.2

Linear-Discontinuous Finite Element Method Spatial Discretization

Linear-discontinuous finite element method (LD-FEM) is a standard spatial discretization scheme for the first-order SN 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 (BL, BR) for cell i is defined as $B_{i\mathrm{,}L}(x)=\frac{x_{i+\frac{1}{2}}-x}{\Delta x_i}\mathrm{,}$ (10) $B_{i\mathrm{,}R}(x)=\frac{x-x_{i-\frac{1}{2}}}{\Delta x_i}\mathrm{,}$ (11) $\Delta x_i=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}\mathrm{.}$ (12)

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: ${\tilde{\psi }}_i^m(x)={\tilde{\Psi }}_{i\mathrm{,}L}^m{B}_{i\mathrm{,}L}(x)+{\tilde{\Psi }}_{i\mathrm{,}R}^m{B}_{i\mathrm{,}R}(x)\mathrm{,}$ (13) where i is the cell-center index, ${\tilde{\Psi }}_{i\mathrm{,}L}^m$ is the solution degree of freedom (DoF) on the left in cell i, and ${\tilde{\Psi }}_{i\mathrm{,}R}^m$ is the solution DoF on the right in cell i.

Fig. 3Full-screen View

Linear-discontinuous basis functions for cell i

The solution at the cell interface is determined by the standard upwinding. For example, ${\tilde{\psi }}_{i+\frac{1}{2}}^m=\{\begin{array}{cc}{\tilde{\Psi }}_{i\mathrm{,}R}^m\mathrm{,}& \text{for}{\mu }_m\ge 0\\ {\tilde{\Psi }}_{i+\mathrm{1,}L}^m\mathrm{,}& \text{for}{\mu }_m<0\end{array}\mathrm{,}$ (14) where $i+\frac{1}{2}$ is the cell surface index. An illustration of the LD spatial representation for μ>0 is shown in Fig. 4

Fig. 4Full-screen View

Linear-discontinuous spatial representations with upwinding

As illustrated above, one of the major advantages of the finite element method is that the solution is defined everywhere in space.

2.3

Residual Generation with LD-FEM

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, ${\tilde{\mathcal{R}}}_i^m(x\mathrm{,}\mu )=\frac{{\Sigma }_{s\mathrm{,}i}}{2}{\tilde{\varphi }}_i(x)+q_i-\mu \frac{{\tilde{\Psi }}_{i\mathrm{,}R}^m-{\tilde{\Psi }}_{i\mathrm{,}L}^m}{\Delta x_i}-{\Sigma }_{t\mathrm{,}i}{\tilde{\Psi }}_i^m(x)\mathrm{,}\text{for}\begin{array}{c}{x}_{i-\frac{1}{2}}<x<x_{i+\frac{1}{2}}\\ -1\le {\mu }_{m-\frac{1}{2}}\le \mu \le {\mu }_{m+\frac{1}{2}}\le 1\end{array}\mathrm{.}$ (15) 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 $i+\frac{1}{2}$, ${\tilde{\mathcal{R}}}_{i+\frac{1}{2}}^m(\mu )=\{\begin{array}{cc}-\mu {\delta }^{+}\left(x-x_{i+\frac{1}{2}}\right)\left({\tilde{\Psi }}_{i+\mathrm{1,}L}^m-{\tilde{\Psi }}_{i\mathrm{,}R}^m\right)\mathrm{,}& \text{for }0\le {\mu }_{m-\frac{1}{2}}\le \mu \le {\mu }_{m+\frac{1}{2}}\le +1\\ -\mu {\delta }^{-}\left(x-x_{i+\frac{1}{2}}\right)\left({\tilde{\Psi }}_{i+\mathrm{1,}L}^m-{\tilde{\Psi }}_{i\mathrm{,}R}^m\right)\mathrm{,}& \text{for }-1\le {\mu }_{m-\frac{1}{2}}\le \mu \le {\mu }_{m+\frac{1}{2}}<0\end{array}\mathrm{,}$ (16) where δ⁺ is the positive one-sided delta that approaches the discontinuity from the right ($x>x_{i+\frac{1}{2}}$), and δ^- is the negative one-sided delta that approaches the discontinuity from the left ($x<x_{i+\frac{1}{2}}$). An illustration of δ⁺ is provided in Fig. 5

Fig. 5Full-screen View

Positive one-sided delta function

Simplified Angular Flux Representation and Residual Generation

Eq. (15) indicates that $\tilde{\mathcal{R}}$ 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 SN codes only store moment values.

3.1

Isotropic and Piecewise Constant Angular Flux Representation

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: ${\overline{\psi }}_i(\mu \mathrm{,}x)\equiv \frac{{\langle \tilde{\varphi }\rangle }_i}{2\cdot \Delta x_i}\mathrm{,}$ (17) where ${\langle \cdot \rangle }_i\equiv {\displaystyle {\int }_{\Delta x_i}}(\cdot )\text{d}x\mathrm{.}$ (18) As a sanity check on the preservation of the cell-integrated scalar flux, ${\langle \overline{\varphi }\rangle }_i={\displaystyle {\int }_{\Delta x_i}}{\displaystyle {\int }_{-1}^1}{\overline{\psi }}_i\text{d}\mu \text{d}x={\langle \tilde{\varphi }\rangle }_i\mathrm{.}$ (19) To show that $\delta \tilde{\psi }=\delta \overline{\psi }$, we present the following proof.

We begin with the following error equations: $\tilde{\psi }+\delta \tilde{\psi }=\psi \mathrm{,}$ (20) $\overline{\psi }+\delta \overline{\psi }=\psi \mathrm{.}$ (21) By subtracting Eq. (22) from Eq. (21), $\delta \tilde{\psi }-\delta \overline{\psi }=\overline{\psi }-\tilde{\psi }\mathrm{.}$ (22) Integrating Eq. (23) over angle and space, we obtain $\delta \langle \tilde{\varphi }\rangle -\delta \langle \overline{\varphi }\rangle =\langle \overline{\varphi }\rangle -\langle \tilde{\varphi }\rangle =0\mathrm{.}$

As in the LD case, the angular flux at the cell interface is also determined by upwinding: ${\overline{\psi }}_{i+\frac{1}{2}}(\mu )=\{\begin{array}{cc}{\overline{\psi }}_i\mathrm{,}& \text{for}\mu \ge 0\\ {\overline{\psi }}_{i+1}\mathrm{,}& \text{for}\mu <0\end{array}\mathrm{.}$ (23) An illustration of the simplified angular flux representation for μ>0 is provided in Fig. 6.

Fig. 6Full-screen View

Isotropic and piecewise constant angular flux representations with upwinding

Source Sampling for Residual Monte Carlo

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, $\left|\mathcal{R}\right|$. More specifically, the residual source probability distribution is $\left|\mathcal{R}\right|$ re-normalized as follows: $Pr(\mu \mathrm{,}x)=\frac{\left|\mathcal{R}\right|}{{\displaystyle {\int }_{-1}^1}{\displaystyle {\int }_{\Delta x}}\left|\mathcal{R}\right|\text{d}x\text{d}\mu }\mathrm{.,}$ (27) 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.

4.2

Sampling Angle-Space Coordinates from within a Cell

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.

4.2.1

Sampling from Cell Interiors

In the case of the standard LD representation with isotropic scattering, ${\tilde{\mathcal{R}}}_i^m$ is linear in both μ and x according to Eq. (16) . As shown in Fig. 7, the maximum value of $\left|{\tilde{\mathcal{R}}}_i^m\right|$ occurs on one of the four phase-space corners. $\max ({\tilde{\mathcal{R}}}_i^m)=\max \left({\tilde{\mathcal{R}}}_i^m({\mu }_{m-\frac{1}{2}}\mathrm{,}{x}_{i-\frac{1}{2}})\mathrm{,}{\tilde{\mathcal{R}}}_i^m({\mu }_{m-\frac{1}{2}}\mathrm{,}{x}_{i+\frac{1}{2}})\mathrm{,}{\tilde{\mathcal{R}}}_i^m({\mu }_{m+\frac{1}{2}}\mathrm{,}{x}_{i-\frac{1}{2}})\mathrm{,}{\tilde{\mathcal{R}}}_i^m({\mu }_{m+\frac{1}{2}}\mathrm{,}{x}_{i+\frac{1}{2}})\right)$ (35)

Fig. 7Full-screen View

LD residual and its absolute value in an angle-space cell

By using $\max ({\tilde{\mathcal{R}}}_i^m)$ as the ceiling, we can sample μ and x from $\left|{\tilde{\mathcal{R}}}_i^m\right|$ using the rejection sampling technique. Note that higher-order finite elements or higher-order scattering results in a higher-order $\left|\tilde{\mathcal{R}}\right|$ 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 $\left|{\overline{\mathcal{R}}}_i^m\right|$ is constant in angle and space, regardless of the actual shape of $\left|\tilde{\mathcal{R}}\right|$ within a given cell.

4.2.2

Sampling from Cell Interfaces

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 $({\mu }_{m-\frac{1}{2}}\mathrm{,}{\mu }_{m+\frac{1}{2}})$. For the simplified representation, μ is sampled from a similar distribution, either from the right or left side of the interface.

4.2.3

Determining Particle Weights

Once the coordinates are sampled, the weight of the sampled particle is determined by $wt=\frac{{\tilde{\mathcal{R}}}_i^m(\mu \mathrm{,}x)}{|{\tilde{\mathcal{R}}}_i^m(\mu \mathrm{,}x)|}\mathrm{,}$ (36) for a standard LD representation or $wt=\frac{{\overline{\mathcal{R}}}_i(\mu \mathrm{,}x)}{|{\overline{\mathcal{R}}}_i(\mu \mathrm{,}x)|}\mathrm{,}$ (37) for the simplified representation.

Preserving the Boundary Flux and Interface Partial Current

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 ${\tilde{\psi }}_{\frac{1}{2}}^m=g_{\text{L}}(\mu )\text{for}{\mu }_m\ge 0\mathrm{,}$ (38) ${\tilde{\psi }}_{I+\frac{1}{2}}^m=g_{\text{R}}(\mu )\text{for}{\mu }_m<0\mathrm{,}$ (39) or ${\overline{\psi }}_{\frac{1}{2}}=g_{\text{L}}(\mu )\text{for}\mu \ge 0\mathrm{,}$ (40) ${\overline{\psi }}_{I+\frac{1}{2}}=g_{\text{R}}(\mu )\text{for}\mu <0\mathrm{,}$ (41) 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 SN partial current, it must preserve the original SN 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 SN partial current. That is, we enforce ${\displaystyle {\int }_0^1}\mu {\overline{\psi }}_{i+\frac{1}{2}}\text{d}\mu =\frac{{\overline{\psi }}_{i+\frac{1}{2}}}{2}\equiv {\displaystyle {\int }_0^1}\mu {\tilde{\psi }}_{i+\frac{1}{2}}(\mu )\text{d}\mu =\frac{1}{2}{\displaystyle \sum _{m|{\mu }_m\ge 0}}{\tilde{\Psi }}_{i\mathrm{,}R}^m\left({\mu }_{m+\frac{1}{2}}^2-{\mu }_{m-\frac{1}{2}}^2\right)\mathrm{,}$ (42) ${\displaystyle {\int }_{-1}^0}\mu {\overline{\psi }}_{i+\frac{1}{2}}\text{d}\mu =\frac{{\overline{\psi }}_{i+\frac{1}{2}}}{2}\equiv {\displaystyle {\int }_{-1}^0}\mu {\tilde{\psi }}_{i+\frac{1}{2}}(\mu )\text{d}\mu =\frac{1}{2}{\displaystyle \sum _{m|{\mu }_m<0}}{\tilde{\Psi }}_{i+\mathrm{1,}L}^m\left({\mu }_{m+\frac{1}{2}}^2-{\mu }_{m-\frac{1}{2}}^2\right)\mathrm{.}$ (43) Therefore: ${\overline{\psi }}_{i+\frac{1}{2}}=\{\begin{array}{cc}{\displaystyle \sum _{m|{\mu }_m\ge 0}}{\tilde{\Psi }}_{i\mathrm{,}R}^m\left({\mu }_{m+\frac{1}{2}}^2-{\mu }_{m-\frac{1}{2}}^2\right)\mathrm{,}& \text{for}\mu \ge 0\\ {\displaystyle \sum _{m|{\mu }_m<0}}{\tilde{\Psi }}_{i+\mathrm{1,}L}^m\left({\mu }_{m+\frac{1}{2}}^2-{\mu }_{m-\frac{1}{2}}^2\right)\mathrm{,}& \text{for}\mu <0\end{array}\mathrm{.}$ (44) A double discontinuity for the case of μ ≥ 0 is illustrated in Fig. 8

Fig. 8Full-screen View

Simplified angular flux representation with double discontinuity for preserving partial current

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

6.1

Manufactured Solution Test

6.1.1

Problem Setting

We verified the efficacy of the proposed residual Monte Carlo method in a test problem defined over a domain of $x\in$[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²s in the interval [0, 5]cm and 1.0n/cm²s in the interval [6,10]cm, with a linear transition region in-between, were taken. We then assumed an S₄ 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: ${\tilde{\psi }}_i^m(x)={\tilde{\varphi }}_i(x)\left[\frac{1}{2}+\frac{{\mu }_m(x-{\overline{x}}_i)}{\Delta x_i}\right]\mathrm{,}\text{for}{x}_{i-\frac{1}{2}}<x<x_{i+\frac{1}{2}}\mathrm{,}$ (45) where ${\overline{x}}_i=\frac{x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}}{2}\mathrm{.}$ (46) 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.

6.1.2

Numerical Results and Remarks

The objective was to compute the error in the cell-integrated scalar flux $\langle \tilde{\varphi }\rangle$. We first computed the error by using the full LD SN angular flux. The results are shown in Fig. 9.

Fig. 9Full-screen View

Error estimation for cell-averaged scalar flux using the full SN angular flux.

This is a dual-y-axis plot. On the left axis, ϕMMS is the manufactured solution, ${\tilde{\varphi }}_{S_N}$ is the LD SN solution, $\langle \tilde{\varphi }\rangle$ is the cell-integrated SN scalar flux, $\delta {\langle \tilde{\varphi }\rangle }_{\text{MC}}$ is the error computed by the residual Monte Carlo method, and $\sigma :\delta {\langle \tilde{\varphi }\rangle }_{\text{MC}}$ is the standard deviation associated with the error estimation. On the right axis, $\mathcal{R}$_cell is the cell-interior residual, $\mathcal{R}$_interface is the cell-interface residual, and Q is the forward source. For this particular test, 10⁶ particles were run in ten batches, with 10⁵ 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⁶ 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 $\overline{\psi }$. The results are shown in Fig. 10, where the quantities plotted are analogous to those shown in Fig. 9.

Fig. 10Full-screen View

Error estimation for cell-averaged scalar flux using the simplified angular flux.

The errors shown in Fig. 10 are statistically identical to those found in 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 SN angular flux.

In both Fig. 9 and Fig. 10, the cell-interior residual ($\mathcal{R}$_cell) is large when the error is large, and the cell-interface residual ($\mathcal{R}$_interface) is nonzero only if a discontinuity is present in the SN 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 plots the standard deviation associated with the residual Monte Carlo calculation against those associated with the standard Monte Carlo calculations. From Fig. 11, we can see that when running with 1×10⁶ 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.

Fig. 11Full-screen View

Standard deviations for various calculations.

6.2

Boundary Source Test

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²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. A zoomed-in plot on [0,2]cm is shown in Fig. 13.

Fig. 12Full-screen View

Scalar flux error estimation in graphite with an incoming neutron source.

Fig. 13Full-screen View

Scalar flux error estimation in graphite with an incoming neutron source (zoomed-in).

In these plots, ϕSN is the piecewise linear scalar flux solution given by an S₂ calculation, ⟨ϕ⟩SN 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₂ solution.

Conclusion and Future Work

In this paper, we present a nonrecursive residual Monte Carlo method for estimating discretization error associated with the SN 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 SN error. We demonstrated the efficacy of our method with a 1-D S₄ test problem. The results show that our methods are effective at estimating the SN errors, and for the case of the angle-space integrated response, the simplified method yields exactly the same results as the full SN method. Moreover, we propose a modified simplified angular flux representation with double discontinuity, which can be employed to preserve the SN 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], 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 SN solver framework to determine whether effective parallelization within the existing infrastructure can be achieved.