MAVRIC Appendix C: Advanced Features

This appendix contains information on several advanced features that are still under development or are non-standard use of the MAVRIC sequence.

Alternate normalization of the importance map and biased source

The importance map and biased source implemented in MAVRIC are only functions of space and energy. The importance for a specific location and energy represents the average over all directions. For applications involving a collimated beam source, a space/energy importance map may not be representative of the true importance of the particles as they stream away from the source.

As an example, consider a 14.1 MeV active interrogation beam source 1 meter from a small spherical boat containing illicit nuclear material. The objective is to compute the fission rate in the nuclear material. To create the biasing parameters, an adjoint source is located within the nuclear material and the resulting importance map is shown in Fig. 48. Note that in both the air and water, the importances change with distance from the ship, but for the beam source, the importance (to causing a fission in the nuclear material) anywhere along the beam should be the same, since there is little chance a 14.1 MeV neutron will interact with the air before striking the ship.

Fig. 48 Importance map computed using standard CADIS.

The CADIS algorithm has done exactly what it was supposed to: it made a space/energy importance map and normalized it such that the target weight where the 14.1 MeV source particles are born is 1. The problem with this is that the source particles will stream towards the ship and strike the hull where the target weight is 0.092. Since source particles have little chance of interacting in the air, the weight windows are not used to split the particle as they travel towards the ship. When source particles cross into the ship, they are split by a factor of 11 to match the target weight. For this example, splitting each particle by a factor of 11 once they strike the ship is not so bad, but for longer distances, this will result in much larger splits. For a polyenergetic source, this could lead to undersampling of the source and could result in higher variances.

To remedy this problem when using beam sources, the normalization of the importance map and biased source should not be done at the source location but instead at the point where the source particles first interact with the ship. The keyword “shiftNormPos Δx Δy Δz” will shift the source normalization position by the amounts Δx, Δy, and Δz when the biased source and importance map are developed. For the Monaco Monte Carlo calculation, the source is returned to its normal position. The source input for the above problem would then be

read sources     src 1         title="14.1 DT neutrons - collimated"         strength=1e30         sphere 0 origin x=-195 y=0 z=0   (true source position)         eDistributionID=1    (a mono-energetic 14.1 MeV distribution)         direction 1.0 0.0 0.0         dDistributionID=2    (a 2° beam )         shiftNormPos 107.7 0.0 0.0    (just inside the hull)     end src end sources

where the shift moves the source position from x = -195 to x = -87.3, just inside the hull. The resulting target weights are shown in Fig. 49 The source particles are born with weight 1 in a location with a target weight 10.9. The particle weight is not checked until the particle crosses into the hull, where the target weight is 1.0.

Fig. 49 Targets weights using the “shiftNormPos” keyword.

Other options to manipulate the importance map for special situations include the “mapMultiplier=f” keyword (in the importanceMap block or the biasing block), which will multiply every target weight by the factor f, and the keyword “noCheckAtBirth” in the parameters block will prevent the weight windows from being applied to source particles when they are started. When used in the MAVRIC sequence, the “shiftNormPos” capability automatically adds “noCheckAtBirth” to the Monaco input that is created.

Importance maps with directional information

In MAVRIC, the CADIS method is implemented in space and energy, but in general, it could also include particle direction as well. This formulation would be the following:

True source:
$q\left( \overrightarrow{r},E,\widehat{\Omega} \right)$
Desired response:
$\sigma\left( \overrightarrow{r},E,\widehat{\Omega}\right)$
Adjoint flux using $$q^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \ \sigma \left( \overrightarrow{r},E,\widehat{\Omega} \right)$$:
$\psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right)$
Estimate of detector response
(25)$R = \iiint_{}^{}{q\left( \overrightarrow{r},E,\widehat{\Omega} \right)\ \psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega}\right)}d\text{Ω } dE \ dV$
Biased source:
(26)$\widehat{q}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \frac{1}{R}q\left( \overrightarrow{r},E,\widehat{\Omega} \right)\ \psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right)$
Target weight windows:
(27)$\overline{w}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \frac{R}{\psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right)}$

For a system using a deterministic method to compute the adjoint fluxes, this completely general, space/energy/angle, approach presents many difficulties in implementation, namely,

1. dealing with the amount of memory required for a $$\left( \overrightarrow{r},E,\widehat{\Omega} \right)$$ importance map in memory,

2. interpolating the importance for particle directions in between quadrature angles, and

3. expressing the biased source in a form suitable for a general MC code since the above biased source is, in general, not separable.

