Variance Reduction

Author

Elliott Biondo, Kalin Kiesling

CADIS Method

The Consistent Adjoint-Driven Importance Sampling (CADIS) method [1] is a Monte Carlo variance reduction method that utilizes a deterministic estimate of the adjoint flux (the importance) to generate a biased source and weight windows that optimize a Monte Carlo simulation relative to a detector response function. One major feature of the scheme is “consistency”, that is, weight windows are chosen such that particles are always born within them.

In the CADIS method the response is defined as:

\[R = \int \, q(P) \, \Psi^+(P) \, dP,\]

where \(q\) is the probability distribution function describing the source strength as a function of the phase space variable \(P\) (which may represent any combination of space, energy, and direction). \(\Phi^+(P)\) is the adjoint flux relative to the detector response function being optimized. The CADIS method defines the biased source distribution as:

\[\hat{q}(P) = \frac{q(P) \, \Psi^+(P)}{R}.\]

The corresponding weight window lower bounds are defined by:

\[ww(P) = \frac{R}{\Psi^+(P) \, \frac{\beta + 1}{2}},\]

where \(\beta\) is the ratio of the weight window upper bound to the weight window lower bound (default of 5 in MCNP5).

PyNE implementation

The PyNE implementation of the CADIS method is a mesh-based implementation and is designed to be used in conjunction with the mesh-based source sampling capabilities in the source_sampling module. This means that the above method, which is continuous in phase space must be adapted for discretization of space (mesh volume elements) and energy (in energy bins).

