Yousry Y. Azmy

Newton's Method for Solving k -Eigenvalue Problems in Neutron Diffusion Theory

Abstract We present an approach to the k-eigenvalue problem in multigroup diffusion theory based on a nonlinear treatment of the generalized eigenvalue problem. A nonlinear function is posed whose roots are equal to solutions of the k-eigenvalue problem; a Newton-Krylov method is used to find these roots. The Jacobian-vector product is found exactly or by using the Jacobian-free Newton-Krylov (JFNK) approximation. Several preconditioners for the Krylov iteration are developed. These preconditioners are based on simple approximations to the Jacobian, with one special instance being the use of power iteration as a preconditioner. Using power iteration as a preconditioner allows for the Newton-Krylov approach to heavily leverage existing power method implementations in production codes. When applied as a left preconditioner, any existing power iteration can be used to form the kernel of a JFNK solution to the k-eigenvalue problem. Numerical results generated for a suite of two-dimensional reactor benchmarks show the feasibility and computational benefits of the Newton formulation as well as examine some of the numerical difficulties potentially encountered with Newton-Krylov methods. The performance of the method is also seen to be relatively insensitive to the dominance ratio for a one-dimensional slab problem.

Newton's Method for the Computation of k-Eigenvalues in SN Transport Applications

Abstract Recently, Jacobian-Free Newton-Krylov (JFNK) methods have been used to solve the k-eigenvalue problem in diffusion and transport theories. We propose an improvement to Newton’s method (NM) for solving the k-eigenvalue problem in transport theory that avoids costly within-group iterations or iterations over energy groups. We present a formulation of the k-eigenvalue problem where a nonlinear function, whose roots are solutions of the k-eigenvalue problem, is written in terms of a generic fixed-point iteration (FPI). In this way any FPI that solves the k-eigenvalue problem can be accelerated using the Newton approach, including our improved formulation. Calculations with a one-dimensional multigroup SN transport implementation in MATLAB provide a proof of principle and show that convergence to the fundamental mode is feasible. Results generated using a three-dimensional Fortran implementation of several formulations of NM for the well-known Takeda and C5G7-MOX benchmark problems confirm the efficiency of NM for realistic k-eigenvalue calculations and highlight numerous advantages over traditional FPI.

Error comparison of diamond difference, nodal, and characteristic methods for solving multidimensional transport problems with the discrete ordinates approximation

Abstract Error norms for three variants of Larsen’s benchmark problem are evaluated using three numerical methods for solving the discrete ordinates approximation of the neutron transport equation in multidimensional Cartesian geometry. The three variants of Larsen’s test problem are concerned with the incoming flux boundary conditions: unit incoming flux on the left and bottom edges (Larsen’s configuration); unit incoming flux only on the left edge; unit incoming flux only on the bottom edge. The three methods considered are the diamond-difference (DD) method, the arbitrarily high order transport (AHOT) method of the nodal type (AHOT-N), and of the characteristic type (AHOT-C). The last two methods are employed in constant, linear, and quadratic orders of spatial approximation. The cell-wise error is computed as the difference between the cell-averaged flux computed by each method and the exact value, then the L1, L2, and L∞ error norms are calculated. The new result of this study is that while integral error norms, i.e., L1 and L2, converge to zero with mesh refinement, the cellwise L∞ norm does not. Via numerical experiments we relate this behavior to solution discontinuity across the singular characteristic. Little difference is observed between the error norm behavior of the methods in spite of the fact that AHOT-C is locally exact, suggesting that numerical diffusion across the singular characteristic is the major source of error on the global scale. Nevertheless, increasing the order of spatial approximation in AHOT methods yields higher accuracy in the integral error norms sense. In general, the characteristic methods possess a given accuracy in a larger fraction of the number of computational cells compared to nodal methods or DD.

Post-convergence automatic differentiation of iterative schemes

A new approach for performing automatic differentiation (AD) of computer codes that embody an iterative procedure, based on differentiating a single additional iteration upon achieving convergence, is described and implemented. This post-convergence automatic differentiation (PAD) technique results in better accuracy of the computed derivatives, as it eliminates part of the derivatives convergence error, and a large reduction in execution time, especially when many iterations are required to achieve convergence. In addition, it provides a way to compute derivatives of the converged solution without having to repeat the entire iterative process every time new parameters are considered. These advantages are demonstrated and the PAD technique is validated via a set of three linear and nonlinear codes used to solve neutron transport and fluid flow problems. The PAD technique reduces the execution time over direct AD by a factor of up to 30 and improves the accuracy of the derivatives by up to two orders of magnitude. The PAD technique`s biggest disadvantage lies in the necessity to compute the iterative map`s Jacobian, which for large problems can be prohibitive. Methods are discussed to alleviate this difficulty.

A posteriori error estimation for the one-dimensional arbitrarily high-order transport-nodal method

Abstract An a posteriori error analysis of the spatial approximation is developed for the one-dimensional Arbitrarily High-Order Transport-Nodal method. The error estimator, based on the L 2 norm of the error, preserves the order of the convergence of the method when the mesh size tends to zero. The proposed error estimator is decomposed into error indicators to allow the quantification of local errors. A simple test problem with isotropic scattering is solved to compare the behavior of the true errors with the estimated errors.

Recent Studies on the Asymptotic Convergence of the Spatial Discretization for Two-Dimensional Discrete Ordinates Solutions

