Variance reduction

MAGIC weight windows and alias-table source sampling in Nucleide.

Nucleide provides two variance-reduction utilities built on MCNP mesh tallies: MAGIC weight-window generation and mesh-based source sampling with alias tables.

MAGIC weight windows

MAGIC converts a forward mesh tally into weight-window lower bounds. The algorithm is simple and deterministic:

For each energy group $g$ , find the maximum flux over all volume elements:

\Phi_g^{\max} = \max_{v} \Phi_{v,g}.

Then the weight-window lower bound for volume element $v$ and group $g$ is

\text{ww}_{v,g} = \begin{cases} \text{null\_value} & \text{if } \varepsilon_{v,g} > \tau, \\[6pt] \dfrac{\Phi_{v,g}}{2\,\Phi_g^{\max}} & \text{otherwise}, \end{cases}

where $\varepsilon_{v,g}$ is the relative error of the tally and $\tau$ is the error tolerance. The factor of $2$ normalizes the peak-flux voxel to a lower bound of $0.5$ ; MAGIC is intended to be iterated so the weight windows refine over successive forward runs (Cooper & Larsen 2001)^cooper-2001.

Where this lives in the code

MAGIC is implemented in crates/vr-tools/src/magic.rs. The public functions are magic() for energy-integrated totals and magic_with() for explicit array selection and parameters.

Total vs. per-group mode

Comparison of MAGIC operating modes.
Aspect	Total	Per-group
Source arrays	`total_result`, `total_rel_error`	`result`, `rel_error`
Lower bounds per voxel	1	1 per energy group
Output tag	`ww_n` / `ww_p`	`ww_n` / `ww_p`
Energy bounds tag	`n_e_upper_bounds` / `p_e_upper_bounds` (single upper bound)	`n_e_upper_bounds` / `p_e_upper_bounds` (`e_bounds[1:]`)

Parameters

MAGIC tuning parameters and their defaults.
Parameter	Default	Meaning
tolerance τ	0.5	Relative error above which a bound is nulled.
null_value	0.0	Value assigned to nulled elements.

Worked example

Consider a tally with four volume elements and total fluxes $\Phi = [1.2,\, 3.3,\, 1.6,\, 1.7]$ and relative errors all below the default tolerance of $0.5$ .

Find the group maximum: $\Phi^{\max} = 3.3$ .
Compute the scale factor: $1 / (2 \times 3.3) = 1 / 6.6$ .
The lower bounds are $\text{ww} = [1.2/6.6,\, 3.3/6.6,\, 1.6/6.6,\, 1.7/6.6] = [0.182,\, 0.5,\, 0.242,\, 0.258]$ .

If the tolerance were tightened to $0.15$ , the element with error $0.19$ would be nulled to null_value instead.

Loading chart…

Alias-table source sampling

The MeshSourceSampler uses a Walker/Vose alias table to draw particle birth voxels in $O(1)$ time per sample after an $O(n)$ construction step.

Alias-table construction

Given a discrete probability mass function $p_i$ over $n$ bins:

Normalize so that $\sum_i p_i = 1$ .
Scale each probability by $n$ , giving $n \bar{p}_i$ .
Partition bins into “small” ( $n\bar{p}_i < 1$ ) and “large” ( $n\bar{p}_i \geq 1$ ).
Pair a small bin $a$ $a$ with a large bin $g$ $g$ :
- Set $\text{prob}[a] = n\bar{p}_a$ .
- Set $\text{alias}[a] = g$ .
- Reduce the large bin’s surplus: $n\bar{p}_g \leftarrow n\bar{p}_g + n\bar{p}_a - 1$ .
Drain any remaining bins with $\text{prob} = 1$ .

Sampling

To draw one index:

Generate $r_1, r_2 \in [0, 1)$ .
Pick column $i = \lfloor n \, r_1 \rfloor$ .
Return $i$ if $r_2 < \text{prob}[i]$ ; otherwise return $\text{alias}[i]$ .

Mesh source sampling modes

MeshSourceSampler bias modes.
Mode	Sampling PDF	Birth weight
Analog	∝ flux × volume	1.0
Uniform	∝ volume	analog_pdf / uniform_pdf
User	∝ user_density × volume	analog_pdf / user_pdf

Voxel-level only

MeshSourceSampler selects a voxel uniformly within the chosen cell. Sub-voxel or fractional-cell sampling modes are not implemented.

Assumptions and limitations

MAGIC requires positive flux somewhere: if every value in an energy group is $\leq 0$ , the normalization divides by zero and Error::ZeroMaxFlux is returned.
MAGIC rejects non-finite flux or error values: NaN or infinite entries in the input arrays produce Error::NonFiniteTally rather than poisoned lower bounds.
Mesh source sampling rejects negative values: negative total_result or user_pdf entries produce Error::NegativeTally; zero values remain valid zero-probability bins.
Error comparison is strict: a relative error exactly equal to the tolerance keeps its scaled value; only errors greater than $\tau$ are nulled.
Alias tables are exact for the given PDF: the construction reproduces the input probabilities to machine precision.
Voxel resolution: birth coordinates are returned as integer voxel indices; users must sample a uniform point inside the voxel themselves if sub-voxel positions are needed.

Adjacent capabilities found in PyNE, OpenMC, and the literature that Nucleide does not implement:

MAGIC iteration: the published method iterates — each forward run refines the windows for the next — and derives upper bounds from a fixed ratio to the lower bounds; Nucleide performs a single lower-bound pass.
FW-CADIS and adjoint-based global variance reduction: hybrid deterministic/Monte Carlo methods that build weight windows from an adjoint importance solve; see Wagner et al. (2014)^wagner-2014.
Built-in generators in OpenMC: WeightWindows and WeightWindowGenerator (MAGIC-like and FW-CADIS modes, updated on the fly during transport), including an internal random-ray adjoint solver.
Source biasing: OpenMC supports energy/space/angle biasing of independent sources; PyNE’s mesh sampler additionally supports sub-voxel sampling modes. Nucleide’s MeshSourceSampler is voxel-level and energy-integrated.