Source density (\(q'\), units: \(time^{-1}length^{-3}\)) is the canonical quantity for representing a mesh-based source within PyNE. This means that the first step of the CADIS method within PyNE is to create a \(q\) PDF from a source density mesh. The total source strength \(q_tot\) is first found by integrating the source density over space (\(i\)) and energy (\(j\)):

\[q_{tot} = \sum_{i \in I} \, \sum_{j \in J} V_i \, q'_{i, j},\]

The \(q\) PDF can then be defined by:

\[q_{i,j } = \frac{V_i \, q'_{i, j}}{q_tot} \, for i \in I, j \in J\]

The response can then be calculated by integrating the product of \(q\) and the adjoint flux over all phase space:

\[R = \sum_{i \in I} \, \sum_{j \in J} \Psi_{i, j}^{+} \, \frac{V_i \, q'_{i, j}}{q_{tot}}\]

The weight window lower bound is then:

\[ww_{i, j} = \frac{R}{\Psi_{i, j}^{+} \, \frac{\beta + 1}{2}}.\]

These values are tagged to the weight window output mesh and can be printed out as an MCNP5 WWINP file. In the event that the adjoint flux is 0 for some \((i, j)\), the \(ww_{i, j}\) value is replaced with 0. MCNP5 will not play the weight window game when a particle enters a region of phase space where the weight window lower bound is 0.

The biased source strength is:

\[\hat{q}_{i, j} = \frac{\Psi_{i, j}^{+} \, q'_{i, j} \, V_i}{R \, q_{tot}}\]

However, the biased source strength is not the quantity of interest, because the source_sampling module is expecting biased source densities. The biased source densities that are tagged to the output mesh are:

\[\hat{q}'_{i, j} = \frac{\Psi_{i, j}^{+} \, q'_{i, j}}{R \, q_{tot}}.\]

Assumptions

The source density mesh and adjoint flux mesh must have the spatial bounds.

Sample Calculations

In this section the expected results for the the test_variancereduction.py unit test “test_cadis_multiple_e” are calculated. Consider a 2D mesh with the following properties.

\(q' = [2.6, 2.5]\)

\(\Phi = [1.3, 1.4]\)

\(V = 2\)

\(q' = [2.9, 0]\)

\(\Phi = [1.7, 1.9]\)

\(V = 2\)

\(q' = [2.9, 2.8]\)

\(\Phi = [1.1, 1.2]\)

\(V = 8\)

\(q' = [2.4, 2.2]\)

\(\Phi = [0, 1.6]\)

\(V = 8\)

Here, the vector quantities represent values at two energy groups. First calculate \(q_{tot}\) and \(R\):

\[\begin{split}q_{tot} & = \sum_{i \in I} \, \sum_{j \in J} V_i \, q'_{i, j} \\ & = 8 \cdot 2.9 + 8 \cdot 2.8 + 2 \cdot 2.6 + 2 \cdot 2.5 \\ & + 8 \cdot 2.4 + 8 \cdot 2.2 + 2 \cdot 2.9 + 2 \cdot 0 \\ & = 98.4\end{split}\]
\[\begin{split}R & = \sum_{i \in I} \, \sum_{j \in J} \Psi_{i, j}^{+} \, \frac{V_i \, q'_{i, j}}{q_{tot}} \\ & = \frac{1}{98.4} ( 1.1 \cdot 8 \cdot 2.9 + 1.2 \cdot 8 \cdot 2.8 + 1.3 \cdot 2 \cdot 2.6 + 1.4 \cdot 2 \cdot 2.5 \\ & \qquad \quad + 0 \cdot 8 \cdot 2.4 + 1.6 \cdot 8 \cdot 2.2 + 1.7 \cdot 2 \cdot 2.9 + 1.9 \cdot 2 \cdot 0 ) \\ & = 1.0587398374\end{split}\]

The expected results are:

\[\begin{split}\hat{q}' &= \frac{\Psi_{i, j}^{+} \, q'_{i, j}}{R \, q_{tot}} \\ &= \frac{1}{98.4 \cdot 1.0587398374} [1.1 \cdot 2.9, 1.2 \cdot 2.8, 1.3 \cdot 2.6, 1.4 \cdot 2.5, \\ & \qquad \qquad \qquad \qquad \qquad 0 \cdot 2.4, 1.6 \cdot 2.2, 1.7 \cdot 2.9, 1.9 \cdot 0] \\ &= [0.0306200806, 0.0322518718, 0.0324438472, 0.0335956998, \\ & \qquad 0.0, 0.0337876752, 0.0473219428, 0.0]\end{split}\]
\[\begin{split}ww &= \frac{R}{\Psi_{i, j}^{+} \, \frac{\beta + 1}{2}} \\ &= [ \frac{1.0587398374}{1.1 \cdot{3}}, \frac{1.0587398374}{1.2 \cdot{3}}, \frac{1.0587398374}{1.3 \cdot{3}}, \frac{1.0587398374}{1.4 \cdot{3}}, \\ & \qquad \frac{1.0587398374}{0 \cdot{3}}, \frac{1.0587398374}{1.6 \cdot{3}}, \frac{1.0587398374}{1.7 \cdot{3}}, \frac{1.0587398374}{1.9 \cdot{3}} ] \\ &= [0.3208302538, 0.2940943993, 0.2714717532, 0.2520809137, \\ & \qquad 0.0, 0.2205707995, 0.2075960465, 0.1857438311]\end{split}\]

Notice that the value in the \(ww\) vector that is a division by 0 has been replaced with 0.

MAGIC Method

The Method of Automatic Generation of Importances by Calculation (MAGIC) is a global variance reduction technique in which an initial particle distribution, in the form of fluxes, populations, or weights is obtained and then used to generate mesh-based weight windows or importances. This method recognizes the initial particle distribution will be poor in some highly attenuated regions but upon iteration of the MAGIC method, the solution will improve. Below are the steps for the MAGIC method. [2]

  1. Run MCNP, in analogue mode, to set up a flux meshtally. Multigroup cross section data and a high energy cut-off, corresponding to a mean-free path no greater than the mesh voxel size, should be used.

  2. Process the resulting meshtally data by normalizing the flux to have a value of 0.5 in the source (or highest) region. Use the normalized flux to create a new weight window file to be used for MCNP.

  3. Modify the original MCNP input to use the generated weight window file and run again.

  4. If results are are sufficient, no further iterations are necessary. Else, repeat starting from step 2 until desired flux results are obtained.

  5. If a high energy cut-off was used, reduce the cut-off energy and repeat iterations until the particle distribution is acceptable. The final iteration should be performed with the appropriate energy cut-off and cross section data.

PyNE implementation

The implementation of MAGIC in PyNE uses a PyNE meshtally object, which is the result of a meshtal file processed by PyNE’s mcnp.Meshtally. Using the results of the meshtal file and a specified tolerance for relative error \(t\) and null value \(\phi_o\), the flux will be normalized for each energy bin and then be used to generate a wwinp file to be used in a subsequent iteration. The steps are as follows:

  1. Read meshtally and determine maximum flux \(\phi_m^k\) for each enery bin \(k\).

  2. Normalize flux \(\phi_i^k\) in every mesh voxel \(i\) for each energy bin \(k\) according to \(\phi_m^k\) to obtain a new \(\phi_i^{'k}\). If the relative error \(e_i^k\) for voxel \(i\) in energy bin \(k\) is larger than the tolerance value \(t\), then set flux \(\phi_i^{'k}\) to the null value \(\phi_o\) instead.

\[ \begin{align}\begin{aligned}\text{If } e_i^k < t \text{ then, } \phi_i^{'k} = \frac{\phi_i^{k}}{2 \, \phi_m^k}\\\text{If } e_i^k > t \text{ then, } \phi_i^{'k} = \phi_o\end{aligned}\end{align} \]

  1. Use new flux values to create a weight window tag on the provide meshtally and use PyNE’s Wwinp class to create a weight window mesh.

Sample Calculations

In this section, the expected results of the test_variancereduction.py unit test “test_magic_multi_bins” are shown. In this test, a 3D 2x2 mesh is given. Each voxel contains flux data corresponding to 2 energy bins. The mesh is described by the following flux and relative error data.

\(\phi_1^{1} = 1.2\)

\(e_1^1 = 0.11\)

\(\phi_1^{2} = 3.3\)

\(e_1^2 = 0.013\)

\(\phi_2^{1} = 1.6\)

\(e_2^1 = 0.14\)

\(\phi_2^{2} = 1.7\)

\(e_2^2 = 0.19\)

\(\phi_3^{1} = 1.5\)

\(e_3^1 = 0.02\)

\(\phi_3^{2} = 1.4\)

\(e_3^2 = 0.16\)

\(\phi_4^{1} = 2.6\)

\(e_4^1 = 0.04\)

\(\phi_4^{2} = 1.0\)

\(e_4^2 = 0.09\)

First, the maximum flux for each energy bin is found. In this case the maximum for energy bin \(k = 1\) occurs in voxel 4 \(\phi_4^{1} = 2.6\) and in voxel 1 \(\phi_1^{2} = 3.3\) for energy bin \(k = 2\). In the first energy bin, the flux values are normalized by \(\phi_m^1 = 2.6\) and in the second \(\phi_m^2 = 3.3\). If the error tolerance is set \(t = 0.15\) and the null value set to \(\phi_o = 0.001\), then voxels 2 and 3 in the second energy bin have errors larger than the tolerance and are therefore set to the null value while everything else is normalized. The following is the result.

\(\phi_1^{'1} = \frac{\phi_1^1}{2 \, \phi_m^1} = \frac{1.2}{2*2.6} = 0.23077\)

\(\phi_1^{'2} = \frac{\phi_1^2}{2 \, \phi_m^2} = \frac{3.3}{2*3.3} = 0.5\)

\(\phi_2^{'1} = \frac{\phi_2^1}{2 \, \phi_m^1} = \frac{1.6}{2*2.6} = 0.30769\)

\(\phi_2^{'2} = \phi_o = 0.001\)

\(\phi_3^{'1} = \frac{\phi_3^1}{2 \, \phi_m^1} = \frac{1.5}{2*2.6} = 0.28846\)

\(\phi_3^{'2} = \phi_o = 0.001\)

\(\phi_4^{'1} = \frac{\phi_4^1}{2 \, \phi_m^1} = \frac{2.6}{2*2.6} = 0.5\)

\(\phi_4^{'2} = \frac{\phi_4^2}{2 \, \phi_m^2} = \frac{1.0}{2*3.3} = 0.12122\)

The values \(\phi_i^{'k}\) are then set as the new weight window values.

References

[1] Haghighat, A. and Wagner, J. C., “Monte Carlo Variance Reduction with

Deterministic Importance Functions,” Progress in Nuclear Energy, Vol. 42, No. 1, pp. 25-53, 2003.

[2] Davis, A. and Turner, A., “Comparison of global variance reduction

techniques for Monte Carlo radiation transport simulations of ITER,” Fusion Engineering and Design, Vol. 86, Issues 9-11, pp. 2698-2700, 2011.