ABSTRACT In this work, four types of quadrature schemes are used to define discrete directions in the solution of a two-dimensional fixed-source discrete ordinates problem in Cartesian geometry. Such schemes enable generating numerical results for averaged scalar fluxes over specified regions of the domain with high number (up to 105) of directions per octant. Two different nodal approaches, the ADO and AHOT-N0 methods, are utilized to obtain the numerical results of interest. The AHOT-N0 solutions on a sequence of refined meshes are then used to develop an asymptotic analysis of the spatial discretization error in order to derive a reference solution. It was more clearly observed that the spatial discretization error converges asymptotically with second order for the source region with all four quadratures employed, while for the other regions refined meshes along with tighter convergence criterion must be applied to evidence the same behavior. However, in that case, some differences among the four quadrature schemes results were found.

Comparison of Spatial Discretization Methods for Solving the S N Equations Using a Three-Dimensional Method of Manufactured Solutions Benchmark Suite with Escalating Order of Nonsmoothness

Abstract A comparison of the accuracy and computational efficiency of spatial discretization methods of the three-dimensional SN equations is conducted, including discontinuous Galerkin finite element methods, the arbitrarily high-order transport method of nodal type (AHOTN), the linear-linear method, the linear-nodal (LN) method, and the higher-order diamond difference method. For this purpose, we have developed a suite of method of manufactured solutions benchmarks that provides an exact solution of the SN equations even in the presence of scattering. Most importantly, our benchmark suite permits the user to set an arbitrary level of smoothness of the exact solution across the singular characteristics. Our study focuses on the computational efficiency of the considered spatial discretization methods. Numerical results indicate that the best-performing method depends on the norm used to measure the discretization error. We employ discrete Lp norms and integral error norms in this work. For configurations with continuous exact angular flux, high-order AHOTNs perform best under Lp error norms, while the LN method performs best when measured by integral error norms. When the angular flux is discontinuous, a new singular-characteristic tracking method for three-dimensional geometries performs best among the considered methods.

A Robust Arbitrarily High-Order Transport Method of the Characteristic Type for Unstructured Grids

Abstract A reformulation of the arbitrarily high-order transport method of the characteristic type (AHOT-C) for unstructured grids (AHOT-C-UG) is presented in this work, which resolves the previous difficulties encountered in the original formalism. A general equivalence between the arbitrary-order neutron balance and arbitrary-order characteristic equations is derived, which improves the numerical computation of the spatial moments of the angular flux and allows a series expansion of the characteristic integral kernel in cases where the medium is optically thin. Numerical results are presented, which verify the convergence behavior of AHOT-C-UG for various expansion orders.

Spatial Convergence Study of Discrete Ordinates Methods Via the Singular Characteristic Tracking Algorithm

Abstract This paper analyzes the spatial discretization of the discrete ordinates (DO) approximation of the transport equation. A new method, the singular characteristics tracking algorithm, is developed to account for potential nonsmoothness across the singular characteristics in the exact solution of the DO approximation to the transport equation. Numerical results in two-dimensional problems show improved rate of convergence of the exact solution of the DO equations in nonscattering and isotropic scattering media. Unlike the standard weighted diamond difference scheme, the new algorithm achieves local convergence in the case of discontinuous angular flux across the singular characteristics. The method also significantly reduces the error for problems where the angular flux presents discontinuous spatial derivatives across these lines. For purposes of testing the performance of the new method, the method of manufactured solutions is used to generate analytical reference solutions that permit accurate estimation of the local error in case of discontinuous flux.

Properties of the SN-Equivalent Integral Transport Operator in Slab Geometry and the Iterative Acceleration of Neutral Particle Transport Methods

Abstract General expressions for the matrix elements of the discrete SN-equivalent integral transport operator are derived in slab geometry. Their asymptotic behavior versus cell optical thickness is investigated both for a homogeneous slab and for a heterogeneous slab characterized by a periodic material discontinuity wherein each optically thick cell is surrounded by two optically thin cells in a repeating pattern. In the case of a homogeneous slab, the asymptotic analysis conducted in the thick-cell limit for a highly scattering medium shows that the discretized integral transport operator approaches a tridiagonal matrix possessing a diffusion-like coupling stencil. It is further shown that this structure is approached at a fast exponential rate with increasing cell thickness when the arbitrarily high order transport method of the nodal type and zero-order spatial approximation (AHOT-N0) formalism is employed to effect the spatial discretization of the discrete ordinates transport operator. In the case of periodically heterogeneous slab configurations, the asymptotic behavior is realized by pushing apart the cells’ optical thicknesses; i.e., the thick cells are made thicker while the thin cells are made thinner at a prescribed rate. We show that in this limit the discretized integral transport operator is approximated by a pentadiagonal structure. Notwithstanding, the discrete operator is amenable to algebraic transformations leading to a matrix representation still asymptotically approaching a tridiagonal structure at a fast exponential rate bearing close resemblance to the diffusive operator. The results of the asymptotic analysis of the integral transport matrix are then used to gain insight into the excellent convergence properties of the adjacent-cell preconditioner (AP) acceleration scheme. Specifically, the AP operator exactly captures the asymptotic structure acquired by the integral transport matrix in the thick-cell limit for homogeneous slabs of pure-scatterer or partial-scatterer material, and for periodically heterogeneous slabs hosting purely scattering materials. In the above limits the integral transport matrix reduces to a diffusive structure consistent with the diffusive matrix template used to construct the AP. In the case of periodically heterogeneous slabs containing absorbing materials, the AP operator partially captures the asymptotic structure acquired by the integral transport matrix. The inexact agreement is due either to discrepancies in the equations for the boundary cells or to the nondiffusive structure acquired by the integral transport matrix. These findings shed light on the immediate convergence, i.e., convergence in two iterations, displayed by the AP acceleration scheme in the asymptotic limit for slabs hosting purely scattering materials, both in the homogeneous and periodically heterogeneous cases. For periodically heterogeneous slabs containing absorbing materials, immediate convergence is achieved by modifying the original recipe for constructing the AP so that the correct asymptotic structure of the integral transport matrix coincides with the AP operator in the asymptotic limit.