Iterative acceleration methods for Monte Carlo and deterministic criticality calculations

GSHMC: An efficient method for molecular simulation

Anomalous fission source convergence in a Monte Carlo criticality calculation for a weakly coupled array of two fissile material units are demonstrated. Introducing coupling coefficients among array units, it is quantitatively explained that this anomaly is caused by an insufficient restoring force to the true distribution and its large statistical uncertainty, especially, in a symmetric system. A new approach for estimating the fission source intensity ratio in an array is proposed by obtaining the eigenvector of a coupling coefficient matrix. This method also gives the uncertainty of the ratio as well as the ratio, which is available for evaluating the accuracy of the obtained ratio. The correlation between a calculated k eff and the fission source intensity ratio is formulated. It is illustrated theoretically and empirically that there is no significant correlation in a symmetric two-unit array system. In general, care should be taken that a calculated k eg may be biased by an incorrect fission source distribution, especially, in a slightly asymmetric system. A regionwise weight adjustment method is developed such that the fission source intensity ratio is forced to converge to a predetermined ratio. Using this method, a satisfactory convergence can be achieved. A larger number of neutrons per generation is recommended for a Monte Carlo criticality calculation of a weakly coupled array of units.

A multidimensional block-based adaptive mesh refinement (BAMR) method for the neutral particle transport equation with diamond and linear discontinuous spatial differencing was developed several years ago. This method was implemented in the PARallel TIme-dependent SN (PARTISN) deterministic transport code. However, the only source acceleration method available with BAMR was transport synthetic acceleration. Although the block-based adaptive mesh is orthogonal, the individual mesh cells may not be simply connected. Because of this lack of simple connectivity, development of a fully consistent diffusion synthetic acceleration (DSA) method has not been possible. This paper describes the development of a DSA method based upon an additive correction to the scalar flux iterate after a transport sweep. This DSA equation is differenced using a vertex-centered diffusion discretization that is diamond-like and may be characterized as “partially” consistent. It does not appear algebraically possible to derive a diffusion discretization that is fully consistent with diamond transport differencing on AMR meshes. The diffusion matrix is symmetric positive definite, and the DSA method is effective for most applications. This BAMR-DSA solver has been implemented and tested in two dimensions for rectangular (X-Y) and cylindrical (R-Z) geometries. As expected, results confirm that a partially consistent BAMR-DSA method will introduce instabilities for extreme cases (e.g., scattering ratios approaching 1.0 with optically thick cells), but for most realistic problems, e.g., the iron-water shielding problem, the BAMR-DSA method provides an effective acceleration method.

Active distribution network is an effective way to accommodate intermittent renewable energy and a realization of smart grid. There are distributed power resources in active distribution network and the load can be adjusted. The traditional deterministic model is no longer able to evaluate the risk accurately. To deal with the randomness and uncertainty in active distribution network, a hybrid Monte Carlo method (HMC) which has the high-dimension insensitivity of Monte Carlo method and the sampling uniformity of quasi Monte Carlo method is proposed and used to form system state matrix. Stochastic load flow calculation and risk assessment of node overvoltage are carried out in active distribution network based on HMC method. An example is taken to prove that the hybrid-Monte Carlo method can make the calculation results more close to the high precision convergence results under the same calculation times.

In criticality safety evaluation, effective multiplication factor and its corresponding source distribution are often evaluated by the Monte Carlo method. For source convergence diagnostics a kind of information entropy (Shannon entropy) has, sometimes, been used. In the last conference (2009), we proposed new type information entropy in which the source distribution is expressed on an eigen-function space in order to evaluate the convergence not only of the effective multiplication factor but also of source distribution without increasing computing resources requirements. However, there arose questions concerning our entropy. Here, we answer some of the questions : (1) how to estimate expansion coefficients, (2) how to estimate the difference between the effective multiplication factor and the corresponding maximum eigenvalue, and (3) how to estimate the deviations of expansion coefficients due to the power iteration method using source distribution calculated by Monte Carlo methods. By applying our entropy to power iteration method, more effective source convergence diagnostics method combined with the “Sandwich Method” is attained.

Deterministic neutral particle transport calculations typically involve iterative solutions of the discrete Boltzmann transport equation. The iterative convergence rate can become unacceptably slow when scattering interactions dominate the particle transport processes. An effective technique to improve convergence in such situations is the diffusion synthetic acceleration (DSA) method, which uses a particle-diffusion equation to estimate the iterative error at each iteration. When group-to-group scattering is present in the neutral particle interaction data (multigroup nuclear cross sections), and the diffusion equations become coupled across energy groups, additional cost of solving the diffusion equation can potentially outweigh the gain in convergence rate. In this article, a new decoupled DSA (D-DSA) method is developed that decouples the energy-dependent diffusion equations such that the equations for each energy group can be solved independently while maintaining the rapid convergence rate of the original DSA method. The D-DSA method makes use of a diagonalization of the group-to-group scattering operator defined for each material in a problem. The D-DSA method is derived for the first-order form of the diffusion equations. Fourier analysis and numerical experiments confirm and illustrate the rapid rate of convergence associated with the less costly D-DSA method compared to the standard fully coupled DSA approach.

New considerations about the acceleration of an iterative process to couple the radiation transport equation and the equation of matter temperature are presented in this paper. Two synthetic acceleration methods [diffusion synthetic acceleration (DSA) and transport synthetic acceleration (TSA)] have been studied and analyzed, showing its strengths and weaknesses. This study is applied to an adaptive mesh refinement (AMR) context, concluding in a better performance of DSA for coarse level resolutions with higher Sn, while TSA is better for finer level boxes and smaller Sn cases. These conclusions are applied to accelerate the resolution of an AMR problem.