Approaches incorporating directional information

Completely general space/energy/angle CADIS is most likely too difficult to implement and may not be necessary for most applications. In most real problems that involve directionally dependent source distributions, the directional dependence is azimuthally symmetric about some reference direction, $$\widehat{d}$$. The angular distribution, $$q_{i}\left( \widehat{\Omega} \right)$$, can be expressed as the product of the uniform azimuthal distribution and a polar distribution about reference direction $${\widehat{d}}_{i}$$ giving $$\frac{1}{2\pi}q_{i}\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right)$$. The geometric size of these sources tends to be small, allowing each source distribution to be expressed as the product of two separable distributions: $$q_{i}\left( \overrightarrow{r},E,\widehat{\Omega} \right) \cong q_{i}\left( \overrightarrow{r},E \right)\ q_{i}\left( \widehat{\Omega} \right)$$.

What is needed is a CADIS method that (1) can account for the importance of a particle traveling in a certain direction; (2) can be cast as a simple modification of the space/energy CADIS method using $$\overline{w}\left( \overrightarrow{r},E \right)$$ and $$\widehat{q}\left( \overrightarrow{r},E \right)$$; and (3) is simpler than the full space/angle/energy approach. This can be done starting with the approximation that the angular component of the adjoint flux $$\psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right)$$ is separable and symmetric about the average adjoint current direction $$\widehat{n}\left( \overrightarrow{r},E \right)$$, such that

$\psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right) \cong \phi^{+}\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}f\left( \widehat{\Omega} \bullet \widehat{n} \right)\text{\ .}$

This is similar to the AVATAR approach [VRUS97] but with explicitly including the azimuthal distribution so that the standard definition $$\int_{}^{}{\phi^{+}\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}f\left( \widehat{\Omega} \bullet \widehat{n} \right)\ d\widehat{\Omega}} = \phi^{+}\left( \overrightarrow{r},E \right)$$ applies. The probability distribution function $$f\left( \mu \right)$$ describing the shape of the azimuthally symmetric current at $$\left( \overrightarrow{r},E \right)$$ has the form of

$f\left( \mu \right) = \frac{\lambda e^{\text{λμ}}}{2\ \mathrm{\sinh}\left( \lambda \right)}\ ,$

with the single parameter $$\lambda\left( \overrightarrow{r},E \right)$$ determined from $$\overline{\mu}\left( \overrightarrow{r},E \right)$$, the average cosine of scatter.

From this, we can propose that weight window targets be developed that are inversely proportional to the approximation of the adjoint angular flux:

(28)$\overline{w}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \frac{2\pi\ k}{\phi^{+}\left( \overrightarrow{r},E \right) \ \ f\left( \widehat{\Omega} \bullet \widehat{n} \right)}\ ,$

where $$k$$ is the constant of proportionality that will be adjusted to make the importance map consistent with the biased source(s). Two methods will be examined here, one without and one with biasing of the source directional dependence.

For both of the methods, the SN code Denovo was modified to report not only the adjoint scalar fluxes, $$\phi^{+}\left( \overrightarrow{r},E \right)$$, but also the adjoint net currents in $$x$$, $$y$$, and $$z$$ directions: $$J_{x}\left( \overrightarrow{r},E \right)$$, $$\ J_{y}\left( \overrightarrow{r},E \right)$$, and $$J_{z}\left( \overrightarrow{r},E \right)$$. These currents are used to find $$\widehat{n}\left( \overrightarrow{r},E \right)$$ and $$\lambda\left( \overrightarrow{r},E \right)$$. The following methods have been developed so that the standard CADIS routines can be used to compute space/energy quantities of the response per unit source $$R$$, the weight window target values $$\overline{w}\left( \overrightarrow{r},E \right)$$, and biased source $$\widehat{q}\left( \overrightarrow{r},E \right)$$ with just the adjoint scalar fluxes. These quantities are then modified by the directional information.

Directionally dependent weight windows without directional source biasing

It is proposed that the biased source $$\widehat{q}\left( \overrightarrow{r},E,\widehat{\Omega} \right)$$ should be proportional to both the true source distribution and the space/energy component of the adjoint flux:

$\widehat{q}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \frac{1}{R}\left\lbrack q\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}q\left( \widehat{\Omega} \bullet \widehat{d} \right) \right\rbrack\ \phi^{+}\left( \overrightarrow{r},E \right)\ ,$

