4.1. The Method of Nearby Problems for Fixed Source Problems

The method of nearby problems [RRH07] works by creating a transport problem with an analytic solution that is approximately equivalent to a numerical solution. This transport problem is “nearby” the transport problem that generated the numerical solution. The nearby problem and the analytical solution can then be compared for discretization error in a manner similar to that between the numerical estimate and an analytical or high-fidelity reference solution.

In order to apply the method of nearby problems to the neutron transport equation, the discrete ordinates multigroup neutron transport equation is taken from (1). The steady-state fixed source neutron transport equation is

(1)\[\bsOmega \cdot \nabla \psi_{m, g} + \sig[g]{t} \psi_{m, g} = \sum_{g'=1}^{G} \sig[g' \rightarrow g]{s} \phi_{g'} + \chi_{g} \sum_{g'=1}^{G} \nu_{g'} \sig[g' \rightarrow g]{f} \phi_{g'} + q_{m,g}\]

with time-independent initial and boundary conditions used to solve for the numerical angular flux \(\psi_{m,g}\). This is solved through source iteration and is explained fully in the discrete ordinates section.

Using the numerical angular flux at each spatial cell center, curve fits (\(S_{m, g}\)) are created for each energy group and angle. While quintic Hermite splines were used in the original presentation for \(C^3\) continuity [RRH07], both cubic Hermite and quintic Hermite splines are used in this analysis. Given that the solution to a transport problem can be non-smooth at material interfaces, separate spline curve fits are created for each material region. The splines at each material interface do not have to be continuous. Since the polynomial for this curve fit is known analytically, it will act as the reference solution to the nearby problem. The formulation for the cubic and quintic Hermite splines is discussed in the curve fit section.

To calculate the difference between the continuous curve fit and the numerical angular flux, a continuous residual \(R_{m,g}\) is created such that

(2)\[R_{m,g} = \bsOmega \cdot \nabla S_{m, g} + \sig[g]{t} S_{m, g} - \left( \sum_{g'=1}^{G} \sig[g' \rightarrow g]{s} \overline{S}_{g'} + \chi_{g} \sum_{g'=1}^{G} \nu_{g'} \sig[g' \rightarrow g]{f} \overline{S}_{g'} + q_{m,g} \right),\]

where \(\overline{S}_{g}\) is the scalar curve fit flux and is calculated in the same manner as the numerical scalar flux in (2). The curve fit must be integrated over each spatial cell and associated splines to calculate the discrete residual \(r_{m,g}\). The process of calculating and discretizing the residual is addressed in the residual section.

The discretized residual is then added to the original external source term \(q_{m, g}\), making the updated neutron transport equation

(3)\[\bsOmega \cdot \nabla \psi_{m, g}\mnp + \sig[g]{t} \psi_{m, g}\mnp = \sum_{g'=1}^{G} \sig[g' \rightarrow g]{s} \phi_{g'}\mnp + \chi_{g} \sum_{g'=1}^{G} \nu_{g'} \sig[g' \rightarrow g]{f} \phi_{g'}\mnp + q_{m,g} + r_{m,g}.\]

The boundary condition must also be updated to agree with the curve fit at the boundaries, such that

(4)\[\psi_{m,g} (\bx) = S_{m, g}(\bx),\]

although this means that the incoming boundary source can be negative. These equations solve for the nearby angular flux \(\psi_{m,g}\mnp\). The nearby angular flux and the curve fit flux are used to estimate the discretization error when an analytical or high-fidelity solution \(\Psi_{m,g}\) is unknown. The assertion being made is that the relative discretization error (RDE) between the exact solution and numerical solution is approximately equivalent to the RDE between the curve fit flux and nearby flux,

(5)\[\text{RDE}_{m,g} = \left| \frac{\psi_{m,g} - \Psi_{m,g}}{\Psi_{m,g}} \right| \approx \left| \frac{\psi_{m,g}\mnp - S_{m,g}}{S_{m,g}} \right| = \text{MNP}_{m,g}.\]

The method of nearby problems does not require multiple spatial grids or analytical solutions.

Previous work [WMHS23] has shown that the method of nearby problems is effective with discrete ordinates multigroup neutron transport problems for one-dimensional problems using a discrete ordinates solver. This analysis includes results for one- and two-dimensional problems using both a discrete ordinates solver and a Monte Carlo based solver, in this instance MC/DC [MVP+24]. To run the method of nearby problems with a Monte Carlo solver, the numerical flux and nearby flux calculations change. The curve fit flux and the residual calculations use the same procedure as the \(S_N\) nearby problems, taking into consideration the Monte Carlo normalization. After applying the residual as an additional source term for the Monte Carlo solver, the same error approximation in (5) is used.