Constant Component

For a constant source distribution, , the associated error is E = 0. This result was verified using our MoC 1D code, which demonstrated that the solution maintained machine precision accuracy, irrespective of the grid resolution. This outcome is expected, as q⁰ precisely represents a constant source. Consequently, the flat source approximation introduced no error. This can be easily proven as follows: (21)

Linear Component

If the source is linear in space, the FS approximation introduces an error of second order with the mesh size. Assuming , the error can be evaluated using Eq. (20). Finally, it is expanded to a Taylor series to obtain the order of accuracy (OoA). (22)where (23)which is the optical thickness seen by a neutron flying in the direction characterized by μ.

As z′ (or τ) approaches zero, namely, as we refine the spatial grid, the above error is expanded near τ = 0 as follows, (24)Therefore, the error introduced by flat source approximation is second order with respect to the mesh size.

Quadratic Component

If the source is quadratic in space, the flat source approximation will converge to the true solution to third order, which is shown below.

Assuming , the error from Eq. (20) can be evaluated, (25)As , the error approaches zero with third order, (26)However, since analytical solutions cannot be obtained directly from pure monomial sources, we instead manufactured corresponding monomial fluxes and generated the source term by substituting them into the left-hand side of Eq. (1). Owing to the presence of the differential operator in Eq. (1), a quadratic neutron flux inevitably produces a source that contains both quadratic and linear components. The inclusion of this linear component leads to a degradation in the observed order of accuracy from third to second order, as the error associated with the linear term dominates the overall error convergence.

Moreover, with induction, it can be shown that the flat source approximation is generally second order accurate in space.

General Polynomial Component

For any source that can be expressed as , the general form of polynomial sources, we can evaluate the error from Eq. (20) as follows: (27)Therefore, (28)Based on the discussions in Sects. 2.1.2 and 2.1.3, it can be inferred that the error behaves as . To prove this, we divide E(n+1) (Eq. (28)) by E(n) (Eq. (27)): (29)

These integrals can be expressed as combinations of the incomplete gamma function and the gamma function, and the behaviors of select incomplete gamma functions as are illustrated in Fig. 3: (30)Given μ and ∑t are both constant, Eq. (29) can be expressed as (31)where is the gamma function, which can be considered a constant. is the incomplete gamma function. When , it can also be viewed as a constant.

Thus Eq. (29) can be viewed as (32)Then Taylor expand this form Eq. (29) at τ=0 yields, (33)Thus the error associated with a source is one order higher than the error associated with a source . By induction, the error follows that , indicating that the flat source approximation is generally second-order accurate in space.

Fig. 3Full-screen View

Behavior of incomplete gamma functions as τ → 0

A more systematic method is presented in the following, lending itself to a more consistent way of proof among FS and LS approximations, which will be shown later.

We take the limit of the ratio of in Eq. (29) as τ approaches zero. Applying L’Hôpital’s rule yields (34)From this equation and the fact that , it follows that (35)which implies that, in general, for , an nth-order source exhibits (n+1) th-order accuracy. Consequently, the flat source approximation achieves at most second-order accuracy.

Constant Component

It is straightforward to demonstrate that the linear source approximation can represent a flat source without errors. Since q¹ vanishes and q⁰ exactly represents a constant source, the proof follows directly from Sect. 2.1.1.

Linear Component

Next, we demonstrate that the linear source approximation can exactly represent a linearly distributed source.

Assuming in q⁰ and q¹ defined in Eq. (38) gives (39)The approximating source distribution is exactly the same as the originally manufactured source distribution, as shown below (40)From Eq. (37), the error are given by (41)Thus, the linear source approximation remains an exact representation of a linear source.

Quadratic Component

Last, if the source distribution is quadratic in space as , we have (42)Constructing the linear approximation with the zeroth and first source moments gives (43)The error introduced via linear source (LS) approximation can be evaluated with Eq. (37) (44)The error approaches zero with fourth order as shown by Taylor expansion near τ = 0. (45)As mentioned in Sect. 2.1.2, we normally do not have a standalone quadratic source; rather, it usually comes with a linear component. Note that the linear component does not affect the accuracy of the linear source approximation, which can be explained by linear superposition. Since the linear source approximation represents the linear component exactly, the error only arises from the higher-order components.

To illustrate this, assume , where C is an arbitrary constant, then we have: (46)Therefore, the linear source approximating form is (47)Now the error associated with this approximation can be evaluated using Eq. (37). (48)This equation equals Eq. (44), the error associated with the source distribution . Specifically, when using the LS approximation for a quadratic source, with or without a linear component, the error is 4th order.

General Polynomial Component

