## Advanced Variance Reduction for Global k-Eigenvalue Simulations in MCNP

### BibTeX

@MISC{Martin_advancedvariance,

author = {William R Martin},

title = { Advanced Variance Reduction for Global k-Eigenvalue Simulations in MCNP},

year = {}

}

### OpenURL

### Abstract

Project Objectives: • Develop advanced variance reduction techniques, based on the Variational Variance Reduction (VVR) methodology, for improving the efficiency of global reactor calculations, and incorporate these results into the production Monte Carlo code (MCNP5). • Implement this new version of MCNP5 on a Linux cluster at the University of Michigan, in order to process these histories faster than presently possible on a typical workstation. • Test this methodology for realistic problems. 1 Background: The estimation of k-eigenvalues and eigenfunctions for practical reactor and other nuclear configurations is an essential computational problem in nuclear engineering. Present-day computer codes for simulating k-eigenvalue problems are either purely deterministic or purely stochastic (Monte Carlo). Deterministic methods employ a discretization of all the independent variables in the Boltzmann equation: space, direction of flight, and energy. This discretization replaces the continuous Boltzmann equation by a (typically) large linear algebraic system of equations, which is then solved. This approach introduces truncation errors in each of the independent variables and is limited by the shape and size of the grids used for each independent variable. Alternatively, the stochastic or Monte Carlo method (i) uses pseudo-random number sequences to simulate the random histories of individual neutrons, and (ii) averages the results over many histories. Of the two approaches, Monte Carlo is often considered to be more accurate because it requires no angular or energy discretization, and in principle it can handle arbitrary geometries. However, Monte Carlo solutions have statistical errors that decrease slowly with the number of histories, i.e. with the computing time. Also, if nonanalog Monte Carlo methods are used, which is often the case, then significant user time can be required to develop -through a potentially lengthy trial-and-error process -adequate biasing parameters. Thus, while Monte Carlo solutions can be more accurate, they are often much more expensive for the user to set up and run. The goal of this project has been to (i) combine concepts from deterministic and stochastic methods to obtain an improved hybrid strategy for obtaining Monte Carlo estimates of k-eigenvalues and their corresponding eigenfunctions, and (ii) implement these methods in MCNP5. This work complements previous recent work in which deterministic and Monte Carlo methods have been merged to enhance the solution of Monte Carlo source-detector problems [1,2]. In all the hybrid methods that we have investigated for both fixed-source and eigenvalue problems, we have preserved the unbiased nature of Monte Carlo simulations. The deterministic elements of our work are designed to (i) increase the efficiency (and reduce the variance) of the Monte Carlo simulations; and (ii) retain the principle that, as the number of particle histories tends to infinity, the Monte Carlo estimate should become exact. For this project, the goal of our work is to develop and implement a new hybrid method that will perform difficult k-eigenvalue calculations more accurately and efficiently, and with less user input. Summary of Work Done for This Project: In work done prior to this project, we explored the use of a variational functional to estimate the criticality k of a fissile system [3,4]. The variational functional is more accurate than the conventional functional used in Monte Carlo k-calculations, but it depends on accurate estimates of both the forward and adjoint eigenfunctions. (The less-accurate conventional functional depends 2 only on an accurate estimate of the forward eigenfunction.) In this preliminary work, we showed that one can estimate the adjoint eigenfunction either by an initial deterministic calculation or by Monte Carlo, and we estimated the forward eigenfunction by Monte Carlo. We then discovered that the Monte Carlo estimation of forward and adjoint eigenfunctions is enhanced by a new correcton procedure in which Monte Carlo is used not to directly estimate the eigenfunction, but rather to estimate the multiplicative correction to an inexpensive deterministic estimate of the eigenfunction [5,6,7]. (For example, this deterministic estimate could be obtained from a diffusion calculation.) Because (i) the variational functional is more accurate than the conventional functional and (ii) the correcton procedure yielded more accurate estimates of the flux for fixed-source problems, we reasoned that for both reasons, this procedure should yield more accurate estimates of the eigenfunction and eigenvalue for eigenvalue problems. The down side is that (i) the variational functional for k requires estimations of both the forward and adjoint eigenfunctions (current methods for estimating eigenvalues use a simpler but less accurate functional that only requires an estimate of the forward eigenfunction) and (ii) evaluating the variational functional requires more computer algebra, and hence more time and expense. In our early work, the plusses greatly outweighed the minuses, giving a significantly more efficient hybrid method. Based on this experience, we originally proposed for this project to numerically estimate k-eigenvalues by (i) using the more accurate variational functional [3,4], and (ii) determining the forward and adjoint eigenfunctions using the more accurate correcton procedure. During the first year of the project, we developed and tested the originally-proposed correcton procedure for estimating a forward (or adjoint) eigenfunction [5,6,7]. This approach was quite successful for reactor shielding problems, and for such applications was written up and published in Nuclear Science and Engineering [7]. Unfortunately, we discovered that the correcton method does not overcome a deficiency that plagues all known Monte Carlo k-eigenvalue simulations: the unreliable estimate of the eigenfunction for "difficult" optically thick problems. This deficiency manifests itself in slightly different ways for two different types of problems: • For optically thick fissile systems (such as a commercial nuclear reactor core), it is difficult for neutrons on one side of the system to communicate with distant parts of the system. For this reason, an unphysical "tilting" or "undersampling" of the Monte Carlo estimate of the eigenfunction occurs, except when an extraordinarily large number of Monte Carlo particles per generation is used. The manifestation of this effect is that instead of acquiring a classic symmetric "cosine" shape, the eigenfunction tilts and slowly -from one generation to the next -wobbles in random and asymmetric ways around the basic cosine shape. When viewed over a sequence of many generations, this wobbling is suggestive of a flag blowing in the wind. • For optically thick systems containing numerous small but weakly-coupled fissile regions (such as spent nuclear fuel storage containers), it is difficult for neutrons in one fissile region to communicate with the other fissile regions. In this case, Monte Carlo estimates of the 3 eigenfunction tend to have the correct shapes within each fissile region, but the amplitude of the eigenfunction estimates -which depends very sensitively on the few neutrons that propagate between fissile regions -is incorrect. These relative amplitudes change slowly and randomly from one fission generation to the next. As the fissile systems become increasingly isolated, the problem for estimating the k-eigenfunction becomes more difficult -the problem becomes increasingly ill-posed. (In contrast, the single thick fissile regions discussed in the first bullet above do not become ill-posed as the system becomes thick.) The only known remedy for these deficiencies -using a sufficient number of Monte Carlo particles per generation -is impractical, due to the potentially huge number of Monte Carlo particles required. During much of the first and second years of the project, we struggled to overcome these difficulties, along with the inherent disadvantage that the variational approach requires separate eigenfunction simulations for the forward and adjoint eigenfunctions. 2. One introduces a spatial grid defined by points x j+1/2 (for 1-D geometry). This grid is used to calculate space-angle moments of the solution, but these various moments will have no truncation errors. 3. One defines functions f j+1/2 (x) and g j+1/2 (x) which are local to x j+1/2 , i.e. which are nonzero only on the interval x j−1/2 < x < x j+3/2 . One operates on the Boltzmann equation by the operator: for n = 0 and 1 and suitable choices of f j+1/2 (x), to obtain an exact system of equations for the space-angle moments of ψ: 5. The odd-n moments φ 1 j+1/2 are eliminated in terms of the even-n moments φ 0 j+1/2 and φ 2 j+1/2 . One multiplies and divides various terms in the resulting equations by to obtain a discrete system of equations for Φ j+1/2 and k, containing nonlinear functionals of the following form: 7. The above process of constructing equations for Φ j+1/2 and k introduces no errors. Therefore, if the functionals F n j+1/2 are known exactly, then the resulting algebraic system of equations determine Φ j+1/2 and k exactly. 8. In the FMC method, (i) Monte Carlo is used to estimate functionals of the form defined by Eq. (1), and then (ii) the resulting discrete system of equations is solved for Φ j+1/2 and k. Because the functionals contain statistical errors and the equations for Φ j+1/2 contain no truncation errors, the resulting estimates of Φ j+1/2 contain only statistical errors. As the statistical errors in the functionals become small, the subsequent statistical errors in Φ j+1/2 will also become small. Thus, as stated above, the FMC method does not require that accurate estimates of the eigenfunction ψ(x, µ) be obtained. Instead, the method requires that accurate estimates of functionals of the form defined in Eq. (2) be obtained. This is a much easier task because these functionals, being nonlinear and local and containing only low-order angular moments of ψ, are much less sensitive to statistical fluctuations than direct estimates of the eigenfunction ψ. 2. For optically thick systems containing isolated fissile regions, such as nuclear waste storage tanks, the FMC method yields accurate estimates of the eigenvalue, but not necessarily the eigenfunction. The reason for the relatively inaccurate estimate of the eigenfunction is that for problems of this type, the eigenfunction can be extremely sensitive to very small changes in the system, and hence the problems can be poorly-posed for the eigenfunction. To explain this last assertion, if two identical fissile spheres are separated and surrounded by an infinite absorbing region, the k-eigenfunction for this infinite system will be symmetric about the plane that symmetrically separates the two spheres. If the distance between the two spheres is large, so that it is rare for neutrons in one sphere to propagate to the other, then a very small change in the value of νΣ f in one sphere will have (i) a major change in the eigenfunction (which will no longer be symmetric but become significantly peaked around the more fissile sphere), but (ii) only a small change on the eigenvalue k. Therefore, problems of this nature will have eigenfunctions that can be very sensitive to small changes in the physical data, but eigenvalues that will be much less sensitive. In the FMC method, the small statistical fluctuations in the functionals can be viewed as corresponding to small statistical fluctuations in the underlying problem data. These small fluctuations in the functionals can yield large changes in the estimates of the eigenfunction, but only small changes in the estimate of k. This reasoning indicates that any numerical method to estimate k for these kinds of systems will at best be able to confidently estimate an accurate k. The problems are inherently illposed for the eigenfunction. Any approximation in the calculation could lead to a major error in the estimate of the eigenfunction. (Fortunately, for these problems, estimates of the eigenvalue are usually much more useful than estimates of the eigenfunction.) 3. A well-known difficulty in standard Monte Carlo estimates of k is that the estimated variance in k can be significantly smaller than the true variance. This can cause the user to erroneously believe that the calculation yields significantly more accurate estimates than it actually does. This underestimation of the variance is caused by strong correlations between the estimates 6 of the fission source in successive fission generations. In our simulations of k using the FMC method, these strong correlations in the fission source are much smaller, so it seems likely that the FMC estimates of the variance are much closer to the true variance than the standard estimates. However, we have not been able to investigate this fully. The current status of this work is described below. (Some of this work was done after the official conclusion of the project in September, 2007.) • The FMC method generalizes to multigroup problems without significant difficulty. We have extended our theory and code to these problems and have performed numerous numerical simulations. The basic advantages of the FMC method observed for one energy group extend with little variation to multiple energy groups. • We have not yet been able to extend the method to 1-D continuous-energy problems, or to energy-dependent problems with anisotropic scattering. The inclusion of continuous energy and anisotropic scattering yields extra terms in the discrete equations that are not as straightforward to treat as in the case of isotropic scattering, and it is not clear how to optimally treat them. This should not be a major hurdle, but it will require time to overcome. • We also have not been able to extend our work to multi-D. Doing this should not be conceptually difficult, but will require time. • We were not able to implement the FMC method in MCNP-5. Unfortunately, this aspect of our original proposal remains unfulfilled. • We presented our 1-D monoenergetic results at the winter 2007 ANS conference [8]. The organizers of this conference were sufficiently impressed with these results that they invited us to expand the abstract to a journal article and submit it to Nuclear Science and Engineering (NS&E), where it was quickly reviewed, and after minor revisions, accepted for publication. It is scheduled to be published in June, 2008 [9]. The ANS abstract and NS&E preprint are copied in appendices to this report. • The FMC method described above yields very accurate estimates of the eigenfunction for optically thick systems. In this sense, the FMC method is global; it inherently yields accurate estimates of the flux everywhere in the system. With very little effort, the FMC method can be adapted to fixed source problems, where it could be used (for example) to generate global solutions for deep penetration shielding calculations. For these problems, some special techniques would have to be used to ensure that a sufficient number of Monte Carlo particles exist in the "deep" regions of the system so that accurate estimates of the functional are obtained in these regions. However, this is not a difficult task. • Currently, we are seeking funding to continue this work, either for eigenvalue or fixed-source problems. We continue to believe that the concept of the FMC method has significant promise. In the following appendix, we include two publications that detail some of the accomplishments of the project. The first publication is a recent 3-page American Nuclear Society conference abstract: • "New 'Monte Carlo Functional' Methods for Estimating k-Eigenvalues and Eigenfunctions," Trans. Am. Nucl. Soc. 97, 469 (2007). The second is a to-be-published 20-page article in Nuclear Science and Engineering: • INTRODUCTION Monte Carlo methods for estimating k-eigenvalues and eigenfunctions can be problematic for optically thick fissile systems. In such problems, Monte Carlo estimates of the eigenfunctions often exhibit an unphysical "tilting," caused by (i) a large dominance ratio, and (ii) the arbitrary amplitude of the eigenfunction and the inability of neutrons in one part of the system to "see" the amplitude of the solution in distant parts. In this abstract, we propose several new Monte Carlo methods that show a strong potential for remedying this difficulty. These methods solve the "high-order" Boltzmann equation by combining it with "low-order" equations obtained by taking space-angle moments of the high-order equation. The low-order equations contain nonlinear functionals of the solution of the highorder equation; these functionals make the high-order and low-order equations consistent. The procedure used here is based on the Quasidiffusion (QD) method for S N problems [1]- [3], but the QD method differs because it uses only angular moments of the transport equation. DESCRIPTION OF THE METHODS We consider for simplicity a planar-geometry, onegroup k-eigenvalue problem with vacuum boundaries: To begin, we follow the QD procedure and operate on Eq. (1) by 1 −1 μ n (·) dμ for n = 1 and 2. Defining φ n (x) = 1 −1 μ n ψ(x, μ) dμ, we get: Solving Eq. (2b) for φ 1 (x) and introducing the result into Eq. (2a), we obtain: Also, operating on Eq. (1b) by 1 0 2μ(·) dμ, we get: Using Eq. (2b) to eliminate φ 1 (0), we obtain: A similar equation holds at x = X, which we will not write here. Eqs. (1). Now we depart from the QD approach and take certain spatial moments of Eq. (3). To simplify the algebra, we assume that the system is homogeneous and impose a uniform grid on the system with J cells, each of width h = X/J. (However, this analysis can easily be extended to inhomogeneous media and nonuniform grids.) We define x j+1/2 = jh for 0 ≤ j ≤ J to be the grid points; x j±1/2 are the edges of the j th cell. For 0 ≤ j ≤ J we define the "tent" functions: and for 1 ≤ j ≤ J − 1, 1.0 Figure 1: The Tent Functions Now for 1 ≤ j ≤ J − 1, we multiply Eq. (3) by f j+1/2 (x) and integrate over x j−1/2 < x < x j+3/2 . Carefully integrating by parts to unfold the derivatives of φ 2 (x), 469 Transport Methods: General we obtain: Also, we multiply Eq. (3) by f 1/2 (x) and integrate over 0 < x < x 1/2 . Carefully integrating by parts to unfold the derivatives of φ 2 and using Eq. (2b), we obtain: A similar equation holds at x = X. Eqs. (1). Now we define functions g j+1/2 (x) for 0 ≤ j ≤ J which are "local" about x = x j+1/2 . In this paper, we consider two such sets of functions: and (Other definitions of g j+1/2 are possible, including for example g j+1/2 = f j+1/2 .) Defining: and the nonlinear functionals: we may rewrite Eqs. with an equation similar to (11b) holding for j = J. In this work, we use Eqs. 1. We run a standard Monte Carlo simulation of Eqs. For each active cycle, we estimate k, and we obtain the standard Monte Carlo estimate of k by averaging the k for each cycle over all active cycles. 2. During each cycle, we also process the histories of the Monte Carlo particles to obtain estimates of the integrals used in Eqs. 3. We also considered an approximate method, in which rather than using Eq. (10b) to define F j+1/2 , we simply set F j+1/2 = 1. [In the system interior, this is the value given by Eq. (10b) when φ 0 (x) is a linear function of x.] Because Eqs. 4. In this work, the Monte Carlo particle histories are not affected by the low-order calculations of Eqs. (10) and (11). The standard Monte Carlo method for simulating particle histories is followed, but the results are processed as described above to obtain the new estimates of k and Φ j+1/2 . A future method could use the results for Φ j+1/2 to improve the actual Monte Carlo process for particle histories. NUMERICAL RESULTS To test the methods described above, we considered a homogeneous slab with X = 200 cm, Σ γ = 0.084 cm −1 , Σ f = 0.060 cm −1 , Σ s = 0.856 cm −1 , and ν = 2.4. We ran this (diffusive) problem with 50 inactive cycles and 200 active cycles, using 50,000 histories/cycle. The MCQD and MCF methods used a grid of h = 1.0 cm = 1.0 mfp. The estimated values of k, with their estimated standard deviations, are given in Method Est. We see that the standard Monte Carlo estimate of φ is significantly in error, and the error in the standard Monte Carlo estimate of k is much larger than the corresponding error in k for the MCQD and MCF methods. (The truncation error in the MCQD method is not apparent for this problem.) Also, the standard Monte Carlo estimate of k has a much larger estimated standard deviation than the new estimates of k. We comment that the MCQD and MCF estimates of the eigenfunction for each cycle appear to the eye to be almost as accurate as the average (over 100 active cycles) shown in These preliminary calculations show that using Monte Carlo simulations to estimate nonlinear functionals of the eigenfunction -which can be much more stable and accurate than estimates of the eigenfunction itself -can indirectly yield much more accurate estimates of the eigenfunction and eigenvalue. CONCLUSIONS We have presented new hybrid (MCQD and MCF) Monte Carlo methods for k eigenvalue/eigenfunction problems. The new methods differ from conventional Monte Carlo methods by the use of estimates of nonlinear functionals of the flux to indirectly obtain estimates of the eigenfunction and eigenvalue. (The MCQD method described here has truncation errors, but the two MCF methods do not.) For the problems tested, the nonlinear functionals are much more accurate than the direct Monte Carlo estimates of the eigenfunction; for this reason, the MCQD and MCF estimates of the eigenfunction and eigenvalue have much smaller variance than in standard Monte Carlo. In future work, we plan to extend the MCF method to realistic multi-D problems with anisotropic scattering and energy dependence. REFERENCES [1] V.Ya. Gol'din, "A Quasi-Diffusion Method for Solving the Kinetic Equation," ZH. Vych. Mat. 4, 1078