Source Sampling¶

Author:	Elliott Biondo

Theory¶

Meshes can be used to represent particle source distributions for Monte Carlo radiation transport. On a mesh, source intensity is discretized spatially (into mesh volume elements) and by energy (into energy bins). In order to randomly sample these distributions to select particle birth parameters (position, energy, statistical weight) a discrete probability density function (PDF) must be created, which can be sampled with pseudo-random variates. It is convenient to create a single PDF to describe all of phase space; in other words, each bin within the PDF represents the probability that a particle is born in a particular energy group within a particular mesh volume element.

In pyne, meshes define volumetric source density $q'$ with units of $\frac{particles}{time \cdot volume}$ . In order to find the source intensity of a single phase space bin (of index $n$ ), the density must be multiplied by the volume of the mesh volume element:

$q(n) = q'(n) \cdot V(n)$

The probability $p$ that a particle is born into a particular phase space bin is given by the normalized PDF:

$p(n) = \frac{q(n)}{\sum_N{\,q(n)}}$

where $N$ is the total number of phase space bins (the number of mesh volume elements and energy groups). Phase-space bins can be selected from this PDF and all particles will have a birth weight of 1. This is known as analog sampling. Alternatively, a biased source density distribution $\hat{q}'$ can be specified yielding a biased PDF $\hat{p}(n)$ . Sampling the biased PDF requires that particles have a statistical weight:

$w(n) = \frac{p(n)}{\hat{p}(n)}$

Once a phase space bin is selected a position must be sampled uniformly within the selected mesh volume element to determine the (x, y, z) birth position, and energy must be uniformly sampled uniformly within the selected energy bin.

Implementation¶

The Sampler class reads $\hat{q}$ and optionally $\hat{q}'$ from a MOAB mesh. PDFs are created using the method described above. In order to efficiently sample these PDFs an alias table is created [1][2]. This data structure requires an $O(n^2)$ setup step, but then allows for $O(1)$ sampling. Monte Carlo radiation transport typically involves the simulation of $10^{6}$ to $10^{12}$ particles, so this expensive setup step is well-justified.

In the analog sampling mode, an alias table is created from $q$ . In the uniform and user-specified sampling modes, an alias table is created from $\hat{q}$ and birth weights are calculated for each phase space bin. In the uniform sampling mode, $\hat{q}$ is created by assigning a total source density of 1 to each mesh volume element, so that all space is sampled equally. Within each mesh volume element, a normalized PDF is created on the basis of source densities at each energy.

The method for uniformly sampling within a mesh volume element of Cartesian mesh is straightforward. A vertex of the hexahedron ( $O$ ) is chosen and three vectors are created: $\vec{x}$ , $\vec{y}$ , and $\vec{z}$ . Each vector points to an adjacent vertex (in the x, y, z, direction respectively) with a magnitude equal to the length of the edge connecting the vertex to the adjacent vertex. Three random variates are chosen ( $v_1$ , $v_2$ , $v_3$ ) in order to randomly select a position ( $P$ ) within the hexahedron:

$P = O + v_1 \cdot \vec{x} + v_2 \cdot \vec{y} + v_3 \cdot \vec{z}$

A similar method is used for uniformly sampling within a tetrahedron, as described in [3].

Assumptions¶

The Sampler class chooses the (x, y, z) position within a mesh volume element with no regard for what geometry cell it lies in. Cell rejection must be implemented within the physic-code-specific wrapper script.

Sample Calculations¶

This section provides the sample calculations to justify the results in the nosetests: test_uniform, test_bias, test_bias_spatial.

Consider a mesh with two mesh volume elements with volumes (3, 0.5). The source on the mesh has two energy groups. The source density distribution is:

$q' = ((2, 1), (9, 3))$