Since the mathematical forms of the error in both Sect. 2.1.4 and this subsection are consistent, it is reasonable to assume that . From Sect. 2.2.2, the error can be estimated as (49)To formally prove this result, it is necessary to derive expressions for a general source. For any source expressed as , its zeroth and first spatial moments can be evaluated using Eq. (38) as follows: (50)LS approximation of presents (51)Thus, the error associated with the LS approximation can be calculated using Eq. (37) as follows: (52)E(n+1) can be easily obtained by replacing n with n+1 in the equation above. Now, we divide E(n+1) by E(n) and take the limit as . Applying L’Hôpital’s rule, i.e., taking the derivative of both the numerator and denominator, yields: (53)From this equation and , we can easily conclude (54)This result implies that, in general, for , an nth-order source exhibits (n+2) th-order accuracy. At least, the linear source approximation achieves fourth-order accuracy.

Testing Suite

The 1D MoC code we developed has tested the formal order of accuracy for FS and LS approximations. Four test cases were devised, each based on an assumed flux shape and its corresponding manufactured source, both defined in the global coordinate system.

Test listings:

Case 1: constant source distribution (55)Case 2: linear source distribution (56)Case 3: quadratic source distribution (57)Case 4: non-polynomial source distribution (58)The zeroth and first spatial moments of the MMS source are calculated using the following equation: (59)

Testing Purely Absorbing Materials

To test the distributed source only, we chose purely absorbing materials, whose scattering cross section was set to zero. Figures 4, 5, 6, 7 present the grid refinement results for four cases, each corresponding to a constant, linear, quadratic, or non-polynomial manufactured source.

Fig. 4Full-screen View

Order of accuracy with constant manufactured source (FS: exact, LS: exact)

Fig. 5Full-screen View

Order of accuracy with linear manufactured source (FS: 2nd order, LS: exact)

Fig. 6Full-screen View

Order of accuracy with quadratic manufactured source (FS: 2nd order, LS: 4th order)

Fig. 7Full-screen View

Order of accuracy with non-polynomial manufactured source (FS: 2nd order, LS: 4th order)

All experimental results align with analytical predictions, as summarized in Table 1, where the observed order of accuracy is presented alongside the formal order of accuracy in the format.

The angular errors inherent in the numerical results were isolated using an error removal technique developed in a previous study [22].

Ganapol benchmark

Benchmark Problem 3.4 from Ganapol’s analytical benchmark book [21] is a robust verification test case that can be implemented in MPACT [1] without special modifications to the code. This benchmark enables the verification of both our order-of-accuracy predictions and the accuracy of the MPACT code itself, as previously published in [15]. It also supports the broader code and solution verification activities for MPACT.

This case is presented again to confirm that the spatial convergence rate established in 1D extends to 2D problems. Below, we restate the problem configuration.

The Ganapol benchmark problem is an analytic, or “semi-analytic,” benchmark based on the exact solution of the singular integral equation that describes the single-group cylindrical transport problem. This solution methodology is a complex sequence of steps, which are described in detail in the book [21].

The benchmark considers a homogeneous right circular cylinder of infinite height (Fig. 8), characterized by:

Fig. 8Full-screen View

Infinite homogeneous cylinder

The benchmark results are provided for several cases:

(a) Uniform isotropic source - the scalar flux is tabulated for selected values of c<1;

(b) Critical rod - the critical radius is tabulated as a function of c>1;

(c) Critical rod - the scalar flux is depicted for critical rods as a function of c>1.

Table 2, sourced from Ganapol [21], provides the benchmark results for the critical rod problem, listing the critical rod radius as a function of c. The values are accurate to the last known significant digit (within eight decimal points). This table also confirms agreement with previously published benchmarks (i.e. [23, 24]), with discrepancies highlighted in bold results.

MPACT results

For MPACT verification, critical rod problems (b) in Sect. 3.2.1 are chosen as benchmarks because they exercise both the 2D MoC solver and the eigenvalue solver.

Because MPACT was used for LWR lattices, special input processing in the code is required to model the isolated cylinder. To avoid this, the benchmark is adapted by modeling the cylinder as a fuel pin within a non-scattering square bounding box (Fig. 9). On the one hand, this approximation ensures the incoming angular flux remains zero outside the rod, with scattering and fission sources contributing only within its boundaries. On the other hand, the addition of a bounding box transforms the 1D cylinder configuration into a 2D fuel-pin problem.

Fig. 9Full-screen View

The material of the cylinder

MPACT determines the eigenvalue k based on the critical rod radii listed in Table 2. Each case should yield k = 1 to high precision, verifying the accuracy of the transport solution. The radial mesh was refined by increasing the number of rings.

All cases were run with the following phase space discretization parameters:

Bounding box side length = 30 cm

Ray spacing = 0.0005 cm

Number of radial rings = [10,160] by increment of 10

Number of azimuthal slices = 32

Quadrature set = CHEBYSHEV GAUSS 32 24

Convergence criterion = 10^-7 for both k and RMS error of ϕ

Table 3 presents the eigenvalues computed by MPACT, all within a few pcm of criticality. This confirms strong agreement with analytical solutions, reinforcing the reliability of MPACT.

Mesh convergence analysis

A radial mesh refinement study was conducted using the c = 1.01 case. The convergence curve in Fig. 10 demonstrates that MPACT achieves second-order accuracy in the radial direction, which is consistent with the FS-MoC expectations.

Fig. 10Full-screen View

Convergence with respect to radial discretization