where the constant of proportionality, $$R$$, is determined by forcing $$\widehat{q}\left( \overrightarrow{r},E,\widehat{\Omega} \right)$$ to be a pdf. Since the angular component of the adjoint flux is not included, the directional distribution of the biased source will be exactly the same as the true source. Note that this approach would be exact for cases where no directional biasing could be applied – beam sources.

For multiple sources (each with a probability distribution function $$q_{i}\left( \overrightarrow{r},E \right)$$ and a strength $$S_{i}$$, giving a total source strength of $$S = \sum_{}^{}S_{i}$$), the user is required to provide one point in phase space $$\left( {\overrightarrow{r}}_{i},E_{i},{\widehat{\Omega}}_{i} \right)$$ for each source $$i$$ that is representative of that entire source where the biased source will match the target weight windows. For each source, a vector $${\widehat{n}}_{i} = \widehat{n}\left( {\overrightarrow{r}}_{i},E_{i} \right)$$ is computed using that point. For the general case of multiple sources, the biased source sampling distribution, the biased source distributions, and the weight windows are computed using

$$R_{i} = \iint_{}^{}{q_{i}\left( \overrightarrow{r},E\right)\ \phi^{+}\left( \overrightarrow{r},E \right)} dE \ dr \ \ \ \ \ \ \ \ \ \text{(estimated response from source} \ i)$$

$$\widehat{p}\left( i \right) = \frac{{S_{i}R}_{i}\ f\left( {\widehat{\Omega}}_{i} \bullet {\widehat{n}}_{i} \right)}{\sum_{}^{}{{S_{i}R}_{i}\ f\left( {\widehat{\Omega}}_{i} \bullet {\widehat{n}}_{i} \right)}} \ \ \ \ \ \ \ \text{(biased sampling of source} \ i)$$

$${\widehat{q}}_{i}\left(\overrightarrow{r},E,\widehat{\Omega} \right) \ = \ \frac{1}{R_{i}}q_{i}\left( \overrightarrow{r},E \right)\ \phi^{+}\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}q_{i}\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right) \ = \ {\widehat{q}}_{i}\left( \overrightarrow{r},E\right)\ \frac{1}{2\pi}q_{i}\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right)$$

$$\overline{w}\left( \overrightarrow{r},E,\widehat{\Omega} \right) \ \ \ = \ \ \ \frac{\sum_{}^{}{{S_{i}R}_{i}\ f\left( {\widehat{\Omega}}_{i} \bullet {\widehat{n}}_{i} \right)}}{S\phi^{+}\left( \overrightarrow{r},E \right)}\frac{1}{\ f\left( \widehat{\Omega} \bullet \widehat{n} \right)} \ \ \ = \ \ \ \frac{\sum_{}^{}{{S_{i}R}_{i}\ f\left( {\widehat{\Omega}}_{i} \bullet {\widehat{n}}_{i} \right)}}{\sum_{}^{}{S_{i}R}_{i}}\overline{w}\left( \overrightarrow{r},E \right)\frac{1}{f\left( \widehat{\Omega} \bullet \widehat{n} \right)}$$

Directionally dependent weight windows with directional source biasing

Here it is proposed that the biased source be proportional to both the true source distribution and the approximation of the adjoint angular flux. With a small geometric source, it is also assumed that there is one vector, $${\widehat{n}}_{0} = \widehat{n}\left( {\overrightarrow{r}}_{0},E_{0} \right),$$ evaluated at a specific location and energy, which represents the adjoint current direction over that source. The biased source then looks like

\begin{align}\begin{aligned} \widehat{q}\left( \overrightarrow{r},E,\widehat{\Omega} \right) & = \frac{1}{\text{Rc}} q\left( \overrightarrow{r},E,\widehat{\Omega} \right) \ \psi^{+}\left( \overrightarrow{r},E,\widehat{\Omega} \right)\\& = \frac{1}{\text{Rc}}\left\lbrack q\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}q\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right) \right\rbrack\ \left\lbrack \phi^{+}\left( \overrightarrow{r},E \right)\ \frac{1}{2\pi}\ f\left( \widehat{\Omega} \bullet {\widehat{n}}_{0} \right) \right\rbrack\ ,\end{aligned}\end{align}

where the constant $$\text{Rc}$$ is used to make $$\widehat{q}$$ a pdf. Note that if either the original source directional distribution $$q\left( \widehat{\Omega} \right)$$ or the adjoint angular flux distribution at the source is isotropic, then $$c = \frac{1}{4\pi}$$.

For the general case of multiple sources, the biased source sampling distribution, the biased source distributions and the weight windows are