The diffusion synthetic acceleration (DSA) method is an iterative procedure for obtaining rapidly convergent numerical solutions of discrete ordinates problems. General P/sub L/ acceleration methods for L>1 are first considered to determine if DSA, which is equivalent to P/sub 1/ acceleration, is optimum. After determining that DSA appears to be the optimum P/sub L/ approach, it is then demonstrated that by also using standard techniques of matrix iterative acceleration, modest improvements can be made is DSA performance.

In this paper a number of variants of the Monte-Carlo method and quasi-Monte-Carlo method are used to evaluate three view factor configurations. The key question asked is “Is the ever-popular Monte-Carlo method combined with ray tracing the best approach for evaluating view factors?”. For two of the view factor configurations a Monte-Carlo method and quasi-Monte-Carlo method based on the numerical integration of the view factor integral are much more accurate than a Monte-Carlo method and quasi-Monte-Carlo method based on ray tracing. An interesting conclusion is that when ray tracing is inefficient, i.e., the view factor is small, numerical integration works very well, exactly when ray tracing has poor accuracy.For the third view factor where the two surfaces have a common edge the view factor integral has a singularity at the common edge and ray tracing works well relative to numerical integration. A new hybrid Monte-Carlo method and quasi-Monte-Carlo method that combines the advantages of ray tracing and numerical integration for the view factor configuration with a common edge is presented. The hybrid quasi-Monte-Carlo method is shown to be the most efficient method of evaluating the view factor integral.

Deterministic and Monte Carlo methods are regularly employed to conduct lattice calculations. Monte Carlo methods can effectively model a large range of complex geometries and, compared to deterministic methods, they have the major advantage of reducing systematic errors and are computationally effective when integral quantities such as effective multiplication factor or reactivity are calculated. In contrast, deterministic methods do introduce discretization approximations but usually require shorter computation times than Monte Carlo methods when detailed flux and reaction-rate solutions are sought. This work compares the results of the deterministic code DRAGON to the Monte Carlo code Serpent in the calculation of the reactivity effects for a pressurized heavy water reactor (PHWR) lattice cell containing a 37-element, natural uranium fuel bundle with heavy water coolant and moderator. The reactivity effects are determined for changes to the coolant, moderator, and fuel temperatures and to the coolant and moderator densities for zero-burnup, mid-burnup [3750 MWd/t(U)] and discharge burnup [7500 MWd/t(U)] fuel. It is found that the overall trend in the reactivity effects calculated using DRAGON match those calculated using Serpent for the burnup cases considered. However, differences that exceed the amount attributable to statistical error have been found for some reactivity effects, particularly for perturbations to coolant and moderator density and fuel temperature.

A deterministic and Monte Carlo coupling method for PWR cavity radiation streaming calculation

A new hybrid Monte Carlo method has been developed that allows the direct simulation of detonation products for ‘real’ explosives such as HMX, TATB, PBX 9501, and PBX 9502. Chemical equilibrium composition, solid carbon, and phase separation are included. The method is called a hybrid because the solid carbon phases (bulk and/or cluster) are encorporated using analytic models, while the fluid components are explicitly included in the Monte Carlo simulation. Chemical reactions are treated by a correlated interchange of atoms between species. Solid carbon enters as virtual particles with an associated Gibbs free energy. Preliminary results are presented for HMX products including the effects of cross potentials on EOS and equilibrium composition.

It is well known that the diffusion synthetic acceleration (DSA) methods for the Sn equations become ineffective in the Fokker-Planck forward-peaked scattering limit. In response to this deficiency, Morel and Manteuffel (1991) developed an angular multigrid method for the 1-D Sn equations. This method is very effective, costing roughly twice as much as DSA per source iteration, and yielding a maximum spectral radius of approximately 0.6 in the Fokker-Planck limit. Pautz, Adams, and Morel (PAM) (1999) later generalized the angular multigrid to 2-D, but it was found that the method was unstable with sufficiently forward-peaked mappings between the angular grids. The method was stabilized via a filtering technique based on diffusion operators, but this filtering also degraded the effectiveness of the overall scheme. The spectral radius was not bounded away from unity in the Fokker-Planck limit, although the method remained more effective than DSA. The purpose of this article is to recast the multidimensional PAM angular multigrid method without the filtering as an Sn preconditioner and use it in conjunction with the Generalized Minimal RESidual (GMRES) Krylov method. The approach ensures stability and our computational results demonstrate that it is also significantly more efficient than an analogous DSA-preconditioned Krylov method.

Adaptation and implementation of the Generalized Shadow Hybrid Monte Carlo (GSHMC) method for molecular simulation at constant pressure in the NPT ensemble are discussed. The resulting method, termed NPT-GSHMC, combines Andersen barostat with GSHMC to enable molecular simulations in the environment natural for biological applications, namely, at constant pressure and constant temperature. Generalized Hybrid Monte Carlo methods are designed to maintain constant temperature and volume and extending their functionality to preserving pressure is not trivial. The theoretical formulation of NPT-GSHMC was previously introduced. Our main contribution is the implementation of this methodology in the GROMACS molecular simulation package and the evaluation of properties of NPT-GSHMC, such as accuracy, performance, effectiveness for real physical systems in comparison with well-established molecular simulation techniques. Benchmarking tests are presented and the obtained preliminary results are promising. For the first time, the generalized hybrid Monte Carlo simulations at constant pressure are available within the popular open source molecular dynamics software package.

In this paper, a simple method to improve sampling efficiency of the generalized hybrid Monte Carlo (GHMC) method is presented. Compared to the standard GHMC method, our method is found to allow us to safely increase the time increment for solving an equation-of-motion in the GHMC calculations by a factor of 4. We have demonstrated various algorithmic parameter dependence on the sampling efficiency of a hydrated alanine dipeptide.