Improving A Posteriori Spatial Discretization Error Estimation for SN Transport by Singular Characteristics Tracking

Previously, we have developed a novel spatial discretization error estimator, the “residual source” estimator, in which an error transport problem, analogous to the discretized transport equation, is solved to acquire an estimate of the error, with a residual term acting as a fixed source. Like all error estimators, the residual source estimator suffers inaccuracy and imprecision in the proximity of singular characteristics, lines across which the solution is irregular. Estimator performance worsens as the irregularities become more pronounced, especially so if the true solution itself is discontinuous. This work introduces a modification to the residual approximation procedure that seeks to reduce the adverse effects of the singular characteristics on the error estimate. A partial singular characteristic tracking scheme is implemented to reduce the portion of the error in the numerical solution born by irregularities in the true solution. This treated numerical solution informs the residual approximations. The partial singular characteristic tracking scheme greatly enhances the numerical solution for a problem with prominent singular characteristics. The residual approximation and resultant residual source error estimate are likewise improved by the scheme, which only incurs the computational cost of an extra inner iteration.


INTRODUCTION
Applying a spatial discretization scheme to the S N neutral particle transport equation incurs a departure from the spatial-continuum solution, known as the spatial discretization error. Traditionally, code users assume the spatial discretization, characterized by mesh size and method order/type, is sufficient to keep the error small. Because knowledge of the true solution predicates knowledge of the true error, there is generally no way of unequivocally assessing the extent of the validity of this assumption. However, the field of a posteriori spatial discretization error estimation seeks to compute this quantity without knowledge of the true solution, generally for the purpose of adaptive mesh refinement (AMR) or error analysis, using the numerical solution in question to construct the estimate.
While LeR exhibited good performance, all estimators were adversely affected by singular characteristics (SCs). SCs are lines (in 2D) extending from sources of irregularity (i.e., source and material discontinuities) across which the solution is irregular, and they are unavoidable in any physically realistic problem configuration. Thus, we seek an improvement on a posteriori spatial discretization error estimation for S N radiation transport problems by capturing the solution's irregularity across SCs in computing the residual with a singular characteristic tracking (SCT) algorithm.
We briefly introduce in Sec. 2 the underlying methods, including the neutron transport equation, the residual source estimator, and the singular characteristic tracking scheme. Then, we demonstrate via numerical experiments in Sec. 3 a partial SCT scheme intended to reduce the negative influence of SCs on the numerical solution, at minimal computational cost, for the purpose of approximating the residual and resultant error estimate for DGFEM-0 spatial discretization. We close with concluding remarks in Sec. 4.

Neutron Transport Equation
We consider the "true solution" to be the solution to the spatially-continuum, 2D Cartesian, onespeed, steady-state, S N transport equation with isotropic scattering [5], and explicit boundary conditions (BCs), where D = (0, X) × (0, Y ) is the spatial domain area, and ∂D is the domain boundary. Note that the assumptions and discretizations, e.g., energy, leading to this equation incur their own error vis-a-vis the full continuum transport equation, and though they are not necessarily separable from the spatial discretization error (e.g., ray effects), they are treated as such. In this work, we generate analytic solutions to Eqs. 1 and 2 using the Method of Manufactured Solutions (MMS) [6].
Equations 1 and 2 typically cannot be solved analytically, so a spatial discretization scheme is employed, DGFEM-0 in this work. The order of regularity of the true solution is highly influential in limiting the numerical solution's accuracy and global and local convergence rates [7]. Note the inability of the DGFEM-0 solution (right) to capture the discontinuities apparent in the MMS solution (left) in Fig. 1, for a problem with N x = N y = 32, X = Y = 1 cm, S 4 Level-Symmetric quadrature, σ t = 0.01 cm −1 , and σ s = 0.001 cm −1 . The reason for this choice of small values for σ t and σ s will become evident in Sec. 2.3. The solution to the spatially-discretized variant of Eqs. 1 and 2, ψ Λ h,m (x, y), where h denotes the mesh characteristic length and Λ denotes the method order, incurs an error, where Π Λ h is an L 2 projection operator onto the h-mesh, method Λ space. An approximation of ε Λ h,m (x, y), i.e., an error estimate, will be denoted with the notation Λ h,m (x, y). Standard notation for assessing an error estimate is the "effectivity", the ratio of local norms of the estimated error to the true error, where is the cell belonging to a N x × N y rectangular mesh with i = 1, . . . , N x , and j = 1, . . . , N y .

Residual Source Estimator
The residual source estimator is implicit in that the residual is used as a fixed source for an auxiliary problem, an analogous error transport problem. In operator notation, the residual is defined as where L Λ h is the discretized transport operator, and S is the discretized scattering operator. By substituting the discretized transport equation, plus BCs (necessarily zero for fixed uniform BCs, by the definition of Eq. 3). The advantage to defining the residual in the manner of Eq. 5 (sometimes referred to as the "discrete residual" or "truncation error" [8]), is that the error transport problem, despite existing in the discrete space, is exact, provided the true residual is known. Thus, we seek accurate representation of the residual.
Because the true solution is not known, the residual cannot be directly computed from Eq. 5. Thus, we have developed the Taylor expansion-based residual with approximated derivatives (TE-AD residual) formalism. By approximating the true solution with a local Taylor series expansion and inserting it into Eq. 5, the residual is approximated by functions of cell widths, problem parameters, and pointwise high-order derivatives. Because the high-order derivatives are outside the space of the numerical solution method, they also require approximation. We achieve this by representing numerical solution spatial moments in cells belonging to a "patch" surrounding the target point with the Taylor expansion. The resultant system of equations is used to acquire approximate values for the desired pointwise high-order derivatives (see [1] for a detailed description of the above process). The approximated residual is used as the fixed source in the error transport problem to attain the a posteriori estimate (LeR/TE-AD).