$R_{i} = \iint_{}^{}{q_{i}\left( \overrightarrow{r},E\right)\ \phi^{+}\left( \overrightarrow{r},E \right)}\text{dE}\ \text{dr}$
$c_{i} = \int_{}^{}{\frac{1}{2\pi}q_{i}\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right)\ \frac{1}{2\pi}f\left( \widehat{\Omega} \bullet {\widehat{n}}_{i} \right)}d\widehat{\Omega}$
$\widehat{p}\left( i \right) = \frac{{S_{i}R}_{i}c_{i}}{\sum_{}^{}{{S_{i}R}_{i}c_{i}}}$
${\widehat{q}}_{i}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \left\lbrack \frac{1}{R_{i}}\ q_{i}\left( \overrightarrow{r},E \right)\ \phi^{+}\left( \overrightarrow{r},E \right) \right\rbrack\ \left\lbrack \frac{1}{c_{i}}\ q_{i}\left( \widehat{\Omega} \right)\ f\left( \widehat{\Omega} \right) \right\rbrack = {\widehat{q}}_{i}\left( \overrightarrow{r},E \right)\ \frac{1}{c_{i}}\ \frac{1}{2\pi}q_{i}\left( \widehat{\Omega} \bullet {\widehat{d}}_{i} \right)\ \frac{1}{2\pi}f\left( \widehat{\Omega} \bullet {\widehat{n}}_{i} \right)$
$\overline{w}\left( \overrightarrow{r},E,\widehat{\Omega} \right) = \frac{\sum_{}^{}{{S_{i}R}_{i}c_{i}}}{S\phi^{+}\left( \overrightarrow{r},E \right)} \ \ \frac{2\pi}{\ f\left( \widehat{\Omega} \bullet \widehat{n} \right)} = \frac{\sum_{}^{}{{S_{i}R}_{i}c_{i}}}{\sum_{}^{}{S_{i}R}_{i}} \ \ \overline{w}\left( \overrightarrow{r},E \right)\ \frac{2\pi}{f\left( \widehat{\Omega} \bullet \widehat{n} \right)} \ .$

More details on the development of these methods and their application for several problems have been presented [PMEW10][PME12].

Using space/energy/angle CADIS in MAVRIC

The two angular CADIS methods that use the AVATAR-type approximation of adjoint flux are specified in MAVRIC with the “angularBiasing=” keyword in the importanceMap block. Values for this keyword are 1 or 2.

Space/Energy/Angle CADIS without directional biasing (for beam sources) – This method uses one specific location, $${\overrightarrow{r}}_{0}$$, energy, $$E_{0}$$, and direction, $${\widehat{\Omega}}_{0},$$ which is the reference direction of the source $$\widehat{d}$$, where the weight of the biased source matches the weight window.

Space/Energy/Angle CADIS with directional biasing (for general sources) – This method uses one specific energy, $$E_{0}$$, to determine the adjoint current vector $${\widehat{n}}_{0}$$ and the $$\lambda_{0}$$ parameter for the biased angular distribution for each source.

With each method, the user must specify at what energy the importance map and the biased sources should be made consistent. The particle type must also be specified. This is done with the keywords “angBiasParType=” (1 for neutron or 2 for photon) and “angBiasEnergy=” (with a value in eV), also in the importanceMap block.

Note that all sources should have a direction $$\widehat{d}$$ set, using “direction u v w” within each source, even if the angular distribution for a given source is isotropic. The direction is used for source biasing and for aligning the weight windows and biased sources. Also note that for either angular biasing method, the Denovo SN calculation must use a Legendre order greater than 0.

With angular biasing, a mesh angular information (*.mai) file is produced which can be visualized with the MeshFileViewer. This file contains the space/energy-dependent $$\lambda\left( \overrightarrow{r},E \right)$$ values and components of the average adjoint current direction $$\widehat{n}\left( \overrightarrow{r},E \right) = \left\langle n_{x},n_{y},n_{z} \right\rangle$$. An existing mesh angular information (*.mai) file can be used in a separate MAVRIC problem by using the “meshAngInfoFile=” keyword in the biasing block.

Example problem

Consider the Ueki shielding problem used as sample problems in the Monaco and MAVRIC manuals. The goal is to calculate the neutron dose on one side of a shield from a partially collimated 252Cf source on the other side of the shield. Both of the angular approaches discussed above can be compared to analog and standard space/energy CADIS calculations. For the analog calculations, no importanceMap block is used. For the other cases, the importance map blocks are shown below.