Singular Characteristic Tracking
SCs penalize LeR/TE-AD in two ways. First, the residual and derivative approximations are agnostic of the SCs, and the Taylor series expansion is invalid due to the discontinuity in the solution or its first derivative. Second, the diffusion of error from the SCs in the numerical solution in the region surrounding the irregularities causes inaccurate derivative approximations in cells with otherwise accurate residual approximation expressions. Furthermore, because the TE-AD acts as a fixed source, large errors in the TE-AD residual are spread in the error transport problem, causing the error estimate to suffer (this is particularly true for DGFEM-1 [2]). The first issue can be addressed by acquiring residual and derivative approximations that are not agnostic of the irregularities across the SCs. However, as long as the numerical solution is affected by diffusion of error from the SCs, the second issue remains dominant.
Because their locations are known a priori, attempts have been made to specially treat the numerical solution in cells intersected by SCs. One method is to align the mesh with the SCs, though this requires unstructured meshing and quickly becomes intractable as quadrature order and sources of discontinuity increase. Another method is Arkuszewski's "Box Explicit Method", which, when encountering a SC-intersected cell during the transport sweep, subdivides the rectangular cell three ways so the SC stretches from corner-to-corner in one cell, a special case which has enhanced convergence properties [9]. Recently, Duo developed the Singular Characteristic Tracking (SCT) algorithm, in which the cell is subdivided into regions above and below the SC, thus creating four separate subregions for computing the angular flux, Fig. 2, and separate solutions in each region are obtained with the Step Characteristics method [6,10]. The outgoing solution on either side of the SC is fed to the downstream cell that is intersected by the SC. It should be noted that the last method does not fully treat SCs, as it does not account for irregularities in the scattering source.
Our modified TE-AD approximation procedure seeks to estimate the error in a numerical solution with no treatment of SCs. However, upon meeting the specified convergence criterion, an additional sweep with Duo's SCT scheme is performed, henceforth referred to as "partial SCT" (pSCT), in an attempt to reduce the effects of the irregularities on the numerical solution without incurring the computational cost of a fully-converged, tracked solution. In our previous parametric study of error estimators, we noted that estimator performance was average in problems with low scattering ratio and a small optical thickness because of the pronounced irregularities in the true solution [1]. Incidentally, these cases also require the fewest source iterations to converge the  Fig. 3, a DGFEM-0/pSCT solution for the same problem as Fig. 1.

NUMERICAL RESULTS
To test our new method, we used the same problem described above. First we computed the true residual, enabled by the MMS, and the TE-AD approximated residual. The true residual (left) and traditionally-computed TE-AD residual with no SCT (right) are plotted in Fig. 4   The resultant LeR/TE-AD error estimates are computed and, the log 10 -scale effectivity indices are plotted over the domain in Fig. 6 (note that the TE-AD residual in a corner cell is dependent on derivatives of the flux on the boundary, which are taken from the BCs to be zero ( [1]); this is incorrect, as the expressions do not take into account the intersection of the SC in the corner cell, as previously mentioned, hence the poor agreement). For the LeR/TE-AD estimate with no SCT, Fig. 6(a), good agreement is seen away from the SCs. Around the SCs, and especially near points where multiple SCs intersect, the agreement is poor. Conversely, for the LeR/TE-AD estimate with pSCT, Fig. 6(b), while the SCs are faintly visible, agreement is generally excellent everywhere.

CONCLUSIONS
Two primary issues related to SCs have been identified in generating an LeR/TE-AD estimate: poor approximation of the TE-AD residual on SC-intersected cells due to deficiencies in the Taylor expansion representation of the true solution, and insufficiently accurate numerical solutions from spread of error from SCs resulting in poor high-order derivative approximations. By implementing a pSCT scheme, in which an additional sweep with SCT is performed after the stopping criterion is reached, and using that solution to build the TE-AD residual, the second issue has been greatly diminished, and, surprisingly, the first issue has as well. For a problem with highly pronounced SCs, the TE-AD residual and resultant LeR/TE-AD estimate have been greatly improved with little additional computational cost (the cost of a single inner iteration).
In the future, we seek to demonstrate and assess this method for other problems with varying regularity and other nuclear properties. Additionally, we have derived residual approximations for SC-intersected cells that are not deficient, and can implement them on the corner cells, using pointwise derivatives acquired by the pSCT scheme to inform the residual approximation. Because the corner cell is the first solution computed in a sweep, its error is spread to the entire problem, so it is possible that we may see even further gains by implementing this relatively simple change. Finally, we wish to extend this method to DGFEM-1 and other error estimators.