crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

adding back files to be reviewed

Browse source
This commit is contained in:
Paul Romano 2019-10-28 11:55:45 -05:00
parent ae28233110
commit bc09d1ef55
1244 changed files with 301904 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

569
docs/source/methods/cmfd.rst Normal file
View file

@ -0,0 +1,569 @@
.. _methods_cmfd:
================================================================
Nonlinear Diffusion Acceleration - Coarse Mesh Finite Difference
================================================================
This page section discusses how nonlinear diffusion acceleration (NDA) using
coarse mesh finite difference (CMFD) is implemented into OpenMC. Before we get
into the theory, general notation for this section is discussed.
Note that the methods discussed in this section are written specifically for
continuous-energy mode but equivalent apply to the multi-group mode if the
particle's energy is replaced with the particle's group
--------
Notation
--------
Before deriving NDA relationships, notation is explained. If a parameter has a
:math:`\overline{\cdot}`, it is surface area-averaged and if it has a
:math:`\overline{\overline\cdot}`, it is volume-averaged. When describing a
specific cell in the geometry, indices :math:`(i,j,k)` are used which correspond
to directions :math:`(x,y,z)`. In most cases, the same operation is performed in
all three directions. To compactly write this, an arbitrary direction set
:math:`(u,v,w)` that corresponds to cell indices :math:`(l,m,n)` is used. Note
that :math:`u` and :math:`l` do not have to correspond to :math:`x` and
:math:`i`. However, if :math:`u` and :math:`l` correspond to :math:`y` and
:math:`j`, :math:`v` and :math:`w` correspond to :math:`x` and :math:`z`
directions. An example of this is shown in the following expression:
.. math::
:label: not1
\sum\limits_{u\in(x,y,z)}\left\langle\overline{J}^{u,g}_{l+1/2,m,n}
\Delta_m^v\Delta_n^w\right\rangle
Here, :math:`u` takes on each direction one at a time. The parameter :math:`J`
is surface area-averaged over the transverse indices :math:`m` and :math:`n`
located at :math:`l+1/2`. Usually, spatial indices are listed as subscripts and
the direction as a superscript. Energy group indices represented by :math:`g`
and :math:`h` are also listed as superscripts here. The group :math:`g` is the
group of interest and, if present, :math:`h` is all groups. Finally, any
parameter surrounded by :math:`\left\langle\cdot\right\rangle` represents a
tally quantity that can be edited from a Monte Carlo (MC) solution.
------
Theory
------
NDA is a diffusion model that has equivalent physics to a transport model. There
are many different methods that can be classified as NDA. The CMFD method is a
type of NDA that represents second order multigroup diffusion equations on a
coarse spatial mesh. Whether a transport model or diffusion model is used to
represent the distribution of neutrons, these models must satisfy the *neutron
balance equation*. This balance is represented by the following formula for a
specific energy group :math:`g` in cell :math:`(l,m,n)`:
.. math::
:label: eq_neut_bal
\sum\limits_{u\in(x,y,z)}\left(\left\langle\overline{J}^{u,g}_{l+1/2,m,n}
\Delta_m^v\Delta_n^w\right\rangle -
\left\langle\overline{J}^{u,g}_{l-1/2,m,n}
\Delta_m^v\Delta_n^w\right\rangle\right)
+
\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle
= \\
\sum\limits_{h=1}^G\left\langle
\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow
g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w
\right\rangle
+
\frac{1}{k_{eff}}\sum\limits_{h=1}^G
\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow
g}\overline{\overline\phi}_{l,m,n}^h
\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle.
In eq. :eq:`eq_neut_bal` the parameters are defined as:
* :math:`\left\langle\overline{J}^{u,g}_{l\pm
1/2,m,n}\Delta_m^v\Delta_n^w\right\rangle` --- surface area-integrated net
current over surface :math:`(l\pm 1/2,m,n)` with surface normal in direction
:math:`u` in energy group :math:`g`. By dividing this quantity by the transverse
area, :math:`\Delta_m^v\Delta_n^w`, the surface area-averaged net current can
be computed.
* :math:`\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle`
--- volume-integrated total reaction rate over energy group :math:`g`.
* :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow
g}
\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle`
--- volume-integrated scattering production rate of neutrons that begin with
energy in group :math:`h` and exit reaction in group :math:`g`. This reaction
rate also includes the energy transfer of reactions (except fission) that
produce multiple neutrons such as (n, 2n); hence, the need for :math:`\nu_s`
to represent neutron multiplicity.
* :math:`k_{eff}` --- core multiplication factor.
* :math:`\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow
g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle`
--- volume-integrated fission production rate of neutrons from fissions in
group :math:`h` that exit in group :math:`g`.
Each quantity in :math:`\left\langle\cdot\right\rangle` represents a scalar value that
is obtained from an MC tally. A good verification step when using an MC code is
to make sure that tallies satisfy this balance equation within statistics. No
NDA acceleration can be performed if the balance equation is not satisfied.
There are three major steps to consider when performing NDA: (1) calculation of
macroscopic cross sections and nonlinear parameters, (2) solving an eigenvalue
problem with a system of linear equations, and (3) modifying MC source
distribution to align with the NDA solution on a chosen mesh. This process is
illustrated as a flow chart below. After a batch of neutrons
is simulated, NDA can take place. Each of the steps described above is described
in detail in the following sections.
.. figure:: ../_images/cmfd_flow.png
:align: center
:figclass: align-center
Flow chart of NDA process. Note "XS" is used for cross section and "DC" is
used for diffusion coefficient.
Calculation of Macroscopic Cross Sections
-----------------------------------------
A diffusion model needs macroscopic cross sections and diffusion coefficients to
solve for multigroup fluxes. Cross sections are derived by conserving reaction
rates predicted by MC tallies. From Eq. :eq:`eq_neut_bal`, total, scattering
production and fission production macroscopic cross sections are needed. They are
defined from MC tallies as follows:
.. math::
:label: xs1
\overline{\overline\Sigma}_{t_{l,m,n}}^g \equiv
\frac{\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}
{\left\langle\overline{\overline\phi}_{l,m,n}^g
\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle},
.. math::
:label: xs2
\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow g} \equiv
\frac{\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow
g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}
{\left\langle\overline{\overline\phi}_{l,m,n}^h
\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}
and
.. math::
:label: xs3
\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow g} \equiv
\frac{\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow
g}\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}
{\left\langle\overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}.
In order to fully conserve neutron balance, leakage rates also need to be
preserved. In standard diffusion theory, leakage rates are represented by
diffusion coefficients. Unfortunately, it is not easy in MC to calculate a
single diffusion coefficient for a cell that describes leakage out of each
surface. Luckily, it does not matter what definition of diffusion coefficient is
used because nonlinear equivalence parameters will correct for this
inconsistency. However, depending on the diffusion coefficient definition
chosen, different convergence properties of NDA equations are observed.
Here, we introduce a diffusion coefficient that is derived for a coarse energy
transport reaction rate. This definition can easily be constructed from
MC tallies provided that angular moments of scattering reaction rates can
be obtained. The diffusion coefficient is defined as follows:
.. math::
:label: eq_transD
\overline{\overline D}_{l,m,n}^g =
\frac{\left\langle\overline{\overline\phi}_{l,m,n}^g
\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle}{3
\left\langle\overline{\overline\Sigma}_{tr_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g
\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle},
where
.. math::
:label: xs4
\left\langle\overline{\overline\Sigma}_{tr_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle
=
\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle
\\ -
\left\langle\overline{\overline{\nu_s\Sigma}}_{s1_{l,m,n}}^g
\overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle.
Note that the transport reaction rate is calculated from the total reaction rate
reduced by the :math:`P_1` scattering production reaction rate. Equation :eq:`eq_transD`
does not represent the best definition of diffusion coefficients from MC;
however, it is very simple and usually fits into MC tally frameworks
easily. Different methods to calculate more accurate diffusion coefficients can
found in [Herman]_.
CMFD Equations
--------------
The first part of this section is devoted to discussing second-order finite
volume discretization of multigroup diffusion equations. This will be followed
up by the formulation of CMFD equations that are used in this NDA
scheme. When performing second-order finite volume discretization of the
diffusion equation, we need information that relates current to flux. In this
numerical scheme, each cell is coupled only to its direct neighbors. Therefore,
only two types of coupling exist: (1) cell-to-cell coupling and (2)
cell-to-boundary coupling. The derivation of this procedure is referred to as
finite difference diffusion equations and can be found in literature such
as [Hebert]_. These current/flux relationships are as follows:
* cell-to-cell coupling
.. math::
:label: eq_cell_cell
\overline{J}^{u,g}_{l\pm1/2,m,n} = -\frac{2\overline{\overline
D}_{l\pm1,m,n}^g\overline{\overline
D}_{l,m,n}^g}{\overline{\overline D}_{l\pm1,m,n}^g\Delta_l^u +
\overline{\overline
D}_{l,m,n}^g\Delta_{l\pm1}^u}
\left(\pm\overline{\overline{\phi}}_{l\pm1,m,n}^g\mp
\overline{\overline{\phi}}_{l,m,n}^g\right),
* cell-to-boundary coupling
.. math::
:label: eq_cell_bound
\overline{J}^{u,g}_{l\pm1/2,m,n} = \pm\frac{2\overline{\overline
D}_{l,m,n}^g\left(1 -
\beta_{l\pm1/2,m,n}^{u,g}\right)}{4\overline{\overline
D}_{l,m,n}^g\left(1 + \beta_{l\pm1/2,m,n}^{u,g}\right) + \left(1 -
\beta_{l\pm1/2,m,n}^{u,g}\right)\Delta_l^u}\overline{\overline{\phi}}_{l,m,n}^{g}.
In Eqs. :eq:`eq_cell_cell` and :eq:`eq_cell_bound`, the :math:`\pm` refers to
left (:math:`-x`) or right (:math:`+x`) surface in the :math:`x` direction,
back (:math:`-y`) or front (:math:`+y`) surface in the :math:`y` direction and
bottom (:math:`-z`) or top (:math:`+z`) surface in the :math:`z` direction. For
cell-to-boundary coupling, a general albedo, :math:`\beta_{l\pm1/2,m,n}^{u,g}`,
is used. The albedo is defined as the ratio of incoming (:math:`-` superscript)
to outgoing (:math:`+` superscript) partial current on any surface represented
as
.. math::
:label: eq_albedo
\beta_{l\pm1/2,m,n}^{u,g} =
\frac{\overline{J}^{u,g-}_{l\pm1/2,m,n}}{\overline{J}^{u,g+}_{l\pm1/2,m,n}}.
Common boundary conditions are: vacuum (:math:`\beta=0`), reflective
(:math:`\beta=1`) and zero flux (:math:`\beta=-1`). Both eq. :eq:`eq_cell_cell`
and eq. :eq:`eq_cell_bound` can be written in this generic form,
.. math::
:label: eq_dtilde
\overline{J}^{u,g}_{l\pm1/2,m,n} = \widetilde{D}_{l,m,n}^{u,g} \left(\dots\right).
The parameter :math:`\widetilde{D}_{l,m,n}^{u,g}` represents the linear
coupling term between current and flux. These current relationships can be
sustituted into eq. :eq:`eq_neut_bal` to produce a linear system of multigroup
diffusion equations for each spatial cell and energy group. However, a solution
to these equations is not consistent with a higher order transport solution
unless equivalence factors are present. This is because both the diffusion
approximation, governed by Fick's Law, and spatial trunction error will produce
differences. Therefore, a nonlinear parameter,
:math:`\widehat{D}_{l,m,n}^{u,g}`, is added to eqs. :eq:`eq_cell_cell` and
:eq:`eq_cell_bound`. These equations are, respectively,
.. math::
:label: eq_dhat_cell
\overline{J}^{u,g}_{l\pm1/2,m,n} = -\widetilde{D}_{l,m,n}^{u,g}
\left(\pm\overline{\overline{\phi}}_{l\pm1,m,n}^g\mp
\overline{\overline{\phi}}_{l,m,n}^g\right) + \widehat{D}_{l,m,n}^{u,g}
\left(\overline{\overline{\phi}}_{l\pm1,m,n}^g +
\overline{\overline{\phi}}_{l,m,n}^g\right)
and
.. math::
:label: eq_dhat_bound
\overline{J}^{u,g}_{l\pm1/2,m,n} = \pm\widetilde{D}_{l,m,n}^{u,g}
\overline{\overline{\phi}}_{l,m,n}^{g} + \widehat{D}_{l,m,n}^{u,g}
\overline{\overline{\phi}}_{l,m,n}^{g}.
The only unknown in each of these equations is the equivalence parameter. The
current, linear coupling term and flux can either be obtained or derived from
MC tallies. Thus, it is called nonlinear because it is dependent on the flux
which is updated on the next iteration.
Equations :eq:`eq_dhat_cell` and :eq:`eq_dhat_bound` can be substituted into
eq. :eq:`eq_neut_bal` to create a linear system of equations that is consistent
with transport physics. One example of this equation is written for an
interior cell,
.. math::
:label: eq_cmfd_sys
\sum_{u\in
x,y,x}\frac{1}{\Delta_l^u}\left[\left(-\tilde{D}_{l-1/2,m,n}^{u,g} -
\hat{D}_{l-1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l-1,m,n}^g\right.
\\ + \left(\tilde{D}_{l-1/2,m,n}^{u,g} +
\tilde{D}_{l+1/2,m,n}^{u,g} - \hat{D}_{l-1/2,m,n}^{u,g} +
\hat{D}_{l+1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l,m,n}^g
\\ +
\left. \left(-\tilde{D}_{l+1/2,m,n}^{u,g} +
\hat{D}_{l+1/2,m,n}^{u,g}\right)\overline{\overline{\phi}}_{l+1,m,n}^g
\right] \\ +
\overline{\overline\Sigma}_{t_{l,m,n}}^g\overline{\overline{\phi}}_{l,m,n}^g
- \sum\limits_{h=1}^G\overline{\overline{\nu_s\Sigma}}^{h\rightarrow
g}_{s_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h =
\frac{1}{k}\sum\limits_{h=1}^G\overline{\overline{\nu_f\Sigma}}^{h\rightarrow
g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h.
It should be noted that before substitution, eq. :eq:`eq_neut_bal` was divided
by the volume of the cell, :math:`\Delta_l^u\Delta_m^v\Delta_n^w`. Equation
:eq:`eq_cmfd_sys` can be represented in operator form as
.. math::
:label: eq_CMFDopers
\mathbb{M}\mathbf{\Phi} = \frac{1}{k}\mathbb{F}\mathbf{\Phi},
where :math:`\mathbb{M}` is the neutron loss matrix operator,
:math:`\mathbb{F}` is the neutron production matrix operator,
:math:`\mathbf{\Phi}` is the multigroup flux vector and :math:`k` is the
eigenvalue. This generalized eigenvalue problem is solved to obtain fundamental
mode multigroup fluxes and eigenvalue. In order to produce consistent results
with transport theory from these equations, the neutron balance equation must
have been satisfied by MC tallies. The desire is that CMFD equations will
produce a more accurate source than MC after each fission source generation.
CMFD Feedback
-------------
Now that a more accurate representation of the expected source distribution is
estimated from CMFD, it needs to be communicated back to MC. The first step
in this process is to generate a probability mass function that provides
information about how probable it is for a neutron to be born in a given cell
and energy group. This is represented as
.. math::
:label: eq_cmfd_psrc
p_{l,m,n}^g =
\frac{\sum_{h=1}^{G}\overline{\overline{\nu_f\Sigma}}^{h\rightarrow
g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h\Delta_l^u\Delta_m^v
\Delta_n^w}{\sum_n\sum_m\sum_l\sum_{h=1}^{G}\overline{
\overline{\nu_f\Sigma}}^{h\rightarrow
g}_{f_{l,m,n}}\overline{\overline{\phi}}_{l,m,n}^h\Delta_l^u\Delta_m^v
\Delta_n^w}.
This equation can be multiplied by the number of source neutrons to obtain an
estimate of the expected number of neutrons to be born in a given cell and
energy group. This distribution can be compared to the MC source distribution
to generate weight adjusted factors defined as
.. math::
:label: eq_waf
f_{l,m,n}^g = \frac{Np_{l,m,n}^g}{\sum\limits_s w_s};\quad s\in
\left(g,l,m,n\right).
The MC source distribution is represented on the same coarse mesh as
CMFD by summing all neutrons' weights, :math:`w_s`, in a given cell and
energy group. MC source weights can then be modified by this weight
adjustment factor so that it matches the CMFD solution on the coarse
mesh,
.. math::
:label: src_mod
w^\prime_s = w_s\times f_{l,m,n}^g;\quad s\in \left(g,l,m,n\right).
It should be noted that heterogeneous information about local coordinates and
energy remain constant throughout this modification process.
------------------------
Implementation in OpenMC
------------------------
The section describes how CMFD was implemented in OpenMC. Before the simulation
begins, a user sets up a CMFD input file that contains the following basic
information:
* CMFD mesh (space and energy),
* boundary conditions at edge of mesh (albedos),
* acceleration region (subset of mesh, optional),
* fission source generation (FSG)/batch that CMFD should begin, and
* whether CMFD feedback should be applied.
It should be noted that for more difficult simulations (e.g., light water
reactors), there are other options available to users such as tally resetting
parameters, effective down-scatter usage, tally estimator, etc. For more
information please see the :class:`openmc.cmfd.CMFDRun` class.
Of the options described above, the optional acceleration subset region is an
uncommon feature. Because OpenMC only has a structured Cartesian mesh, mesh
cells may overlay regions that don't contain fissionable material and may be so
far from the core that the neutron flux is very low. If these regions were
included in the CMFD solution, bad estimates of diffusion parameters may result
and affect CMFD feedback. To deal with this, a user can carve out an active
acceleration region from their structured Cartesian mesh. This is illustrated
in diagram below. When placing a CMFD mesh over a geometry, the boundary
conditions must be known at the global edges of the mesh. If the geometry is
complex like the one below, one may have to cover the whole geometry including
the reactor pressure vessel because we know that there is a zero incoming
current boundary condition at the outer edge of the pressure vessel. This is
not viable in practice because neutrons in simulations may not reach mesh cells
that are near the pressure vessel. To circumvent this, one can shrink the mesh
to cover just the core region as shown in the diagram. However, one must still
estimate the boundary conditions at the global boundaries, but at these
locations, they are not readily known. In OpenMC, one can carve out the active
core region from the entire structured Cartesian mesh. This is shown in the
diagram below by the darkened region over the core. The albedo boundary
conditions at the active core/reflector boundary can be tallied indirectly
during the MC simulation with incoming and outgoing partial currents. This
allows the user to not have to worry about neutrons producing adequate tallies
in mesh cells far away from the core.
.. figure:: ../_images/meshfig.png
:align: center
:figclass: align-center
Diagram of CMFD acceleration mesh
During an MC simulation, CMFD tallies are accumulated. The basic tallies needed
are listed in Table :ref:`tab_tally`. Each tally is performed on a spatial and
energy mesh basis. The surface area-integrated net current is tallied on every
surface of the mesh. OpenMC tally objects are created by the CMFD code
internally, and cross sections are calculated at each CMFD feedback iteration.
The first CMFD iteration, controlled by the user, occurs just after tallies are
communicated to the master processor. Once tallies are collapsed, cross
sections, diffusion coefficients and equivalence parameters are calculated. This
is performed only on the acceleration region if that option has been activated
by the user. Once all diffusion parameters are calculated, CMFD matrices are
formed where energy groups are the inner most iteration index. In OpenMC,
compressed row storage sparse matrices are used due to the sparsity of CMFD
operators. An example of this sparsity is shown for the 3-D BEAVRS model in
figures :num:`fig-loss` and :num:`fig-prod` [BEAVRS]_. These matrices represent
an assembly radial mesh, 24 cell mesh in the axial direction and two energy
groups. The loss matrix is 99.92% sparse and the production matrix is 99.99%
sparse. Although the loss matrix looks like it is tridiagonal, it is really a
seven banded matrix with a block diagonal matrix for scattering. The production
matrix is a :math:`2\times 2` block diagonal; however, zeros are present because
no fission neutrons appear with energies in the thermal group.
.. _tab_tally:
.. table:: OpenMC CMFD tally list
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| tally | score | filter |
+============================================================================================+================+===========================+
| \ :math:`\left\langle\overline{\overline\phi}_{l,m,n}^g | flux | mesh, energy |
| \Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| \ :math:`\left\langle\overline{\overline\Sigma}_{t_{l,m,n}}^g | total | mesh, energy |
| \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| \ :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s1_{l,m,n}}^g | nu-scatter-1 | mesh, energy |
| \overline{\overline\phi}_{l,m,n}^g\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| \ :math:`\left\langle\overline{\overline{\nu_s\Sigma}}_{s_{l,m,n}}^{h\rightarrow g} | nu-scatter | mesh, energy, energyout |
| \overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| \ :math:`\left\langle\overline{\overline{\nu_f\Sigma}}_{f_{l,m,n}}^{h\rightarrow g} | nu-fission | mesh, energy, energyout |
| \overline{\overline\phi}_{l,m,n}^h\Delta_l^u\Delta_m^v\Delta_n^w\right\rangle` | | |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
| \ :math:`\left\langle\overline{J}^{u,g}_{l\pm 1/2,m,n}\Delta_m^v\Delta_n^w\right\rangle` | current | mesh, energy |
+--------------------------------------------------------------------------------------------+----------------+---------------------------+
.. _fig-loss:
.. figure:: ../_images/loss.png
:scale: 50
Sparsity of Neutron Loss Operator
.. _fig-prod:
.. figure:: ../_images/prod.png
:scale: 50
Sparsity of Neutron Production Operator
To solve the eigenvalue problem with these matrices, different source iteration
and linear solvers can be used. The most common source iteration solver used is
standard power iteration as described in [Gill]_. To accelerate these source
iterations, a Wielandt shift scheme can be used as discussed in [Park]_. PETSc
solvers were first implemented to perform the linear solution in parallel that
occurs once per source iteration. When using PETSc, different types of parallel
linear solvers and preconditioners can be used. By default, OpenMC uses an
incomplete LU preconditioner and a GMRES Krylov solver. After some initial
studies of parallelization with PETSc, it was observed that because CMFD
matrices are very sparse, solution times do not scale well. An additional
Gauss-Seidel linear solver with Chebyshev acceleration was added that is
similar to the one used for CMFD in CASMO [Rhodes]_ and [Smith]_. This solver
was implemented with a custom section for two energy groups. Because energy
group is the inner most index, a block diagonal is formed when using more than
one group. For two groups, it is easy to invert this diagonal analytically
inside the Gauss-Seidel iterative solver. For more than two groups, this
analytic inversion can still be performed, but with more computational effort.
A standard Gauss-Seidel solver is used for more than two groups.
Besides a power iteration, a Jacobian-free Newton-Krylov method was also
implemented to obtain eigenvalue and multigroup fluxes as described in [Gill]_
and [Knoll]_. This method is not the primary one used, but has gotten recent
attention due to its coupling advantages to other physics such as thermal
hydraulics. Once multigroup fluxes are obtained, a normalized fission source is
calculated in the code using eq. :eq:`eq_cmfd_psrc` directly.
The next step in the process is to compute weight adjustment factors. These are
calculated by taking the ratio of the expected number of neutrons from the CMFD
source distribution to the current number of neutrons in each mesh. It is
straightforward to compute the CMFD number of neutrons because it is the
product between the total starting initial weight of neutrons and the CMFD
normalized fission source distribution. To compute the number of neutrons from
the current MC source, OpenMC sums the statistical
weights of neutrons from the source bank on a given spatial and energy mesh.
Once weight adjustment factors were calculated, each neutron's statistical
weight in the source bank was modified according to its location and energy.
Examples of CMFD simulations using OpenMC can be found in [HermanThesis]_.
.. only:: html
.. rubric:: References
.. [BEAVRS] Nick Horelik, Bryan Herman. *Benchmark for Evaluation And Verification of Reactor
Simulations*. Massachusetts Institute of Technology, https://crpg.mit.edu/research/beavrs
, 2013.
.. [Gill] Daniel F. Gill. *Newton-Krylov methods for the solution of the k-eigenvalue problem in
multigroup neutronics calculations*. Ph.D. thesis, Pennsylvania State University, 2010.
.. [Hebert] Alain Hebert. *Applied reactor physics*. Presses Internationales Polytechnique,
Montreal, 2009.
.. [Herman] Bryan R. Herman, Benoit Forget, Kord Smith, and Brian N. Aviles. Improved
diffusion coefficients generated from Monte Carlo codes. In *Proceedings of M&C
2013*, Sun Valley, ID, USA, May 5 - 9, 2013.
.. [HermanThesis] Bryan R. Herman. *Monte Carlo and Thermal Hydraulic Coupling using
Low-Order Nonlinear Diffusion Acceleration*. Sc.D. thesis,
Massachusetts Institute of Technology, 2014.
.. [Knoll] D.A. Knoll, H. Park, and C. Newman. *Acceleration of k-eigenvalue/criticality
calculations using the Jacobian-free Newton-Krylov method*. Nuclear Science and
Engineering, 167:133–140, 2011.
.. [Park] H. Park, D.A. Knoll, and C.K. Newman. *Nonlinear acceleration of transport
criticality problems*. Nuclear Science and Engineering, 172:52–65, 2012.
.. [Rhodes] Joel Rhodes and Malte Edenius. *CASMO-4 --- A Fuel Assembly Burnup Program.
User’s Manual*. Studsvik of America, ssp-09/443-u rev 0, proprietary edition, 2001.
.. [Smith] Kord S Smith and Joel D Rhodes III. *Full-core, 2-D, LWR core calculations with
CASMO-4E*. In Proceedings of PHYSOR 2002, Seoul, Korea, October 7 - 10, 2002.

280
docs/source/methods/cross_sections.rst Normal file
View file

@ -0,0 +1,280 @@
.. _methods_cross_sections:
=============================
Cross Section Representations
=============================
----------------------
Continuous-Energy Data
----------------------
In OpenMC, the data governing the interaction of neutrons with various nuclei
for continous-energy problems are represented using an HDF5 format that can be
produced by converting files in the ACE format, which is used by MCNP_ and
Serpent_. ACE-format data can be generated with the NJOY_ nuclear data
processing system, which converts raw `ENDF/B data`_ into linearly-interpolable
data as required by most Monte Carlo codes. Since ACE-format data can be
converted into OpenMC's HDF5 format, it is possible to perform direct comparison
of OpenMC with other codes using the same underlying nuclear data library.
The ACE format contains continuous-energy cross sections for the following types
of reactions: elastic scattering, fission (or first-chance fission,
second-chance fission, etc.), inelastic scattering, :math:`(n,xn)`,
:math:`(n,\gamma)`, and various other absorption reactions. For those reactions
with one or more neutrons in the exit channel, secondary angle and energy
distributions may be provided. In addition, fissionable nuclides have total,
prompt, and/or delayed :math:`\nu` as a function of energy and neutron precursor
distributions. Many nuclides also have probability tables to be used for
accurate treatment of self-shielding in the unresolved resonance range. For
bound scatterers, separate tables with :math:`S(\alpha,\beta,T)` scattering law
data can be used.
Energy Grid Methods
-------------------
The method by which continuous-energy cross sections for each nuclide in a
problem are stored as a function of energy can have a substantial effect on the
performance of a Monte Carlo simulation. Since the ACE format is based on
linearly-interpolable cross sections, each nuclide has cross sections tabulated
over a wide range of energies. Some nuclides may only have a few points
tabulated (e.g. H-1) whereas other nuclides may have hundreds or thousands of
points tabulated (e.g. U-238).
At each collision, it is necessary to sample the probability of having a
particular type of interaction whether it be elastic scattering, :math:`(n,2n)`,
level inelastic scattering, etc. This requires looking up the microscopic cross
sections for these reactions for each nuclide within the target material. Since
each nuclide has a unique energy grid, it would be necessary to search for the
appropriate index for each nuclide at every collision. This can become a very
time-consuming process, especially if there are many nuclides in a problem as
there would be for burnup calculations. Thus, there is a strong motive to
implement a method of reducing the number of energy grid searches in order to
speed up the calculation.
Logarithmic Mapping
+++++++++++++++++++
To speed up energy grid searches, OpenMC uses a `logarithmic mapping technique`_
to limit the range of energies that must be searched for each nuclide. The
entire energy range is divided up into equal-lethargy segments, and the bounding
energies of each segment are mapped to bounding indices on each of the nuclide
energy grids. By default, OpenMC uses 8000 equal-lethargy segments as
recommended by Brown.
Other Methods
+++++++++++++
A good survey of other energy grid techniques, including unionized energy grids,
can be found in a paper by Leppanen_.
.. _windowed_multipole:
Windowed Multipole Representation
---------------------------------
In addition to the usual pointwise representation of cross sections, OpenMC
offers support for a data format called windowed multipole (WMP). This data
format requires less memory than pointwise cross sections, and it allows
on-the-fly Doppler broadening to arbitrary temperature.
The multipole method was introduced by Hwang_ and the faster windowed multipole
method by Josey_. In the multipole format, cross section resonances are
represented by poles, :math:`p_j`, and residues, :math:`r_j`, in the complex
plane. The 0K cross sections in the resolved resonance region can be computed
by summing up a contribution from each pole:
.. math::
\sigma(E, T=0\text{K}) = \frac{1}{E} \sum_j \text{Re} \left[
\frac{i r_j}{\sqrt{E} - p_j} \right]
Assuming free-gas thermal motion, cross sections in the multipole form can be
analytically Doppler broadened to give the form:
.. math::
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[r_j
\sqrt{\pi} W_i(z) - \frac{r_j}{\sqrt{\pi}} C \left(\frac{p_j}{\sqrt{\xi}},
\frac{u}{2 \sqrt{\xi}}\right)\right]
.. math::
W_i(z) = \frac{i}{\pi} \int_{-\infty}^\infty dt \frac{e^{-t^2}}{z - t}
.. math::
C \left(\frac{p_j}{\sqrt{\xi}},\frac{u}{2 \sqrt{\xi}}\right) =
2p_j \int_0^\infty du' \frac{e^{-(u + u')^2/4\xi}}{p_j^2 - u'^2}
.. math::
z = \frac{\sqrt{E} - p_j}{2 \sqrt{\xi}}
.. math::
\xi = \frac{k_B T}{4 A}
.. math::
u = \sqrt{E}
where :math:`T` is the temperature of the resonant scatterer, :math:`k_B` is the
Boltzmann constant, :math:`A` is the mass of the target nucleus. For
:math:`E \gg k_b T/A`, the :math:`C` integral is approximately zero, simplifying
the cross section to:
.. math::
\sigma(E, T) = \frac{1}{2 E \sqrt{\xi}} \sum_j \text{Re} \left[i r_j
\sqrt{\pi} W_i(z)\right]
The :math:`W_i` integral simplifies down to an analytic form. We define the
Faddeeva function, :math:`W` as:
.. math::
W(z) = e^{-z^2} \text{Erfc}(-iz)
Through this, the integral transforms as follows:
.. math::
\text{Im} (z) > 0 : W_i(z) = W(z)
.. math::
\text{Im} (z) < 0 : W_i(z) = -W(z^*)^*
There are freely available algorithms_ to evaluate the Faddeeva function. For
many nuclides, the Faddeeva function needs to be evaluated thousands of times to
calculate a cross section. To mitigate that computational cost, the WMP method
only evaluates poles within a certain energy "window" around the incident
neutron energy and accounts for the effect of resonances outside that window
with a polynomial fit. This polynomial fit is then broadened exactly. This
exact broadening can make up for the removal of the :math:`C` integral, as
typically at low energies, only curve fits are used.
Note that the implementation of WMP in OpenMC currently assumes that inelastic
scattering does not occur in the resolved resonance region. This is usually,
but not always the case. Future library versions may eliminate this issue.
The data format used by OpenMC to represent windowed multipole data is specified
in :ref:`io_data_wmp` with a publicly available `WMP library`_.
.. _temperature_treatment:
Temperature Treatment
---------------------
At the beginning of a simulation, OpenMC collects a list of all temperatures
that are present in a model. It then uses this list to determine what cross
sections to load. The data that is loaded depends on what temperature method has
been selected. There are three methods available:
:Nearest: Cross sections are loaded only if they are within a specified
tolerance of the actual temperatures in the model.
:Interpolation: Cross sections are loaded at temperatures that bound the actual
temperatures in the model. During transport, cross sections for
each material are calculated using statistical linear-linear
interpolation between bounding temperature. Suppose cross
sections are available at temperatures :math:`T_1, T_2, ...,
T_n` and a material is assigned a temperature :math:`T` where
:math:`T_i < T < T_{i+1}`. Statistical interpolation is applied
as follows: a uniformly-distributed random number of the unit
interval, :math:`\xi`, is sampled. If :math:`\xi < (T -
T_i)/(T_{i+1} - T_i)`, then cross sections at temperature
:math:`T_{i+1}` are used. Otherwise, cross sections at
:math:`T_i` are used. This procedure is applied for pointwise
cross sections in the resolved resonance range, unresolved
resonance probability tables, and :math:`S(\alpha,\beta)`
thermal scattering tables.
:Multipole: Resolved resonance cross sections are calculated on-the-fly using
techniques/data described in :ref:`windowed_multipole`. Cross
section data is loaded for a single temperature and is used in the
unresolved resonance and fast energy ranges.
----------------
Multi-Group Data
----------------
The data governing the interaction of particles with various nuclei or materials
are represented using a multi-group library format specific to the OpenMC code.
The format is described in the :ref:`mgxs_lib_spec`. The data itself can be
prepared via traditional paths or directly from a continuous-energy OpenMC
calculation by use of the Python API as is shown in the
:ref:`notebook_mg_mode_part_i` example notebook. This multi-group library
consists of meta-data (such as the energy group structure) and multiple `xsdata`
objects which contains the required microscopic or macroscopic multi-group data.
At a minimum, the library must contain the absorption cross section
(:math:`\sigma_{a,g}`) and a scattering matrix. If the problem is an eigenvalue
problem then all fissionable materials must also contain either a fission
production matrix cross section (:math:`\nu\sigma_{f,g\rightarrow g'}`), or both
the fission spectrum data (:math:`\chi_{g'}`) and a fission production cross
section (:math:`\nu\sigma_{f,g}`), or, . The library must also contain the
fission cross section (:math:`\sigma_{f,g}`) or the fission energy release cross
section (:math:`\kappa\sigma_{f,g}`) if the associated tallies are required by
the model using the library.
After a scattering collision, the outgoing particle experiences a change in both
energy and angle. The probability of a particle resulting in a given outgoing
energy group (`g'`) given a certain incoming energy group (`g`) is provided by
the scattering matrix data. The angular information can be expressed either via
Legendre expansion of the particle's change-in-angle (:math:`\mu`), a tabular
representation of the probability distribution function of :math:`\mu`, or a
histogram representation of the same PDF. The formats used to represent these
are described in the :ref:`mgxs_lib_spec`.
Unlike the continuous-energy mode, the multi-group mode does not explicitly
track particles produced from scattering multiplication (i.e., :math:`(n,xn)`)
reactions. These are instead accounted for by adjusting the weight of the
particle after the collision such that the correct total weight is maintained.
The weight adjustment factor is optionally provided by the `multiplicity` data
which is required to be provided in the form of a group-wise matrix. This data
is provided as a group-wise matrix since the probability of producing multiple
particles in a scattering reaction depends on both the incoming energy, `g`, and
the sampled outgoing energy, `g'`. This data represents the average number of
particles emitted from a scattering reaction, given a scattering reaction has
occurred:
.. math::
multiplicity_{g \rightarrow g'} = \frac{\nu_{scatter}\sigma_{s,g \rightarrow g'}}{
\sigma_{s,g \rightarrow g'}}
If this scattering multiplication information is not provided in the library
then no weight adjustment will be performed. This is equivalent to neglecting
any additional particles produced in scattering multiplication reactions.
However, this assumption will result in a loss of accuracy since the total
particle population would not be conserved. This reduction in accuracy due to
the loss in particle conservation can be mitigated by reducing the absorption
cross section as needed to maintain particle conservation. This adjustment can
be done when generating the library, or by OpenMC. To have OpenMC perform the
adjustment, the total cross section (:math:`\sigma_{t,g}`) must be provided.
With this information, OpenMC will then adjust the absorption cross section as
follows:
.. math::
\sigma_{a,g} = \sigma_{t,g} - \sum_{g'}\nu_{scatter}\sigma_{s,g \rightarrow g'}
The above method is the same as is usually done with most deterministic solvers.
Note that this method is less accurate than using the scattering multiplication
weight adjustment since simply reducing the absorption cross section does not
include any information about the outgoing energy of the particles produced in
these reactions.
All of the data discussed in this section can be provided to the code
independent of the particle's direction of motion (i.e., isotropic), or the data
can be provided as a tabular distribution of the polar and azimuthal particle
direction angles. The isotropic representation is the most commonly used,
however inaccuracies are to be expected especially near material interfaces
where a material has a very large cross sections relative to the other material
(as can be expected in the resonance range). The angular representation can be
used to minimize this error.
Finally, the above options for representing the physics do not have to be
consistent across the problem. The number of groups and the structure, however,
does have to be consistent across the data sets. That is to say that each
microscopic or macroscopic data set does not have to apply the same scattering
expansion, treatment of multiplicity or angular representation of the cross
sections. This allows flexibility for the model to use highly anisotropic
scattering information in the water while the fuel can be simulated with linear
or even isotropic scattering.
.. _logarithmic mapping technique:
https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-14-24530.pdf
.. _Hwang: http://www.ans.org/pubs/journals/nse/a_16381
.. _Josey: https://doi.org/10.1016/j.jcp.2015.08.013
.. _WMP Library: https://github.com/mit-crpg/WMP_Library
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi
.. _NJOY: http://t2.lanl.gov/codes.shtml
.. _ENDF/B data: http://www.nndc.bnl.gov/endf
.. _Leppanen: https://doi.org/10.1016/j.anucene.2009.03.019
.. _algorithms: http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package

161
docs/source/methods/eigenvalue.rst Normal file
View file

@ -0,0 +1,161 @@
.. _methods_eigenvalue:
=======================
Eigenvalue Calculations
=======================
An eigenvalue calculation, also referred to as a criticality calculation, is a
transport simulation wherein the source of neutrons includes a fissionable
material. Some common eigenvalue calculations include the simulation of nuclear
reactors, spent fuel pools, nuclear weapons, and other fissile systems. The
reason they are called *eigenvalue* calculations is that the transport equation
becomes an eigenvalue equation if a fissionable source is present since then the
source of neutrons will depend on the flux of neutrons itself. Eigenvalue
simulations using Monte Carlo methods are becoming increasingly common with the
advent of high-performance computing.
This section will explore the theory behind and implementation of eigenvalue
calculations in a Monte Carlo code.
.. _method-successive-generations:
--------------------------------
Method of Successive Generations
--------------------------------
The method used to converge on the fission source distribution in an eigenvalue
calculation, known as the method of successive generations, was first introduced
by [Lieberoth]_. In this method, a finite number of neutron histories,
:math:`N`, are tracked through their lifetime iteratively. If fission occurs,
rather than tracking the resulting fission neutrons, the spatial coordinates of
the fission site, the sampled outgoing energy and direction of the fission
neutron, and the weight of the neutron are stored for use in the subsequent
generation. In OpenMC, the array used for storing the fission site information
is called the *fission bank*. At the end of each fission generation, :math:`N`
source sites for the next generation must be randomly sampled from the :math:`M`
fission sites that were stored to ensure that the neutron population does not
grow exponentially. The sampled source sites are stored in an array called the
*source bank* and can be retrieved during the subsequent generation.
It's important to recognize that in the method of successive generations, we
must start with some assumption on how the fission source sites are distributed
since the distribution is not known *a priori*. Typically, a user will make a
guess as to what the distribution is -- this guess could be a uniform
distribution over some region of the geometry or simply a point
source. Fortunately, regardless of the choice of initial source distribution,
the method is guaranteed to converge to the true source distribution. Until the
source distribution converges, tallies should not be scored to since they will
otherwise include contributions from an unconverged source distribution.
The method by which the fission source iterations are parallelized can have a
large impact on the achievable parallel scaling. This topic is discussed at length
in :ref:`fission-bank-algorithms`.
-------------------------
Source Convergence Issues
-------------------------
Diagnosing Convergence with Shannon Entropy
-------------------------------------------
As discussed earlier, it is necessary to converge both :math:`k_{eff}` and the
source distribution before any tallies can begin. Moreover, the convergence rate
of the source distribution is in general slower than that of
:math:`k_{eff}`. One should thus examine not only the convergence of
:math:`k_{eff}` but also the convergence of the source distribution in order to
make decisions on when to start active batches.
However, the representation of the source distribution makes it a bit more
difficult to analyze its convergence. Since :math:`k_{eff}` is a scalar
quantity, it is easy to simply look at a line plot of :math:`k_{eff}` versus the
number of batches and this should give the user some idea about whether it has
converged. On the other hand, the source distribution at any given batch is a
finite set of coordinates in Euclidean space. In order to analyze the
convergence, we would either need to use a method for assessing convergence of
an N-dimensional quantity or transform our set of coordinates into a scalar
metric. The latter approach has been developed considerably over the last decade
and a method now commonly used in Monte Carlo eigenvalue calculations is to use
a metric called the `Shannon entropy`_, a concept borrowed from information
theory.
To compute the Shannon entropy of the source distribution, we first need to
discretize the source distribution rather than having a set of coordinates in
Euclidean space. This can be done by superimposing a structured mesh over the
geometry (containing at least all fissionable materials). Then, the fraction of
source sites that are present in each mesh element is counted:
.. math::
:label: fraction-source
S_i = \frac{\text{Source sites in $i$-th mesh element}}{\text{Total number of
source sites}}
The Shannon entropy is then computed as
.. math::
:label: shannon-entropy
H = - \sum_{i=1}^N S_i \log_2 S_i
where :math:`N` is the number of mesh elements. With equation
:eq:`shannon-entropy`, we now have a scalar metric that we can use to assess the
convergence of the source distribution by observing line plots of the Shannon
entropy versus the number of batches.
In recent years, researchers have started looking at ways of automatically
assessing source convergence to relieve the burden on the user of having to look
at plots of :math:`k_{eff}` and the Shannon entropy. A number of methods have
been proposed (see e.g. [Romano]_, [Ueki]_), but each of these is not without
problems.
---------------------------
Uniform Fission Site Method
---------------------------
Generally speaking, the variance of a Monte Carlo tally will be inversely
proportional to the number of events that score to the tally. In a reactor
problem, this implies that regions with low relative power density will have
higher variance that regions with high relative power density. One method to
circumvent the uneven distribution of relative errors is the uniform fission
site (UFS) method introduced by [Sutton]_. In this method, the portion of the
problem containing fissionable material is subdivided into a number of cells
(typically using a structured mesh). Rather than producing
.. math::
m = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t}
fission sites at each collision where :math:`w` is the weight of the neutron,
:math:`k` is the previous-generation estimate of the neutron multiplication
factor, :math:`\nu\Sigma_f` is the neutron production cross section, and
:math:`\Sigma_t` is the total cross section, in the UFS method we produce
.. math::
m_{UFS} = \frac{w}{k} \frac{\nu\Sigma_f}{\Sigma_t} \frac{v_i}{s_i}
fission sites at each collision where :math:`v_i` is the fraction of the total
volume occupied by cell :math:`i` and :math:`s_i` is the fraction of the fission
source contained in cell :math:`i`. To ensure that no bias is introduced, the
weight of each fission site stored in the fission bank is :math:`s_i/v_i` rather
than unity. By ensuring that the expected number of fission sites in each mesh
cell is constant, the collision density across all cells, and hence the variance
of tallies, is more uniform than it would be otherwise.
.. _Shannon entropy: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-06-3737.pdf
.. [Lieberoth] J. Lieberoth, "A Monte Carlo Technique to Solve the Static
Eigenvalue Problem of the Boltzmann Transport Equation," *Nukleonik*, **11**,
213-219 (1968).
.. [Romano] Paul K. Romano, "Application of the Stochastic Oscillator to Assess
Source Convergence in Monte Carlo Criticality Calculations,"
*Proc. International Conference on Mathematics, Computational Methods, and
Reactor Physics*, Saratoga Springs, New York (2009).
.. [Sutton] Daniel J. Kelly, Thomas M. Sutton, and Stephen C. Wilson, "MC21
Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark
Problem," *Proc. PHYSOR 2012*, Knoxville, Tennessee, Apr. 15--20 (2012).
.. [Ueki] Taro Ueki, "On-the-Fly Judgments of Monte Carlo Fission Source
Convergence," *Trans. Am. Nucl. Soc.*, **98**, 512 (2008).

149
docs/source/methods/energy_deposition.rst Normal file
View file

@ -0,0 +1,149 @@
.. _methods_heating:
=============================
Heating and Energy Deposition
=============================
As particles traverse a problem, some portion of their energy is deposited at
collision sites. This energy is deposited when charged particles, including
electrons and recoil nuclei, undergo electromagnetic interactions with
surrounding electons and ions. The information describing how much energy
is deposited for a specific reaction is referred to as
"heating numbers" and can be computed using a program like NJOY with the
``heatr`` module.
These heating rate is the product of reaction-specific coefficients and
a reaction cross section
.. math::
H(E) = \phi(E)\sum_i\rho_i\sum_rk_{i, r}(E),
and has units energy per time, typically eV / s.
Here, :math:`k_{i, r}` are the KERMA (Kinetic Energy Release in Materials)
[Mack97]_ coefficients for reaction :math:`r` of isotope :math:`i`.
The KERMA coefficients have units energy :math:`\times` cross-section, e.g.
eV-barn, and can be used much like a reaction cross section for the purpose
of tallying energy deposition.
KERMA coefficients can be computed using the energy-balance method with
a nuclear data processing code like NJOY, which performs the following
iteration over all reactions :math:`r` for all isotopes :math:`i`
requested
.. math::
k_{i, r}(E) = \left(E + Q_{i, r} - \bar{E}_{i, r, n}
- \bar{E}_{i, r, \gamma}\right)\sigma_{i, r}(E),
removing the energy of neutral particles (neutrons and photons) that are
transported away from the reaction site :math:`\bar{E}`, and the reaction
:math:`Q` value.
-------
Fission
-------
During a fission event, there are potentially many secondary particles, and all
must be considered. The total energy released in a fission event is typically
broken up into the following categories:
- :math:`E_{fr}` - kinetic energy of fission fragments
- :math:`E_{n,p}` - energy of prompt fission neutrons
- :math:`E_{n,d}` - energy of delayed fission neutrons
- :math:`E_{\gamma,p}` - energy of prompt fission photons
- :math:`E_{\gamma,d}` - energy of delayed fission photons
- :math:`E_{\beta}` - energy of released :math:`\beta` particles
- :math:`E_{\nu}` - energy of neutrinos
These components are defined in MF=1,MT=458 data in a standard ENDF/B-6 formatted
file. All these quantities may depend upon incident neutron energy,
but this dependence is not shown to make the following demonstrations cleaner.
As neutrinos scarcely interact with matter, the recoverable energy from
fission is defined as
.. math::
E_r\equiv E_{fr} + E_{n,p} + E_{n, d} + E_{\gamma, p}
+ E_{\gamma, d} + E_{\beta}
Furthermore, the energy of the secondary neutrons and photons is given as
:math:`E_{n, p}` and :math:`E_{\gamma, p}`, respectively.
NJOY computes the fission KERMA coefficient using this energy-balance method to be
.. math::
k_{i, f}(E) = \left[E + Q(E) - \bar{E}(E)\right]\sigma_{i, f}(E)
= \left[E_{fr} + E_{\gamma, p}\right]\sigma_{i, j}(E)
.. note::
The energy from delayed neutrons and photons and beta particles is intentionally
left out from the NJOY calculations.
---------------------
OpenMC Implementation
---------------------
For fissile isotopes, OpenMC makes modifications to the heating reaction to
include all relevant components of fission energy release. These modifications
are made to the total heating reaction, MT=301. Breaking the total heating
KERMA into a fission and non-fission section, one can write
.. math::
k_i(E) = k_{i, nf}(E) + \left[E_{fr}(E) + E_{\gamma, p}\right]\sigma_{i, f}(E)
OpenMC seeks to modify the total heating data to include energy from
:math:`\beta` particles and, conditionally, delayed photons. This conditional
inclusion depends on the simulation mode: neutron transport, or coupled
neutron-photon transport. The heating due to fission is removed using MT=318
data, and then re-built using the desired components of fission energy release
from MF=1,MT=458 data.
Neutron Transport
-----------------
For this case, OpenMC instructs ``heatr`` to produce heating coefficients
assuming that energy from photons, :math:`E_{\gamma, p}` and
:math:`E_{\gamma, d}`, is deposited at the fission site.
Let :math:`N901` represent the total heating number returned from this ``heatr``
run with :math:`N918` reflecting fission heating computed from NJOY.
:math:`M901` represent the following modification
.. math::
M901_{i}(E)\equiv N901_{i}(E) - N918_{i}(E)
+ \left[E_{i, fr} + E_{i, \beta} + E_{i, \gamma, p}
+ E_{i, \gamma, d}\right]\sigma_{i, f}(E).
This modified heating data is stored as the MT=901 reaction and will be scored
if ``heating-local`` is included in :attr:`openmc.Tally.scores`.
Coupled neutron-photon transport
--------------------------------
Here, OpenMC instructs ``heatr`` to assume that energy from photons is not
deposited locally. However, the definitions provided in the NJOY manual
indicate that, regardless of this mode, the prompt photon energy is still
included in :math:`k_{i, f}`, and therefore must be manually removed.
Let :math:`N301` represent the total heating number returned from this
``heatr`` run and :math:`M301` be
.. math::
M301_{i}(E)\equiv N301_{i}(E) - N318_{i}(E)
+ \left[E_{i, fr}(E) + E_{i, \beta}(E)\right]\sigma_{i, f}(E).
This modified heating data is stored as the MT=301 reaction and will be scored
if ``heating`` is included in :attr:`openmc.Tally.scores`.
----------
References
----------
.. [Mack97] Abdou, M.A., Maynard, C.W., and Wright, R.Q. MACK: computer
program to calculate neutron energy release parameters (fluence-to-kerma
factors) and multigroup neutron reaction cross sections from nuclear data
in ENDF Format. Oak Ridge National Laboratory report ORNL-TM-3994.

967
docs/source/methods/geometry.rst Normal file
View file

@ -0,0 +1,967 @@
.. _methods_geometry:
========
Geometry
========
---------------------------
Constructive Solid Geometry
---------------------------
OpenMC uses a technique known as `constructive solid geometry`_ (CSG) to build
arbitrarily complex three-dimensional models in Euclidean space. In a CSG model,
every unique object is described as the union and/or intersection of
*half-spaces* created by bounding `surfaces`_. Every surface divides all of
space into exactly two half-spaces. We can mathematically define a surface as a
collection of points that satisfy an equation of the form :math:`f(x,y,z) = 0`
where :math:`f(x,y,z)` is a given function. All coordinates for which
:math:`f(x,y,z) < 0` are referred to as the negative half-space (or simply the
*negative side*) and coordinates for which :math:`f(x,y,z) > 0` are referred to
as the positive half-space.
Let us take the example of a sphere centered at the point :math:`(x_0,y_0,z_0)`
with radius :math:`R`. One would normally write the equation of the sphere as
.. math::
:label: sphere-equation
(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2
By subtracting the right-hand term from both sides of equation
:eq:`sphere-equation`, we can then write the surface equation for the sphere:
.. math::
:label: surface-equation-sphere
f(x,y,z) = (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0
One can confirm that any point inside this sphere will correspond to
:math:`f(x,y,z) < 0` and any point outside the sphere will correspond to
:math:`f(x,y,z) > 0`.
In OpenMC, every surface defined by the user is assigned an integer to uniquely
identify it. We can then refer to either of the two half-spaces created by a
surface by a combination of the unique ID of the surface and a positive/negative
sign. Figure :num:`fig-halfspace` shows an example of an ellipse with unique ID 1
dividing space into two half-spaces.
.. _fig-halfspace:
.. figure:: ../_images/halfspace.svg
:align: center
:figclass: align-center
Example of an ellipse and its associated half-spaces.
References to half-spaces created by surfaces are used to define regions of
space of uniform composition, which are then assigned to *cells*. OpenMC allows
regions to be defined using union, intersection, and complement operators. As in
MCNP_, the intersection operator is implicit as doesn't need to be written in a
region specification. A defined region is then associated with a material
composition in a cell. Figure :num:`fig-union` shows an example of a cell region
defined as the intersection of an ellipse and two planes.
.. _fig-union:
.. figure:: ../_images/union.svg
:align: center
:figclass: align-center
The shaded region represents a cell bounded by three surfaces.
The ability to form regions based on bounding quadratic surfaces enables OpenMC
to model arbitrarily complex three-dimensional objects. In practice, one is
limited only by the different surface types available in OpenMC. The following
table lists the available surface types, the identifier used to specify them in
input files, the corresponding surface equation, and the input parameters needed
to fully define the surface.
.. table:: Surface types available in OpenMC.
+----------------------+------------+------------------------------+-------------------------+
| Surface | Identifier | Equation | Parameters |
+======================+============+==============================+=========================+
| Plane perpendicular | x-plane | :math:`x - x_0 = 0` | :math:`x_0` |
| to :math:`x`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Plane perpendicular | y-plane | :math:`y - y_0 = 0` | :math:`y_0` |
| to :math:`y`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Plane perpendicular | z-plane | :math:`z - z_0 = 0` | :math:`z_0` |
| to :math:`z`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Arbitrary plane | plane | :math:`Ax + By + Cz = D` | :math:`A\;B\;C\;D` |
+----------------------+------------+------------------------------+-------------------------+
| Infinite cylinder | x-cylinder | :math:`(y-y_0)^2 + (z-z_0)^2 | :math:`y_0\;z_0\;R` |
| parallel to | | = R^2` | |
| :math:`x`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Infinite cylinder | y-cylinder | :math:`(x-x_0)^2 + (z-z_0)^2 | :math:`x_0\;z_0\;R` |
| parallel to | | = R^2` | |
| :math:`y`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Infinite cylinder | z-cylinder | :math:`(x-x_0)^2 + (y-y_0)^2 | :math:`x_0\;y_0\;R` |
| parallel to | | = R^2` | |
| :math:`z`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Sphere | sphere | :math:`(x-x_0)^2 + (y-y_0)^2 | :math:`x_0 \; y_0 \; |
| | | + (z-z_0)^2 = R^2` | z_0 \; R` |
+----------------------+------------+------------------------------+-------------------------+
| Cone parallel to the | x-cone | :math:`(y-y_0)^2 + (z-z_0)^2 | :math:`x_0 \; y_0 \; |
| :math:`x`-axis | | = R^2(x-x_0)^2` | z_0 \; R^2` |
+----------------------+------------+------------------------------+-------------------------+
| Cone parallel to the | y-cone | :math:`(x-x_0)^2 + (z-z_0)^2 | :math:`x_0 \; y_0 \; |
| :math:`y`-axis | | = R^2(y-y_0)^2` | z_0 \; R^2` |
+----------------------+------------+------------------------------+-------------------------+
| Cone parallel to the | z-cone | :math:`(x-x_0)^2 + (y-y_0)^2 | :math:`x_0 \; y_0 \; |
| :math:`z`-axis | | = R^2(z-z_0)^2` | z_0 \; R^2` |
+----------------------+------------+------------------------------+-------------------------+
| General quadric | quadric | :math:`Ax^2 + By^2 + Cz^2 + | :math:`A \; B \; C \; D |
| surface | | Dxy + Eyz + Fxz + Gx + Hy + | \; E \; F \; G \; H \; |
| | | Jz + K` | J \; K` |
+----------------------+------------+------------------------------+-------------------------+
.. _universes:
Universes
---------
OpenMC supports universe-based geometry similar to the likes of MCNP_ and
Serpent_. This capability enables user to model any identical repeated
structures once and then fill them in various spots in the geometry. A
prototypical example of a repeated structure would be a fuel pin within a fuel
assembly or a fuel assembly within a core.
Each cell in OpenMC can either be filled with a normal material or with a
universe. If the cell is filled with a universe, only the region of the universe
that is within the defined boundaries of the parent cell will be present in the
geometry. That is to say, even though a collection of cells in a universe may
extend to infinity, not all of the universe will be "visible" in the geometry
since it will be truncated by the boundaries of the cell that contains it.
When a cell is filled with a universe, it is possible to specify that the
universe filling the cell should be rotated and translated. This is done through
the ``rotation`` and ``translation`` attributes on a cell (note though that
these can only be specified on a cell that is filled with another universe, not
a material).
It is not necessary to use or assign universes in a geometry if there are no
repeated structures. Any cell in the geometry that is not assigned to a
specified universe is automatically part of the *base universe* whose
coordinates are just the normal coordinates in Euclidean space.
Lattices
--------
Often times, repeated structures in a geometry occur in a regular pattern such
as a rectangular or hexagonal lattice. In such a case, it would be cumbersome
for a user to have to define the boundaries of each of the cells to be filled
with a universe. Thus, OpenMC provides a lattice capability similar to that used
in MCNP_ and Serpent_.
The implementation of lattices is similar in principle to universes --- instead
of a cell being filled with a universe, the user can specify that it is filled
with a finite lattice. The lattice is then defined by a two-dimensional array of
universes that are to fill each position in the lattice. A good example of the
use of lattices and universes can be seen in the OpenMC model for the `Monte
Carlo Performance benchmark`_.
------------------------------------------
Computing the Distance to Nearest Boundary
------------------------------------------
One of the most basic algorithms in any Monte Carlo code is determining the
distance to the nearest surface within a cell. Since each cell is defined by
the surfaces that bound it, if we compute the distance to all surfaces bounding
a cell, we can determine the nearest one.
With the possibility of a particle having coordinates on multiple levels
(universes) in a geometry, we must exercise care when calculating the distance
to the nearest surface. Each different level of geometry has a set of boundaries
with which the particle's direction of travel may intersect. Thus, it is
necessary to check the distance to the surfaces bounding the cell in each
level. This should be done starting the highest (most global) level going down
to the lowest (most local) level. That ensures that if two surfaces on different
levels are coincident, by default the one on the higher level will be selected
as the nearest surface. Although they are not explicitly defined, it is also
necessary to check the distance to surfaces representing lattice boundaries if a
lattice exists on a given level.
The following procedure is used to calculate the distance to each bounding
surface. Suppose we have a particle at :math:`(x_0,y_0,z_0)` traveling in the
direction :math:`u_0,v_0,w_0`. To find the distance :math:`d` to a surface
:math:`f(x,y,z) = 0`, we need to solve the equation:
.. math::
:label: dist-to-boundary-1
f(x_0 + du_0, y_0 + dv_0, z_0 + dw_0) = 0
If no solutions to equation :eq:`dist-to-boundary-1` exist or the only solutions
are complex, then the particle's direction of travel will not intersect the
surface. If the solution to equation :eq:`dist-to-boundary-1` is negative, this
means that the surface is "behind" the particle, i.e. if the particle continues
traveling in its current direction, it will not hit the surface. The complete
derivation for different types of surfaces used in OpenMC will be presented in
the following sections.
Since :math:`f(x,y,z)` in general is quadratic in :math:`x`, :math:`y`, and
:math:`z`, this implies that :math:`f(x_0 + du_0, y + dv_0, z + dw_0)` is
quadratic in :math:`d`. Thus we expect at most two real solutions to
:eq:`dist-to-boundary-1`. If no solutions to :eq:`dist-to-boundary-1` exist or
the only solutions are complex, then the particle's direction of travel will not
intersect the surface. If the solution to :eq:`dist-to-boundary-1` is negative,
this means that the surface is "behind" the particle, i.e. if the particle
continues traveling in its current direction, it will not hit the surface.
Once a distance has been computed to a surface, we need to check if it is closer
than previously-computed distances to surfaces. Unfortunately, we cannot just
use the minimum function because some of the calculated distances, which should
be the same in theory (e.g. coincident surfaces), may be slightly different due
to the use of floating-point arithmetic. Consequently, we should first check for
floating-point equality of the current distance calculated and the minimum found
thus far. This is done by checking if
.. math::
:label: fp-distance
\frac{| d - d_{min} |}{d_{min}} < \epsilon
where :math:`d` is the distance to a surface just calculated, :math:`d_{min}` is
the minimum distance found thus far, and :math:`\epsilon` is a small number. In
OpenMC, this parameter is set to :math:`\epsilon = 10^{-14}` since all floating
calculations are done on 8-byte floating point numbers.
Plane Perpendicular to an Axis
------------------------------
The equation for a plane perpendicular to, for example, the x-axis is simply
:math:`x - x_0 = 0`. As such, we need to solve :math:`x + du - x_0 = 0`. The
solution for the distance is
.. math::
:label: dist-xplane
d = \frac{x_0 - x}{u}
Note that if the particle's direction of flight is parallel to the x-axis,
i.e. :math:`u = 0`, the distance to the surface will be infinity. While the
example here was for a plane perpendicular to the x-axis, the same formula can
be applied for the surfaces :math:`y = y_0` and :math:`z = z_0`.
Generic Plane
-------------
The equation for a generic plane is :math:`Ax + By + Cz = D`. Thus, we need to
solve the equation :math:`A(x + du) + B(y + dv) + C(z + dw) = D`. The solution
to this equation for the distance is
.. math::
:label: dist-plane
d = \frac{D - Ax - By - Cz}{Au + Bv + Cw}
Again, we need to check whether the denominator is zero. If so, this means that
the particle's direction of flight is parallel to the plane and it will
therefore never hit the plane.
.. _cylinder_distance:
Cylinder Parallel to an Axis
----------------------------
The equation for a cylinder parallel to, for example, the x-axis is :math:`(y -
y_0)^2 + (z - z_0)^2 = R^2`. Thus, we need to solve :math:`(y + dv - y_0)^2 +
(z + dw - z_0)^2 = R^2`. Let us define :math:`\bar{y} = y - y_0` and
:math:`\bar{z} = z - z_0`. We then have
.. math::
:label: dist-xcylinder-1
(\bar{y} + dv)^2 + (\bar{z} + dw)^2 = R^2
Expanding equation :eq:`dist-xcylinder-1` and rearranging terms, we obtain
.. math::
:label: dist-xcylinder-2
(v^2 + w^2) d^2 + 2 (\bar{y}v + \bar{z}w) d + (\bar{y}^2 + \bar{z}^2 - R^2)
= 0
This is a quadratic equation for :math:`d`. To simplify notation, let us define
:math:`a = v^2 + w^2`, :math:`k = \bar{y}v + \bar{z}w`, and :math:`c =
\bar{y}^2 + \bar{z}^2 - R^2`. Thus, the distance is just the solution to
:math:`ad^2 + 2kd + c = 0`:
.. math::
:label: dist-xcylinder-3
d = \frac{-k \pm \sqrt{k^2 - ac}}{a}
A few conditions must be checked for. If :math:`a = 0`, this means the particle
is parallel to the cylinder and will thus never intersect it. Also, if
:math:`k^2 - ac < 0`, this means that both solutions to the quadratic are
complex. In physical terms, this means that the ray along which the particle is
traveling does not make any intersections with the cylinder.
If we do have intersections and :math:`c < 0`, this means that the particle is
inside the cylinder. Thus, one solution should be positive and one should be
negative. Clearly, the positive distance will occur when the sign on the
square root of the discriminant is positive since :math:`a > 0`.
If we have intersections and :math:`c > 0` this means that the particle is
outside the cylinder. Thus, the solutions to the quadratic are either both
positive or both negative. If they are both positive, the smaller (closer) one
will be the solution with a negative sign on the square root of the
discriminant.
The same equations and logic here can be used for cylinders that are parallel to
the y- or z-axis with appropriate substitution of constants.
Sphere
------
The equation for a sphere is :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 =
R^2`. Thus, we need to solve the equation
.. math::
:label: dist-sphere-1
(x + du - x_0)^2 + (y + dv - y_0)^2 + (z + dw - z_0)^2 = R^2
Let us define :math:`\bar{x} = x - x_0`, :math:`\bar{y} = y - y_0`, and
:math:`\bar{z} = z - z_0`. We then have
.. math::
:label: dist-sphere-2
(\bar{x} + du)^2 + (\bar{y} + dv)^2 + (\bar{z} - dw)^2 = R^2
Expanding equation :eq:`dist-sphere-2` and rearranging terms, we obtain
.. math::
:label: dist-sphere-3
d^2 + 2 (\bar{x}u + \bar{y}v + \bar{z}w) d + (\bar{x}^2 + \bar{y}^2 +
\bar{z}^2 - R^2) = 0
This is a quadratic equation for :math:`d`. To simplify notation, let us define
:math:`k = \bar{x}u + \bar{y}v + \bar{z}w` and :math:`c = \bar{x}^2 +
\bar{y}^2 + \bar{z}^2 - R^2`. Thus, the distance is just the solution to
:math:`d^2 + 2kd + c = 0`:
.. math::
:label: dist-sphere-4
d = -k \pm \sqrt{k^2 - c}
If the discriminant :math:`k^2 - c < 0`, this means that both solutions to the
quadratic are complex. In physical terms, this means that the ray along which
the particle is traveling does not make any intersections with the sphere.
If we do have intersections and :math:`c < 0`, this means that the particle is
inside the sphere. Thus, one solution should be positive and one should be
negative. The positive distance will occur when the sign on the square root of
the discriminant is positive. If we have intersections but :math:`c > 0` this
means that the particle is outside the sphere. The solutions to the quadratic
will then be either both positive or both negative. If they are both positive,
the smaller (closer) one will be the solution with a negative sign on the square
root of the discriminant.
Cone Parallel to an Axis
------------------------
The equation for a cone parallel to, for example, the x-axis is :math:`(y -
y_0)^2 + (z - z_0)^2 = R^2(x - x_0)^2`. Thus, we need to solve :math:`(y + dv -
y_0)^2 + (z + dw - z_0)^2 = R^2(x + du - x_0)^2`. Let us define :math:`\bar{x} =
x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. We then have
.. math::
:label: dist-xcone-1
(\bar{y} + dv)^2 + (\bar{z} + dw)^2 = R^2(\bar{x} + du)^2
Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain
.. math::
:label: dist-xcone-2
(v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d +
(\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0
Defining the terms
.. math::
:label: dist-xcone-terms
a = v^2 + w^2 - R^2u^2
k = \bar{y}v + \bar{z}w - R^2\bar{x}u
c = \bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2
we then have the simple quadratic equation :math:`ad^2 + 2kd + c = 0` which can
be solved as described in :ref:`cylinder_distance`.
General Quadric
---------------
The equation for a general quadric surface is :math:`Ax^2 + By^2 + Cz^2 + Dxy +
Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, we need to solve the equation
.. math::
:label: dist-quadric-1
A(x+du)^2 + B(y+dv)^2 + C(z+dw)^2 + D(x+du)(y+dv) + E(y+dv)(z+dw) + \\
F(x+du)(z+dw) + G(x+du) + H(y+dv) + J(z+dw) + K = 0
Expanding equation :eq:`dist-quadric-1` and rearranging terms, we obtain
.. math::
:label: dist-quadric-2
d^2(uv + vw + uw) + 2d(Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + \\
F(wx + uz))/2) + (x(Ax + Dy) + y(By + Ez) + z(Cz + Fx)) = 0
Defining the terms
.. math::
:label: dist-quadric-terms
a = uv + vw + uw
k = Aux + Bvy + Cwx + (D(uv + vx) + E(vz + wy) + F(wx + uz))/2
c = x(Ax + Dy) + y(By + Ez) + z(Cz + Fx)
we then have the simple quadratic equation :math:`ad^2 + 2kd + c = 0` which can
be solved as described in :ref:`cylinder_distance`.
.. _find-cell:
----------------------------
Finding a Cell Given a Point
----------------------------
Another basic algorithm is to determine which cell contains a given point in the
global coordinate system, i.e. if the particle's position is :math:`(x,y,z)`,
what cell is it currently in. This is done in the following manner in
OpenMC. With the possibility of multiple levels of coordinates, we must perform
a recursive search for the cell. First, we start in the highest (most global)
universe, which we call the base universe, and loop over each cell within
that universe. For each cell, we check whether the specified point is inside the
cell using the algorithm described in :ref:`cell-contains`. If the cell is
filled with a normal material, the search is done and we have identified the
cell containing the point. If the cell is filled with another universe, we then
search all cells within that universe to see if any of them contain the
specified point. If the cell is filled with a lattice, the position within the
lattice is determined, and then whatever universe fills that lattice position is
recursively searched. The search ends once a cell containing a normal material
is found that contains the specified point.
.. _cell-contains:
----------------------
Finding a Lattice Tile
----------------------
If a particle is inside a lattice, its position inside the lattice must be
determined before assigning it to a cell. Throughout this section, the
volumetric units of the lattice will be referred to as "tiles". Tiles are
identified by thier indices, and the process of discovering which tile contains
the particle is referred to as "indexing".
Rectilinear Lattice Indexing
----------------------------
Indices are assigned to tiles in a rectilinear lattice based on the tile's
position along the :math:`x`, :math:`y`, and :math:`z` axes. Figure
:num:`fig-rect-lat` maps the indices for a 2D lattice. The indices, (1, 1),
map to the lower-left tile. (5, 1) and (5, 5) map to the lower-right and
upper-right tiles, respectively.
.. _fig-rect-lat:
.. figure:: ../_images/rect_lat.svg
:align: center
:figclass: align-center
:width: 400px
Rectilinear lattice tile indices.
In general, a lattice tile is specified by the three indices,
:math:`(i_x, i_y, i_z)`. If a particle's current coordinates are
:math:`(x, y, z)` then the indices can be determined from these formulas:
.. math::
:label: rect_indexing
i_x = \left \lceil \frac{x - x_0}{p_0} \right \rceil
i_y = \left \lceil \frac{y - y_0}{p_1} \right \rceil
i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil
where :math:`(x_0, y_0, z_0)` are the coordinates to the lower-left-bottom
corner of the lattice, and :math:`p_0, p_1, p_2` are the pitches along the
:math:`x`, :math:`y`, and :math:`z` axes, respectively.
.. _hexagonal_indexing:
Hexagonal Lattice Indexing
--------------------------
A skewed coordinate system is used for indexing hexagonal lattice tiles.
Rather than a :math:`y`-axis, another axis is used that is rotated 30 degrees
counter-clockwise from the :math:`y`-axis. This axis is referred to as the
:math:`\alpha`-axis. Figure :num:`fig-hex-lat` shows how 2D hexagonal tiles
are mapped with the :math:`(x, \alpha)` basis. In this system, (0, 0) maps to
the center tile, (0, 2) to the top tile, and (2, -1) to the middle tile on the
right side.
.. _fig-hex-lat:
.. figure:: ../_images/hex_lat.svg
:align: center
:figclass: align-center
:width: 400px
Hexagonal lattice tile indices.
Unfortunately, the indices cannot be determined with one simple formula as
before. Indexing requires a two-step process, a coarse step which determines a
set of four tiles that contains the particle and a fine step that determines
which of those four tiles actually contains the particle.
In the first step, indices are found using these formulas:
.. math::
:label: hex_indexing
\alpha = -\frac{x}{\sqrt{3}} + y
i_x^* = \left \lfloor \frac{x}{p_0 \sqrt{3} / 2} \right \rfloor
i_\alpha^* = \left \lfloor \frac{\alpha}{p_0} \right \rfloor
where :math:`p_0` is the lattice pitch (in the :math:`x`-:math:`y` plane). The
true index of the particle could be :math:`(i_x^*, i_\alpha^*)`,
:math:`(i_x^* + 1, i_\alpha^*)`, :math:`(i_x^*, i_\alpha^* + 1)`, or
:math:`(i_x^* + 1, i_\alpha^* + 1)`.
The second step selects the correct tile from that neighborhood of 4. OpenMC
does this by calculating the distance between the particle and the centers of
each of the 4 tiles, and then picking the closest tile. This works because
regular hexagonal tiles form a Voronoi tessellation which means that all of the
points within a tile are closest to the center of that same tile.
Indexing along the :math:`z`-axis uses the same method from rectilinear
lattices, i.e.
.. math::
:label: hex_indexing_z
i_z = \left \lceil \frac{z - z_0}{p_2} \right \rceil
----------------------------------------
Determining if a Coordinate is in a Cell
----------------------------------------
To determine which cell a particle is in given its coordinates, we need to be
able to check whether a given cell contains a point. The algorithm for
determining if a cell contains a point is as follows. For each surface that
bounds a cell, we determine the particle's sense with respect to the surface. As
explained earlier, if we have a point :math:`(x_0,y_0,z_0)` and a surface
:math:`f(x,y,z) = 0`, the point is said to have negative sense if
:math:`f(x_0,y_0,z_0) < 0` and positive sense if :math:`f(x_0,y_0,z_0) > 0`. If
for all surfaces, the sense of the particle with respect to the surface matches
the specified sense that defines the half-space within the cell, then the point
is inside the cell. Note that this algorithm works only for *simple cells*
defined as intersections of half-spaces.
It may help to illustrate this algorithm using a simple example. Let's say we
have a cell defined as
.. code-block:: xml
<surface id="1" type="sphere" coeffs="0 0 0 10" />
<surface id="2" type="x-plane" coeffs="-3" />
<surface id="3" type="y-plane" coeffs="2" />
<cell id="1" surfaces="-1 2 -3" />
This means that the cell is defined as the intersection of the negative half
space of a sphere, the positive half-space of an x-plane, and the negative
half-space of a y-plane. Said another way, any point inside this cell must
satisfy the following equations
.. math::
:label: cell-contains-example
x^2 + y^2 + z^2 - 10^2 < 0 \\
x - (-3) > 0 \\
y - 2 < 0
In order to determine if a point is inside the cell, we would substitute its
coordinates into equation :eq:`cell-contains-example`. If the inequalities are
satisfied, than the point is indeed inside the cell.
--------------------------
Handling Surface Crossings
--------------------------
A particle will cross a surface if the distance to the nearest surface is closer
than the distance sampled to the next collision. A number of things happen when
a particle hits a surface. First, we need to check if a non-transmissive
boundary condition has been applied to the surface. If a vacuum boundary
condition has been applied, the particle is killed and any surface current
tallies are scored to as needed. If a reflective boundary condition has been
applied to the surface, surface current tallies are scored to and then the
particle's direction is changed according to the procedure in :ref:`reflection`.
Note that the white boundary condition can be considered as the special case of
reflective boundary condition, where the same processing method will be applied to
deal with the surface current tallies scoring, except for determining the
changes of particle's direction according to the procedures in :ref:`white`.
Next, we need to determine what cell is beyond the surface in the direction of
travel of the particle so that we can evaluate cross sections based on its
material properties. At initialization, a list of neighboring cells is created
for each surface in the problem as described in :ref:`neighbor-lists`. The
algorithm outlined in :ref:`find-cell` is used to find a cell containing the
particle with one minor modification; rather than searching all cells in the
base universe, only the list of neighboring cells is searched. If this search is
unsuccessful, then a search is done over every cell in the base universe.
.. _neighbor-lists:
-----------------------
Building Neighbor Lists
-----------------------
After the geometry has been loaded and stored in memory from an input file,
OpenMC builds a list for each surface containing any cells that are bounded by
that surface in order to speed up processing of surface crossings. The algorithm
to build these lists is as follows. First, we loop over all cells in the
geometry and count up how many times each surface appears in a specification as
bounding a negative half-space and bounding a positive half-space. Two arrays
are then allocated for each surface, one that lists each cell that contains the
negative half-space of the surface and one that lists each cell that contains
the positive half-space of the surface. Another loop is performed over all cells
and the neighbor lists are populated for each surface.
.. _reflection:
------------------------------
Reflective Boundary Conditions
------------------------------
If the velocity of a particle is :math:`\mathbf{v}` and it crosses a surface of
the form :math:`f(x,y,z) = 0` with a reflective boundary condition, it can be
shown based on geometric arguments that the velocity vector will then become
.. math::
:label: reflection-v
\mathbf{v'} = \mathbf{v} - 2 (\mathbf{v} \cdot \hat{\mathbf{n}})
\hat{\mathbf{n}}
where :math:`\hat{\mathbf{n}}` is a unit vector normal to the surface at the
point of the surface crossing. The rationale for this can be understood by
noting that :math:`(\mathbf{v} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}` is the
projection of the velocity vector onto the normal vector. By subtracting two
times this projection, the velocity is reflected with respect to the surface
normal. Since the magnitude of the velocity of the particle will not change as
it undergoes reflection, we can work with the direction of the particle instead,
simplifying equation :eq:`reflection-v` to
.. math::
:label: reflection-omega
\mathbf{\Omega'} = \mathbf{\Omega} - 2 (\mathbf{\Omega} \cdot
\hat{\mathbf{n}}) \hat{\mathbf{n}}
where :math:`\mathbf{v} = || \mathbf{v} || \mathbf{\Omega}`. The direction of
the surface normal will be the gradient of the surface at the point of crossing,
i.e. :math:`\mathbf{n} = \nabla f(x,y,z)`. Substituting this into equation
:eq:`reflection-omega`, we get
.. math::
:label: reflection-omega-2
\mathbf{\Omega'} = \mathbf{\Omega} - \frac{2 ( \mathbf{\Omega} \cdot \nabla
f )}{|| \nabla f ||^2} \nabla f
If we write the initial and final directions in terms of their vector
components, :math:`\mathbf{\Omega} = (u,v,w)` and :math:`\mathbf{\Omega'} = (u',
v', w')`, this allows us to represent equation :eq:`reflection-omega` as a
series of equations:
.. math::
:label: reflection-system
u' = u - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2}
\frac{\partial f}{\partial x} \\
v' = v - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2}
\frac{\partial f}{\partial y} \\
w' = w - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2}
\frac{\partial f}{\partial z}
One can then use equation :eq:`reflection-system` to develop equations for
transforming a particle's direction given the equation of the surface.
Plane Perpendicular to an Axis
------------------------------
For a plane that is perpendicular to an axis, the rule for reflection is almost
so simple that no derivation is needed at all. Nevertheless, we will proceed
with the derivation to confirm that the rules of geometry agree with our
intuition. The gradient of the surface :math:`f(x,y,z) = x - x_0 = 0` is simply
:math:`\nabla f = (1, 0, 0)`. Note that this vector is already normalized,
i.e. :math:`|| \nabla f || = 1`. The second two equations in
:eq:`reflection-system` tell us that :math:`v` and :math:`w` do not change and
the first tell us that
.. math::
:label: reflection-xplane
u' = u - 2u = -u
We see that reflection for a plane perpendicular to an axis only entails
negating the directional cosine for that axis.
Generic Plane
-------------
A generic plane has the form :math:`f(x,y,z) = Ax + By + Cz - D = 0`. Thus, the
gradient to the surface is simply :math:`\nabla f = (A,B,C)` whose norm squared
is :math:`A^2 + B^2 + C^2`. This implies that
.. math::
:label: reflection-plane-constant
\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{2(Au +
Bv + Cw)}{A^2 + B^2 + C^2}
Substituting equation :eq:`reflection-plane-constant` into equation
:eq:`reflection-system` gives us the form of the solution. For example, the
x-component of the reflected direction will be
.. math::
:label: reflection-plane
u' = u - \frac{2A(Au + Bv + Cw)}{A^2 + B^2 + C^2}
Cylinder Parallel to an Axis
----------------------------
A cylinder parallel to, for example, the x-axis has the form :math:`f(x,y,z) =
(y - y_0)^2 + (z - z_0)^2 - R^2 = 0`. Thus, the gradient to the surface is
.. math::
:label: reflection-cylinder-grad
\nabla f = 2 \left ( \begin{array}{c} 0 \\ y - y_0 \\ z - z_0 \end{array}
\right ) = 2 \left ( \begin{array}{c} 0 \\ \bar{y} \\ \bar{z} \end{array}
\right )
where we have introduced the constants :math:`\bar{y}` and
:math:`\bar{z}`. Taking the square of the norm of the gradient, we find that
.. math::
:label: reflection-cylinder-norm
|| \nabla f ||^2 = 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2
This implies that
.. math::
:label: reflection-cylinder-constant
\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} =
\frac{\bar{y}v + \bar{z}w}{R^2}
Substituting equations :eq:`reflection-cylinder-constant` and
:eq:`reflection-cylinder-grad` into equation :eq:`reflection-system` gives us
the form of the solution. In this case, the x-component will not change. The y-
and z-components of the reflected direction will be
.. math::
:label: reflection-cylinder
v' = v - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{y}}{R^2} \\
w' = w - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{z}}{R^2}
Sphere
------
The surface equation for a sphere has the form :math:`f(x,y,z) = (x - x_0)^2 +
(y - y_0)^2 + (z - z_0)^2 - R^2 = 0`. Thus, the gradient to the surface is
.. math::
:label: reflection-sphere-grad
\nabla f = 2 \left ( \begin{array}{c} x - x_0 \\ y - y_0 \\ z - z_0
\end{array} \right ) = 2 \left ( \begin{array}{c} \bar{x} \\ \bar{y} \\
\bar{z} \end{array} \right )
where we have introduced the constants :math:`\bar{x}, \bar{y}, \bar{z}`. Taking
the square of the norm of the gradient, we find that
.. math::
:label: reflection-sphere-norm
|| \nabla f ||^2 = 4 \bar{x}^2 + 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2
This implies that
.. math::
:label: reflection-sphere-constant
\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} =
\frac{\bar{x}u + \bar{y}v + \bar{z}w}{R^2}
Substituting equations :eq:`reflection-sphere-constant` and
:eq:`reflection-sphere-grad` into equation :eq:`reflection-system` gives us the
form of the solution:
.. math::
:label: reflection-sphere
u' = u - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{x} }{R^2} \\
v' = v - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{y} }{R^2} \\
w' = w - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{z} }{R^2}
Cone Parallel to an Axis
------------------------
A cone parallel to, for example, the z-axis has the form :math:`f(x,y,z) = (x -
x_0)^2 + (y - y_0)^2 - R^2(z - z_0)^2 = 0`. Thus, the gradient to the surface is
.. math::
:label: reflection-cone-grad
\nabla f = 2 \left ( \begin{array}{c} x - x_0 \\ y - y_0 \\ -R^2(z - z_0)
\end{array} \right ) = 2 \left ( \begin{array}{c} \bar{x} \\ \bar{y} \\
-R^2\bar{z} \end{array} \right )
where we have introduced the constants :math:`\bar{x}`, :math:`\bar{y}`, and
:math:`\bar{z}`. Taking the square of the norm of the gradient, we find that
.. math::
:label: reflection-cone-norm
|| \nabla f ||^2 = 4 \bar{x}^2 + \bar{y}^2 + 4 R^4 \bar{z}^2 \\ = 4 R^2
\bar{z}^2 + 4 R^4 \bar{z}^2 \\ = 4 R^2 (1 + R^2) \bar{z}^2
This implies that
.. math::
:label: reflection-cone-constant
\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} =
\frac{\bar{x}u + \bar{y}v - R^2\bar{z}w}{R^2 (1 + R^2) \bar{z}^2}
Substituting equations :eq:`reflection-cone-constant` and
:eq:`reflection-cone-grad` into equation :eq:`reflection-system` gives us the
form of the solution:
.. math::
:label: reflection-cone
u' = u - \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w) \bar{x}}{R^2 (1 + R^2)
\bar{z}^2}
v' = v - \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w) \bar{y}}{R^2 (1 + R^2)
\bar{z}^2}
w' = w + \frac{2 (\bar{x}u + \bar{y}v - R^2\bar{z}w)}{R^2 (1 + R^2) \bar{z}}
General Quadric
---------------
A general quadric surface has the form :math:`f(x,y,z) = Ax^2 + By^2 + Cz^2 +
Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is
.. math::
:label: reflection-quadric-grad
\nabla f = \left ( \begin{array}{c} 2Ax + Dy + Fz + G \\ 2By + Dx + Ez + H
\\ 2Cz + Ey + Fx + J \end{array} \right ).
.. _white:
-------------------------
White Boundary Conditions
-------------------------
The `white boundary condition <https://doi.org/10.1016/j.anucene.2019.05.006>`_
is usually applied in deterministic codes, where the particle will hit the
surface and travel back with isotropic angular distribution. The change in
particle's direction is sampled from a cosine distribution instead of uniform.
Figure :num:`fig-cosine-dist` shows an example of cosine-distribution reflection
on the arbitrary surface relative to the surface normal.
.. _fig-cosine-dist:
.. figure:: ../_images/cosine-dist.png
:align: center
:figclass: align-center
Cosine-distribution reflection on an arbitrary surface.
The probability density function (pdf) for the reflected direction can be
expressed as follows,
.. math::
:label: white-reflection-pdf
f(\mu, \phi) d\mu d\phi = \frac{\mu}{\pi} d\mu d\phi = 2\mu d\mu \frac{d\phi}{2\pi}
where :math:`\mu = \cos \theta` is the cosine of the polar angle between
reflected direction and the normal to the surface; and :math:`\theta` is the
azimuthal angle in :math:`[0,2\pi]`. We can separate the multivariate
probability density into two separate univariate density functions, one for
the cosine of the polar angle,
.. math::
:label: white-reflection-cosine
f(\mu) = 2\mu
and one for the azimuthal angle,
.. math::
:label: white-reflection-uniform
f(\phi) = \frac{1}{2\pi}.
Each of these density functions can be sampled by analytical inversion of the
cumulative distribution distribution, resulting in the following sampling
scheme:
.. math::
:label: white-reflection-sqrt-prn
\mu = \sqrt{\xi_1} \\
\phi = 2\pi\xi_2
where :math:`\xi_1` and :math:`\xi_2` are uniform random numbers on
:math:`[0,1)`. With the sampled values of :math:`\mu` and :math:`\phi`, the
final reflected direction vector can be computed via rotation of the surface
normal using the equations from :ref:`transform-coordinates`. The white boundary
condition can be applied to any kind of surface, as long as the normal to the
surface is known as in :ref:`reflection`.
.. _constructive solid geometry: https://en.wikipedia.org/wiki/Constructive_solid_geometry
.. _surfaces: https://en.wikipedia.org/wiki/Surface
.. _MCNP: http://mcnp.lanl.gov
.. _Serpent: http://montecarlo.vtt.fi
.. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc

21
docs/source/methods/index.rst Normal file
View file

@ -0,0 +1,21 @@
.. _methods:
======================
Theory and Methodology
======================
.. toctree::
:numbered:
:maxdepth: 3
introduction
geometry
cross_sections
random_numbers
neutron_physics
photon_physics
tallies
eigenvalue
parallelization
cmfd
energy_deposition

147
docs/source/methods/introduction.rst Normal file
View file

@ -0,0 +1,147 @@
.. _methods_introduction:
============
Introduction
============
The physical process by which a population of particles evolves over time is
governed by a number of `probability distributions`_. For instance, given a
particle traveling through some material, there is a probability distribution
for the distance it will travel until its next collision (an exponential
distribution). Then, when it collides with a nucleus, there is an associated
probability of undergoing each possible reaction with that nucleus. While the
behavior of any single particle is unpredictable, the average behavior of a
large population of particles originating from the same source is well defined.
If the probability distributions that govern the transport of a particle are
known, the process of single particles randomly streaming and colliding with
nuclei can be simulated directly with computers using a technique known as
`Monte Carlo`_ simulation. If enough particles are simulated this way, the
average behavior can be determined to within arbitrarily small statistical
error, a fact guaranteed by the `central limit theorem`_. To be more precise,
the central limit theorem tells us that the variance of the sample mean of some
physical parameter being estimated with Monte Carlo will be inversely
proportional to the number of realizations, i.e. the number of particles we
simulate:
.. math::
\sigma^2 \propto \frac{1}{N}.
where :math:`\sigma^2` is the variance of the sample mean and :math:`N` is the
number of realizations.
------------------------
Overview of Program Flow
------------------------
OpenMC performs a Monte Carlo simulation one particle at a time -- at no point
is more than one particle being tracked on a single program instance. Before any
particles are tracked, the problem must be initialized. This involves the
following steps:
- Read input files and building data structures for the geometry, materials,
tallies, and other associated variables.
- Initialize the pseudorandom number generator.
- Read the contiuous-energy or multi-group cross section data specified in
the problem.
- If using a special energy grid treatment such as a union energy grid or
lethargy bins, that must be initialized as well in a continuous-energy
problem.
- In a multi-group problem, individual nuclide cross section information is
combined to produce material-specific cross section data.
- In a fixed source problem, source sites are sampled from the specified
source. In an eigenvalue problem, source sites are sampled from some
initial source distribution or from a source file. The source sites
consist of coordinates, a direction, and an energy.
Once initialization is complete, the actual transport simulation can
proceed. The life of a single particle will proceed as follows:
1. The particle's properties are initialized from a source site previously
sampled.
2. Based on the particle's coordinates, the current cell in which the particle
resides is determined.
3. The energy-dependent cross sections for the material that the particle is
currently in are determined. Note that this includes the total
cross section, which is not pre-calculated.
4. The distance to the nearest boundary of the particle's cell is determined
based on the bounding surfaces to the cell.
5. The distance to the next collision is sampled. If the total material
cross section is :math:`\Sigma_t`, this can be shown to be
.. math::
d = -\frac{\ln \xi}{\Sigma_t}
where :math:`\xi` is a `pseudorandom number`_ sampled from a uniform
distribution on :math:`[0,1)`.
6. If the distance to the nearest boundary is less than the distance to the next
collision, the particle is moved forward to this boundary. Then, the process
is repeated from step 2. If the distance to collision is closer than the
distance to the nearest boundary, then the particle will undergo a collision.
7. The material at the collision site may consist of multiple nuclides. First,
the nuclide with which the collision will happen is sampled based on the
total cross sections. If the total cross section of material :math:`i` is
:math:`\Sigma_{t,i}`, then the probability that any nuclide is sampled is
.. math::
P(i) = \frac{\Sigma_{t,i}}{\Sigma_t}.
Note that the above selection of collided nuclide only applies to
continuous-energy simulations as multi-group simulations use nuclide
data which has already been combined in to material-specific data.
8. Once the specific nuclide is sampled, the random samples a reaction for
that nuclide based on the microscopic cross sections. If the microscopic
cross section for some reaction :math:`x` is :math:`\sigma_x` and the total
microscopic cross section for the nuclide is :math:`\sigma_t`, then the
probability that reaction :math:`x` will occur is
.. math::
P(x) = \frac{\sigma_x}{\sigma_t}.
Since multi-group simulations use material-specific data, the above is
performed with those material multi-group cross sections (i.e.,
macroscopic cross sections for the material) instead of microscopic
cross sections for the nuclide).
9. If the sampled reaction is elastic or inelastic scattering, the outgoing
energy and angle is sampled from the appropriate distribution. In
continuous-energy simulation, reactions of type :math:`(n,xn)` are treated
as scattering and any additional particles which may be created are added
to a secondary particle bank to be tracked later. In a multi-group
simulation, this secondary bank is not used but the particle weight is
increased accordingly. The original particle then continues from step 3.
If the reaction is absorption or fission, the particle dies and if
necessary, fission sites are created and stored in the fission bank.
After all particles have been simulated, there are a few final tasks that must
be performed before the run is finished. This include the following:
- With the accumulated sum and sum of squares for each tally, the sample mean
and its variance is calculated.
- All tallies and other results are written to disk.
- If requested, a source file is written to disk.
- Dynamically-allocated memory should be freed.
.. _probability distributions: https://en.wikipedia.org/wiki/Probability_distribution
.. _Monte Carlo: https://en.wikipedia.org/wiki/Monte_Carlo_method
.. _central limit theorem: https://en.wikipedia.org/wiki/Central_limit_theorem
.. _pseudorandom number: https://en.wikipedia.org/wiki/Pseudorandom_number_generator

1677
docs/source/methods/neutron_physics.rst Normal file
View file

File diff suppressed because it is too large Load diff

650
docs/source/methods/parallelization.rst Normal file
View file

@ -0,0 +1,650 @@
.. _methods_parallel:
===============
Parallelization
===============
Due to the computationally-intensive nature of Monte Carlo methods, there has
been an ever-present interest in parallelizing such simulations. Even in the
`first paper`_ on the Monte Carlo method, John Metropolis and Stanislaw Ulam
recognized that solving the Boltzmann equation with the Monte Carlo method could
be done in parallel very easily whereas the deterministic counterparts for
solving the Boltzmann equation did not offer such a natural means of
parallelism. With the introduction of `vector computers`_ in the early 1970s,
general-purpose parallel computing became a reality. In 1972, Troubetzkoy et
al. designed a Monte Carlo code to be run on the first vector computer, the
ILLIAC-IV [Troubetzkoy]_. The general principles from that work were later
refined and extended greatly through the `work of Forrest Brown`_ in the
1980s. However, as Brown's work shows, the `single-instruction multiple-data`_
(SIMD) parallel model inherent to vector processing does not lend itself to the
parallelism on particles in Monte Carlo simulations. Troubetzkoy et
al. recognized this, remarking that "the order and the nature of these physical
events have little, if any, correlation from history to history," and thus
following independent particle histories simultaneously using a SIMD model is
difficult.
The difficulties with vector processing of Monte Carlo codes led to the adoption
of the `single program multiple data`_ (SPMD) technique for parallelization. In
this model, each different process tracks a particle independently of other
processes, and between fission source generations the processes communicate data
through a `message-passing interface`_. This means of parallelism was enabled by
the introduction of message-passing standards in the late 1980s and early 1990s
such as PVM_ and MPI_. The SPMD model proved much easier to use in practice and
took advantage of the inherent parallelism on particles rather than
instruction-level parallelism. As a result, it has since become ubiquitous for
Monte Carlo simulations of transport phenomena.
Thanks to the particle-level parallelism using SPMD techniques, extremely high
parallel efficiencies could be achieved in Monte Carlo codes. Until the last
decade, even the most demanding problems did not require transmitting large
amounts of data between processors, and thus the total amount of time spent on
communication was not significant compared to the amount of time spent on
computation. However, today's computing power has created a demand for
increasingly large and complex problems, requiring a greater number of particles
to obtain decent statistics (and convergence in the case of criticality
calculations). This results in a correspondingly higher amount of communication,
potentially degrading the parallel efficiency. Thus, while Monte Carlo
simulations may seem `embarrassingly parallel`_, obtaining good parallel scaling
with large numbers of processors can be quite difficult to achieve in practice.
.. _fission-bank-algorithms:
-----------------------
Fission Bank Algorithms
-----------------------
Master-Slave Algorithm
----------------------
Monte Carlo particle transport codes commonly implement a SPMD model by having
one master process that controls the scheduling of work and the remaining
processes wait to receive work from the master, process the work, and then send
their results to the master at the end of the simulation (or a source iteration
in the case of an eigenvalue calculation). This idea is illustrated in
:ref:`figure-master-slave`.
.. _figure-master-slave:
.. figure:: ../_images/master-slave.png
:align: center
:figclass: align-center
Communication pattern in master-slave algorithm.
Eigenvalue calculations are slightly more difficult to parallelize than fixed
source calculations since it is necessary to converge on the fission source
distribution and eigenvalue before tallying. In the
:ref:`method-successive-generations`, to ensure that the results are
reproducible, one must guarantee that the process by which fission sites are
randomly sampled does not depend on the number of processors. What is typically
done is the following:
1. Each compute node sends_ its fission bank sites to a master process;
2. The master process sorts or orders the fission sites based on a unique
identifier;
3. The master process samples :math:`N` fission sites from the ordered array
of :math:`M` sites; and
4. The master process broadcasts_ all the fission sites to the compute
nodes.
The first and last steps of this process are the major sources of communication
overhead between cycles. Since the master process must receive :math:`M` fission
sites from the compute nodes, the first step is necessarily serial. This step
can be completed in :math:`O(M)` time. The broadcast step can benefit from
parallelization through a tree-based algorithm. Despite this, the communication
overhead is still considerable.
To see why this is the case, it is instructive to look at a hypothetical
example. Suppose that a calculation is run with :math:`N = 10,000,000` neutrons
across 64 compute nodes. On average, :math:`M = 10,000,000` fission sites will
be produced. If the data for each fission site consists of a spatial location
(three 8 byte real numbers) and a unique identifier (one 4 byte integer), the
memory required per site is 28 bytes. To broadcast 10,000,000 source sites to 64
nodes will thus require transferring 17.92 GB of data. Since each compute node
does not need to keep every source site in memory, one could modify the
algorithm from a broadcast to a scatter_. However, for practical reasons
(e.g. work self-scheduling), this is normally not done in production Monte Carlo
codes.
.. _nearest-neighbors-algorithm:
Nearest Neighbors Algorithm
---------------------------
To reduce the amount of communication required in a fission bank synchronization
algorithm, it is desirable to move away from the typical master-slave algorithm
to an algorithm whereby the compute nodes communicate with one another only as
needed. This concept is illustrated in :ref:`figure-nearest-neighbor`.
.. _figure-nearest-neighbor:
.. figure:: ../_images/nearest-neighbor.png
:align: center
:figclass: align-center
Communication pattern in nearest neighbor algorithm.
Since the source sites for each cycle are sampled from the fission sites banked
from the previous cycle, it is a common occurrence for a fission site to be
banked on one compute node and sent back to the master only to get sent back to
the same compute node as a source site. As a result, much of the communication
inherent in the algorithm described previously is entirely unnecessary. By
keeping the fission sites local, having each compute node sample fission sites,
and sending sites between nodes only as needed, one can cut down on most of the
communication. One algorithm to achieve this is as follows:
1. An exclusive scan is performed on the number of sites banked, and the
total number of fission bank sites is broadcasted to all compute nodes. By
picturing the fission bank as one large array distributed across multiple
nodes, one can see that this step enables each compute node to determine the
starting index of fission bank sites in this array. Let us call the starting
and ending indices on the :math:`i`-th node :math:`a_i` and :math:`b_i`,
respectively;
2. Each compute node samples sites at random from the fission bank using the
same starting seed. A separate array on each compute node is created that
consists of sites that were sampled local to that node, i.e. if the index of
the sampled site is between :math:`a_i` and :math:`b_i`, it is set aside;
3. If any node sampled more than :math:`N/p` fission sites where :math:`p`
is the number of compute nodes, the extra sites are put in a separate array
and sent to all other compute nodes. This can be done efficiently using the
allgather_ collective operation;
4. The extra sites are divided among those compute nodes that sampled fewer
than :math:`N/p` fission sites.
However, even this algorithm exhibits more communication than necessary since
the allgather will send fission bank sites to nodes that don't necessarily
need any extra sites.
One alternative is to replace the allgather with a series of sends. If
:math:`a_i` is less than :math:`iN/p`, then send :math:`iN/p - a_i` sites to the
left adjacent node. Similarly, if :math:`a_i` is greater than :math:`iN/p`, then
receive :math:`a_i - iN/p` from the left adjacent node. This idea is applied to
the fission bank sites at the end of each node's array as well. If :math:`b_i`
is less than :math:`(i+1)N/p`, then receive :math:`(i+1)N/p - b_i` sites from
the right adjacent node. If :math:`b_i` is greater than :math:`(i+1)N/p`, then
send :math:`b_i - (i+1)N/p` sites to the right adjacent node. Thus, each compute
node sends/receives only two messages under normal circumstances.
The following example illustrates how this algorithm works. Let us suppose we
are simulating :math:`N = 1000` neutrons across four compute nodes. For this
example, it is instructive to look at the state of the fission bank and source
bank at several points in the algorithm:
1. The beginning of a cycle where each node has :math:`N/p` source sites;
2. The end of a cycle where each node has accumulated fission sites;
3. After sampling, where each node has some amount of source sites usually
not equal to :math:`N/p`;
4. After redistribution, each node again has :math:`N/p` source sites for
the next cycle;
At the end of each cycle, each compute node needs 250 fission bank sites to
continue on the next cycle. Let us suppose that :math:`p_0` produces 270 fission
banks sites, :math:`p_1` produces 230, :math:`p_2` produces 290, and :math:`p_3`
produces 250. After each node samples from its fission bank sites, let's assume
that :math:`p_0` has 260 source sites, :math:`p_1` has 215, :math:`p_2` has 280,
and :math:`p_3` has 245. Note that the total number of sampled sites is 1000 as
needed. For each node to have the same number of source sites, :math:`p_0` needs
to send its right-most 10 sites to :math:`p_1`, and :math:`p_2` needs to send
its left-most 25 sites to :math:`p_1` and its right-most 5 sites to
:math:`p_3`. A schematic of this example is shown in
:ref:`figure-neighbor-example`. The data local to each node is given a different
hatching, and the cross-hatched regions represent source sites that are
communicated between adjacent nodes.
.. _figure-neighbor-example:
.. figure:: ../_images/nearest-neighbor-example.png
:align: center
:figclass: align-center
Example of nearest neighbor algorithm.
.. _master-slave-cost:
Cost of Master-Slave Algorithm
------------------------------
While the prior considerations may make it readily apparent that the novel
algorithm should outperform the traditional algorithm, it is instructive to look
at the total communication cost of the novel algorithm relative to the
traditional algorithm. This is especially so because the novel algorithm does
not have a constant communication cost due to stochastic fluctuations. Let us
begin by looking at the cost of communication in the traditional algorithm
As discussed earlier, the traditional algorithm is composed of a series of sends
and typically a broadcast. To estimate the communication cost of the algorithm,
we can apply a simple model that captures the essential features. In this model,
we assume that the time that it takes to send a message between two nodes is
given by :math:`\alpha + (sN)\beta`, where :math:`\alpha` is the time it takes
to initiate the communication (commonly called the latency_), :math:`\beta` is
the transfer time per unit of data (commonly called the bandwidth_), :math:`N`
is the number of fission sites, and :math:`s` is the size in bytes of each
fission site.
The first step of the traditional algorithm is to send :math:`p` messages to the
master node, each of size :math:`sN/p`. Thus, the total time to send these
messages is
.. math::
:label: t-send
t_{\text{send}} = p\alpha + sN\beta.
Generally, the best parallel performance is achieved in a weak scaling scheme
where the total number of histories is proportional to the number of
processors. However, we see that when :math:`N` is proportional to :math:`p`,
the time to send these messages increases proportionally with :math:`p`.
Estimating the time of the broadcast is complicated by the fact that different
MPI implementations may use different algorithms to perform collective
communications. Worse yet, a single implementation may use a different algorithm
depending on how many nodes are communicating and the size of the message. Using
multiple algorithms allows one to minimize latency for small messages and
minimize bandwidth for long messages.
We will focus here on the implementation of broadcast in the MPICH2_
implementation. For short messages, MPICH2 uses a `binomial tree`_ algorithm. In
this algorithm, the root process sends the data to one node in the first step,
and then in the subsequent, both the root and the other node can send the data
to other nodes. Thus, it takes a total of :math:`\lceil \log_2 p \rceil` steps
to complete the communication. The time to complete the communication is
.. math::
:label: t-short
t_{\text{short}} = \lceil \log_2 p \rceil \left ( \alpha + sN\beta \right ).
This algorithm works well for short messages since the latency term scales
logarithmically with the number of nodes. However, for long messages, an
algorithm that has lower bandwidth has been proposed by Barnett_ and implemented
in MPICH2. Rather than using a binomial tree, the broadcast is divided into a
scatter and an allgather. The time to complete the scatter is :math:` \log_2 p
\: \alpha + \frac{p-1}{p} N\beta` using a binomial tree algorithm. The allgather
is performed using a ring algorithm that completes in :math:`p-1) \alpha +
\frac{p-1}{p} N\beta`. Thus, together the time to complete the broadcast is
.. math::
:label: t-broadcast
t_{\text{long}} = \left ( \log_2 p + p - 1 \right ) \alpha + 2 \frac{p-1}{p}
sN\beta.
The fission bank data will generally exceed the threshold for switching from
short to long messages (typically 8 kilobytes), and thus we will use the
equation for long messages. Adding equations :eq:`t-send` and :eq:`t-broadcast`,
the total cost of the series of sends and the broadcast is
.. math::
:label: t-old
t_{\text{old}} = \left ( \log_2 p + 2p - 1 \right ) \alpha + \frac{3p-2}{p}
sN\beta.
Cost of Nearest Neighbor Algorithm
----------------------------------
With the communication cost of the traditional fission bank algorithm
quantified, we now proceed to discuss the communicatin cost of the proposed
algorithm. Comparing the cost of communication of this algorithm with the
traditional algorithm is not trivial due to fact that the cost will be a
function of how many fission sites are sampled on each node. If each node
samples exactly :math:`N/p` sites, there will not be communication between nodes
at all. However, if any one node samples more or less than :math:`N/p` sites,
the deviation will result in communication between logically adjacent nodes. To
determine the expected deviation, one can analyze the process based on the
fundamentals of the Monte Carlo process.
The steady-state neutron transport equation for a multiplying medium can be
written in the form of an eigenvalue problem,
.. math::
:label: NTE
S(\mathbf{r})= \frac{1}{k} \int F(\mathbf{r}' \rightarrow
\mathbf{r})S(\mathbf{r}')\: d\mathbf{r},
where :math:`\mathbf{r}` is the spatial coordinates of the neutron,
:math:`S(\mathbf{r})` is the source distribution defined as the expected number
of neutrons born from fission per unit phase-space volume at :math:`\mathbf{r}`,
:math:`F( \mathbf{r}' \rightarrow \mathbf{r})` is the expected number of
neutrons born from fission per unit phase space volume at :math:`\mathbf{r}`
caused by a neutron at :math:`\mathbf{r}`, and :math:`k` is the eigenvalue. The
fundamental eigenvalue of equation :eq:`NTE` is known as :math:`k_{eff}`, but
for simplicity we will simply refer to it as :math:`k`.
In a Monte Carlo criticality simulation, the power iteration method is applied
iteratively to obtain stochastic realizations of the source distribution and
estimates of the :math:`k`-eigenvalue. Let us define :math:`\hat{S}^{(m)}` to be
the realization of the source distribution at cycle :math:`m` and
:math:`\hat{\epsilon}^{(m)}` be the noise arising from the stochastic nature of
the tracking process. We can write the stochastic realization in terms of the
fundamental source distribution and the noise component as (see `Brissenden and
Garlick`_):
.. math::
:label: source
\hat{S}^{(m)}(\mathbf{r})= N S(\mathbf{r}) + \sqrt{N}
\hat{\epsilon}^{(m)}(\mathbf{r}),
where :math:`N` is the number of particle histories per cycle. Without loss of
generality, we shall drop the superscript notation indicating the cycle as it is
understood that the stochastic realization is at a particular cycle. The
expected value of the stochastic source distribution is simply
.. math::
:label: expected-value-source
E \left[ \hat{S}(\mathbf{r})\right] = N S (\mathbf{r})
since :math:`E \left[ \hat{\epsilon}(\mathbf{r})\right] = 0`. The noise in the
source distribution is due only to :math:`\hat{\epsilon}(\mathbf{r})` and thus
the variance of the source distribution will be
.. math::
:label: var-source
\text{Var} \left[ \hat{S}(\mathbf{r})\right] = N \text{Var} \left[
\hat{\epsilon}(\mathbf{r}) \right].
Lastly, the stochastic and true eigenvalues can be written as integrals over all
phase space of the stochastic and true source distributions, respectively, as
.. math::
:label: k-to-source
\hat{k} = \frac{1}{N} \int \hat{S}(\mathbf{r}) \: d\mathbf{r} \quad
\text{and} \quad k = \int S(\mathbf{r}) \: d\mathbf{r},
noting that :math:`S(\mathbf{r})` is :math:`O(1)`. One should note that the
expected value :math:`k` calculated by Monte Carlo power iteration (i.e. the
method of successive generations) will be biased from the true fundamental
eigenvalue of equation :eq:`NTE` by :math:`O(1/N)` (see `Brissenden and
Garlick`_), but we will assume henceforth that the number of particle histories
per cycle is sufficiently large to neglect this bias.
With this formalism, we now have a framework within which we can determine the
properties of the distribution of expected number of fission sites. The explicit
form of the source distribution can be written as
.. math::
:label: source-explicit
\hat{S}(\mathbf{r}) = \sum_{i=1}^{M} w_i \delta( \mathbf{r} - \mathbf{r}_i )
where :math:`\mathbf{r}_i` is the spatial location of the :math:`i`-th fission
site, :math:`w_i` is the statistical weight of the fission site at
:math:`\mathbf{r}_i`, and :math:`M` is the total number of fission sites. It is
clear that the total weight of the fission sites is simply the integral of the
source distribution. Integrating equation :eq:`source` over all space, we obtain
.. math::
:label: source-integrated
\int \hat{S}(\mathbf{r}) \: d\mathbf{r} = N \int S(\mathbf{r}) \:
d\mathbf{r} + \sqrt{N} \int \hat{\epsilon}(\mathbf{r}) \: d\mathbf{r} .
Substituting the expressions for the stochastic and true eigenvalues from
equation :eq:`k-to-source`, we can relate the stochastic eigenvalue to the
integral of the noise component of the source distribution as
.. math::
:label: noise-integeral
N\hat{k} = Nk + \sqrt{N} \int \hat{\epsilon}(\mathbf{r}) \: d\mathbf{r}.
Since the expected value of :math:`\hat{\epsilon}` is zero, the expected value
of its integral will also be zero. We thus see that the variance of the integral
of the source distribution, i.e. the variance of the total weight of fission
sites produced, is directly proportional to the variance of the integral of the
noise component. Let us call this term :math:`\sigma^2` for simplicity:
.. math::
:label: variance-sigma2
\text{Var} \left[ \int \hat{S}(\mathbf{r}) \right ] = N \sigma^2.
The actual value of :math:`\sigma^2` will depend on the physical nature of the
problem, whether variance reduction techniques are employed, etc. For instance,
one could surmise that for a highly scattering problem, :math:`\sigma^2` would
be smaller than for a highly absorbing problem since more collisions will lead
to a more precise estimate of the source distribution. Similarly, using implicit
capture should in theory reduce the value of :math:`\sigma^2`.
Let us now consider the case where the :math:`N` total histories are divided up
evenly across :math:`p` compute nodes. Since each node simulates :math:`N/p`
histories, we can write the source distribution as
.. math::
:label: source-node
\hat{S}_i(\mathbf{r})= \frac{N}{p} S(\mathbf{r}) + \sqrt{\frac{N}{p}}
\hat{\epsilon}_i(\mathbf{r}) \quad \text{for} \quad i = 1, \dots, p
Integrating over all space and simplifying, we can obtain an expression for the
eigenvalue on the :math:`i`-th node:
.. math::
:label: k-i-hat
\hat{k}_i = k + \sqrt{\frac{p}{N}} \int \hat{\epsilon}_i(\mathbf{r}) \:
d\mathbf{r}.
It is easy to show from this expression that the stochastic realization of the
global eigenvalue is merely the average of these local eigenvalues:
.. math::
:label: average-k-as-sum
\hat{k} = \frac{1}{p} \sum_{i=1}^p \hat{k}_i.
As was mentioned earlier, at the end of each cycle one must sample :math:`N`
sites from the :math:`M` sites that were created. Thus, the source for the next
cycle can be seen as the fission source from the current cycle divided by the
stochastic realization of the eigenvalue since it is clear from equation
:eq:`k-to-source` that :math:`\hat{k} = M/N`. Similarly, the number of sites
sampled on each compute node that will be used for the next cycle is
.. math::
:label: sites-per-node
M_i = \frac{1}{\hat{k}} \int \hat{S}_i(\mathbf{r}) \: d\mathbf{r} =
\frac{N}{p} \frac{\hat{k}_i}{\hat{k}}.
While we know conceptually that each compute node will under normal
circumstances send two messages, many of these messages will overlap. Rather
than trying to determine the actual communication cost, we will instead attempt
to determine the maximum amount of data being communicated from one node to
another. At any given cycle, the number of fission sites that the :math:`j`-th
compute node will send or receive (:math:`\Lambda_j`) is
.. math::
:label: Lambda
\Lambda_j = \left | \sum_{i=1}^j M_i - \frac{jN}{p} \right |.
Noting that :math:`jN/p` is the expected value of the summation, we can write
the expected value of :math:`\Lambda_j` as the mean absolute deviation of the
summation:
.. math::
:label: mean-dev-lambda
E \left [ \Lambda_j \right ] = E \left [ \left | \sum_{i=1}^j M_i -
\frac{jN}{p} \right | \right ] = \text{MD} \left [ \sum_{i=1}^j M_i \right ]
where :math:`\text{MD}` indicates the mean absolute deviation of a random
variable. The mean absolute deviation is an alternative measure of variability.
In order to ascertain any information about the mean deviation of :math:`M_i`,
we need to know the nature of its distribution. Thus far, we have said nothing
of the distributions of the random variables in question. The total number of
fission sites resulting from the tracking of :math:`N` neutrons can be shown to
be normally distributed via the :ref:`central-limit-theorem` (provided that
:math:`N` is sufficiently large) since the fission sites resulting from each
neutron are "sampled" from independent, identically-distributed random
variables. Thus, :math:`\hat{k}` and :math:`\int \hat{S} (\mathbf{r}) \:
d\mathbf{r}` will be normally distributed as will the individual estimates of
these on each compute node.
Next, we need to know what the distribution of :math:`M_i` in equation
:eq:`sites-per-node` is or, equivalently, how :math:`\hat{k}_i / \hat{k}` is
distributed. The distribution of a ratio of random variables is not easy to
calculate analytically, and it is not guaranteed that the ratio distribution is
normal if the numerator and denominator are normally distributed. For example,
if :math:`X` is a standard normal distribution and :math:`Y` is also standard
normal distribution, then the ratio :math:`X/Y` has the standard `Cauchy
distribution`_. The reader should be reminded that the Cauchy distribution has
no defined mean or variance. That being said, Geary_ has shown that, for the
case of two normal distributions, if the denominator is unlikely to assume
values less than zero, then the ratio distribution is indeed approximately
normal. In our case, :math:`\hat{k}` absolutely cannot assume a value less than
zero, so we can be reasonably assured that the distribution of :math:`M_i` will
be normal.
For a normal distribution with mean :math:`\mu` and distribution function
:math:`f(x)`, it can be shown that
.. math::
:label: mean-dev-to-stdev
\int_{-\infty}^{\infty} f(x) \left | x - \mu \right | \: dx =
\sqrt{\frac{2}{\pi} \int_{-\infty}^{\infty} f(x) \left ( x - \mu \right )^2
\: dx}
and thus the mean absolute deviation is :math:`\sqrt{2/\pi}` times the standard
deviation. Therefore, to evaluate the mean absolute deviation of :math:`M_i`, we
need to first determine its variance. Substituting equation
:eq:`average-k-as-sum` into equation :eq:`sites-per-node`, we can rewrite
:math:`M_i` solely in terms of :math:`\hat{k}_1, \dots, \hat{k}_p`:
.. math::
:label: M-i
M_i = \frac{N \hat{k}_i}{\sum\limits_{j=1}^p \hat{k}_j}.
Since we know the variance of :math:`\hat{k}_i`, we can use the error
propagation law to determine the variance of :math:`M_i`:
.. math::
:label: M-variance
\text{Var} \left [ M_i \right ] = \sum_{j=1}^p \left ( \frac{\partial
M_i}{\partial \hat{k}_j} \right )^2 \text{Var} \left [ \hat{k}_j \right ] +
\sum\limits_{j \neq m} \sum\limits_{m=1}^p \left ( \frac{\partial
M_i}{\partial \hat{k}_j} \right ) \left ( \frac{\partial M_i}{\partial
\hat{k}_m} \right ) \text{Cov} \left [ \hat{k}_j, \hat{k}_m \right ]
where the partial derivatives are evaluated at :math:`\hat{k}_j = k`. Since
:math:`\hat{k}_j` and :math:`\hat{k}_m` are independent if :math:`j \neq m`,
their covariance is zero and thus the second term cancels out. Evaluating the
partial derivatives, we obtain
.. math::
:label: M-variance-2
\text{Var} \left [ M_i \right ] = \left ( \frac{N(p-1)}{kp^2} \right )^2
\frac{p\sigma^2}{N} + \sum_{j \neq i} \left ( \frac{-N}{kp^2} \right )^2
\frac{p\sigma^2}{N} = \frac{N(p-1)}{k^2p^2} \sigma^2.
Through a similar analysis, one can show that the variance of
:math:`\sum_{i=1}^j M_i` is
.. math::
:label: sum-M-variance
\text{Var} \left [ \sum_{i=1}^j M_i \right ] = \frac{Nj(p-j)}{k^2p^2}
\sigma^2
Thus, the expected amount of communication on node :math:`j`, i.e. the mean
absolute deviation of :math:`\sum_{i=1}^j M_i` is proportional to
.. math::
:label: communication-cost
E \left [ \Lambda_j \right ] = \sqrt{\frac{2Nj(p-j)\sigma^2}{\pi k^2p^2}}.
This formula has all the properties that one would expect based on intuition:
1. As the number of histories increases, the communication cost on each node
increases as well;
2. If :math:`p=1`, i.e. if the problem is run on only one compute node, the
variance will be zero. This reflects the fact that exactly :math:`N` sites
will be sampled if there is only one node.
3. For :math:`j=p`, the variance will be zero. Again, this says that when
you sum the number of sites from each node, you will get exactly :math:`N`
sites.
We can determine the node that has the highest communication cost by
differentiating equation :eq:`communication-cost` with respect to :math:`j`,
setting it equal to zero, and solving for :math:`j`. Doing so yields
:math:`j_{\text{max}} = p/2`. Interestingly, substituting :math:`j = p/2` in
equation :eq:`communication-cost` shows us that the maximum communication cost
is actually independent of the number of nodes:
.. math::
:label: maximum-communication
E \left [ \Lambda_{j_{\text{max}}} \right ] = \sqrt{ \frac{N\sigma^2}{2\pi
k^2}}.
.. only:: html
.. rubric:: References
.. [Troubetzkoy] E. Troubetzkoy, H. Steinberg, and M. Kalos, "Monte Carlo
Radiation Penetration Calculations on a Parallel Computer,"
*Trans. Am. Nucl. Soc.*, **17**, 260 (1973).
.. _first paper: https://doi.org/10.2307/2280232
.. _work of Forrest Brown: http://hdl.handle.net/2027.42/24996
.. _Brissenden and Garlick: https://doi.org/10.1016/0306-4549(86)90095-2
.. _MPICH2: http://www.mcs.anl.gov/mpi/mpich
.. _binomial tree: https://www.mcs.anl.gov/~thakur/papers/ijhpca-coll.pdf
.. _Geary: https://doi.org/10.2307/2342070
.. _Barnett: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.51.7772
.. _single-instruction multiple-data: https://en.wikipedia.org/wiki/SIMD
.. _vector computers: https://en.wikipedia.org/wiki/Vector_processor
.. _single program multiple data: https://en.wikipedia.org/wiki/SPMD
.. _message-passing interface: https://en.wikipedia.org/wiki/Message_Passing_Interface
.. _PVM: http://www.csm.ornl.gov/pvm/pvm_home.html
.. _MPI: http://www.mcs.anl.gov/research/projects/mpi/
.. _embarrassingly parallel: https://en.wikipedia.org/wiki/Embarrassingly_parallel
.. _sends: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Send.html
.. _broadcasts: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Bcast.html
.. _scatter: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Scatter.html
.. _allgather: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Allgather.html
.. _Cauchy distribution: https://en.wikipedia.org/wiki/Cauchy_distribution
.. _latency: https://en.wikipedia.org/wiki/Latency_(engineering)#Packet-switched_networks
.. _bandwidth: https://en.wikipedia.org/wiki/Bandwidth_(computing)

1072
docs/source/methods/photon_physics.rst Normal file
View file

File diff suppressed because it is too large Load diff

72
docs/source/methods/random_numbers.rst Normal file
View file

@ -0,0 +1,72 @@
.. _methods_random_numbers:
========================
Random Number Generation
========================
In order to sample probability distributions, one must be able to produce random
numbers. The standard technique to do this is to generate numbers on the
interval :math:`[0,1)` from a deterministic sequence that has properties that
make it appear to be random, e.g. being uniformly distributed and not exhibiting
correlation between successive terms. Since the numbers produced this way are
not truly "random" in a strict sense, they are typically referred to as
pseudorandom numbers, and the techniques used to generate them are pseudorandom
number generators (PRNGs). Numbers sampled on the unit interval can then be
transformed for the purpose of sampling other continuous or discrete probability
distributions.
------------------------------
Linear Congruential Generators
------------------------------
There are a great number of algorithms for generating random numbers. One of the
simplest and commonly used algorithms is called a `linear congruential
generator`_. We start with a random number *seed* :math:`\xi_0` and a sequence
of random numbers can then be generated using the following recurrence relation:
.. math::
:label: lcg
\xi_{i+1} = g \xi_i + c \mod M
where :math:`g`, :math:`c`, and :math:`M` are constants. The choice of these
constants will have a profound effect on the quality and performance of the
generator, so they should not be chosen arbitrarily. As Donald Knuth stated in
his seminal work *The Art of Computer Programming*, "random numbers should not
be generated with a method chosen at random. Some theory should be used."
Typically, :math:`M` is chosen to be a power of two as this enables :math:`x
\mod M` to be performed using the bitwise AND operator with a bit mask. The
constants for the linear congruential generator used by default in OpenMC are
:math:`g = 2806196910506780709`, :math:`c = 1`, and :math:`M = 2^{63}` (see
`L'Ecuyer`_).
Skip-ahead Capability
---------------------
One of the important capabilities for a random number generator is to be able to
skip ahead in the sequence of random numbers. Without this capability, it would
be very difficult to maintain reproducibility in a parallel calculation. If we
want to skip ahead :math:`N` random numbers and :math:`N` is large, the cost of
sampling :math:`N` random numbers to get to that position may be prohibitively
expensive. Fortunately, algorithms have been developed that allow us to skip
ahead in :math:`O(\log_2 N)` operations instead of :math:`O(N)`. One algorithm
to do so is described in a paper by Brown_. This algorithm relies on the following
relationship:
.. math::
:label: lcg-skipahead
\xi_{i+k} = g^k \xi_i + c \frac{g^k - 1}{g - 1} \mod M
Note that equation :eq:`lcg-skipahead` has the same general form as equation :eq:`lcg`, so
the idea is to determine the new multiplicative and additive constants in
:math:`O(\log_2 N)` operations.
.. only:: html
.. rubric:: References
.. _L'Ecuyer: https://doi.org/10.1090/S0025-5718-99-00996-5
.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl-rn-arb-stride.pdf
.. _linear congruential generator: https://en.wikipedia.org/wiki/Linear_congruential_generator

515
docs/source/methods/tallies.rst Normal file
View file

@ -0,0 +1,515 @@
.. _methods_tallies:
=======
Tallies
=======
Note that the methods discussed in this section are written specifically for
continuous-energy mode but equivalent apply to the multi-group mode if the
particle's energy is replaced with the particle's group
------------------
Filters and Scores
------------------
The tally capability in OpenMC takes a similar philosophy as that employed in
the MC21_ Monte Carlo code to give maximum flexibility in specifying tallies
while still maintaining scalability. Any tally in a Monte Carlo simulation can
be written in the following form:
.. math::
:label: tally-integral
X = \underbrace{\int d\mathbf{r} \int d\mathbf{\Omega} \int
dE}_{\text{filters}} \underbrace{f(\mathbf{r}, \mathbf{\Omega},
E)}_{\text{scores}} \psi (\mathbf{r}, \mathbf{\Omega}, E)
A user can specify one or more filters which identify which regions of phase
space should score to a given tally (the limits of integration as shown in
equation :eq:`tally-integral`) as well as the scoring function (:math:`f` in
equation :eq:`tally-integral`). For example, if the desired tally was the
:math:`(n,\gamma)` reaction rate in a fuel pin, the filter would specify the
cell which contains the fuel pin and the scoring function would be the radiative
capture macroscopic cross section. The following quantities can be scored in
OpenMC: flux, total reaction rate, scattering reaction rate, neutron production
from scattering, higher scattering moments, :math:`(n,xn)` reaction rates,
absorption reaction rate, fission reaction rate, neutron production rate from
fission, and surface currents. The following variables can be used as filters:
universe, material, cell, birth cell, surface, mesh, pre-collision energy,
post-collision energy, polar angle, azimuthal angle, and the cosine of the
change-in-angle due to a scattering event.
With filters for pre- and post-collision energy and scoring functions for
scattering and fission production, it is possible to use OpenMC to generate
cross sections with user-defined group structures. These multigroup cross
sections can subsequently be used in deterministic solvers such as coarse mesh
finite difference (CMFD) diffusion.
------------------------------
Using Maps for Filter-Matching
------------------------------
Some Monte Carlo codes suffer severe performance penalties when tallying a large
number of quantities. Care must be taken to ensure that a tally system scales
well with the total number of tally bins. In OpenMC, a mapping technique is used
that allows for a fast determination of what tally/bin combinations need to be
scored to a given particle's phase space coordinates. For each discrete filter
variable, a list is stored that contains the tally/bin combinations that could
be scored to for each value of the filter variable. If a particle is in cell
:math:`n`, the mapping would identify what tally/bin combinations specify cell
:math:`n` for the cell filter variable. In this manner, it is not necessary to
check the phase space variables against each tally. Note that this technique
only applies to discrete filter variables and cannot be applied to energy,
angle, or change-in-angle bins. For these filters, it is necessary to perform
a binary search on the specified energy grid.
-----------------------------------------
Volume-Integrated Flux and Reaction Rates
-----------------------------------------
One quantity we may wish to compute during the course of a Monte Carlo
simulation is the flux or a reaction rate integrated over a finite volume. The
volume may be a particular cell, a collection of cells, or the entire
geometry. There are various methods by which we can estimate reaction rates
Analog Estimator
----------------
The analog estimator is the simplest type of estimator for reaction rates. The
basic idea is that we simply count the number of actual reactions that take
place and use that as our estimate for the reaction rate. This can be written
mathematically as
.. math::
:label: analog-estimator
R_x = \frac{1}{W} \sum_{i \in A} w_i
where :math:`R_x` is the reaction rate for reaction :math:`x`, :math:`i` denotes
an index for each event, :math:`A` is the set of all events resulting in
reaction :math:`x`, and :math:`W` is the total starting weight of the particles,
and :math:`w_i` is the pre-collision weight of the particle as it enters event
:math:`i`. One should note that equation :eq:`analog-estimator` is
volume-integrated so if we want a volume-averaged quantity, we need to divided
by the volume of the region of integration. If survival biasing is employed, the
analog estimator cannot be used for any reactions with zero neutrons in the exit
channel.
Collision Estimator
-------------------
While the analog estimator is conceptually very simple and easy to implement, it
can suffer higher variance due to the fact low probability events will not occur
often enough to get good statistics if they are being tallied. Thus, it is
desirable to use a different estimator that allows us to score to the tally more
often. One such estimator is the collision estimator. Instead of tallying a
reaction only when it happens, the idea is to make a contribution to the tally
at every collision.
We can start by writing a formula for the collision estimate of the flux. Since
:math:`R = \Sigma_t \phi` where :math:`R` is the total reaction rate,
:math:`\Sigma_t` is the total macroscopic cross section, and :math:`\phi` is the
scalar flux, it stands to reason that we can estimate the flux by taking an
estimate of the total reaction rate and dividing it by the total macroscopic
cross section. This gives us the following formula:
.. math::
:label: collision-estimator-flux
\phi = \frac{1}{W} \sum_{i \in C} \frac{w_i}{\Sigma_t (E_i)}
where :math:`W` is again the total starting weight of the particles, :math:`C`
is the set of all events resulting in a collision with a nucleus, and
:math:`\Sigma_t (E)` is the total macroscopic cross section of the target
material at the incoming energy of the particle :math:`E_i`.
If we multiply both sides of equation :eq:`collision-estimator-flux` by the
macroscopic cross section for some reaction :math:`x`, then we get the collision
estimate for the reaction rate for that reaction:
.. math::
:label: collision-estimator
R_x = \frac{1}{W} \sum_{i \in C} \frac{w_i \Sigma_x (E_i)}{\Sigma_t (E_i)}
where :math:`\Sigma_x (E_i)` is the macroscopic cross section for reaction
:math:`x` at the incoming energy of the particle :math:`E_i`. In comparison to
equation :eq:`analog-estimator`, we see that the collision estimate will result
in a tally with a larger number of events that score to it with smaller
contributions (since we have multiplied it by :math:`\Sigma_x / \Sigma_t`).
Track-length Estimator
----------------------
One other method we can use to increase the number of events that scores to
tallies is to use an estimator the scores contributions to a tally at every
track for the particle rather than every collision. This is known as a
track-length estimator, sometimes also called a path-length estimator. We first
start with an expression for the volume integrated flux, which can be written as
.. math::
:label: flux-integrated
V \phi = \int d\mathbf{r} \int dE \int d\mathbf{\Omega} \int dt \,
\psi(\mathbf{r}, \mathbf{\hat{\Omega}}, E, t)
where :math:`V` is the volume, :math:`\psi` is the angular flux,
:math:`\mathbf{r}` is the position of the particle, :math:`\mathbf{\hat{\Omega}}`
is the direction of the particle, :math:`E` is the energy of the particle, and
:math:`t` is the time. By noting that :math:`\psi(\mathbf{r},
\mathbf{\hat{\Omega}}, E, t) = v n(\mathbf{r}, \mathbf{\hat{\Omega}}, E, t)`
where :math:`n` is the angular neutron density, we can rewrite equation
:eq:`flux-integrated` as
.. math::
:label: flux-integrated-2
V \phi = \int d\mathbf{r} \int dE \int dt v \int d\mathbf{\Omega} \, n(\mathbf{r},
\mathbf{\hat{\Omega}}, E, t)).
Using the relations :math:`N(\mathbf{r}, E, t) = \int d\mathbf{\Omega}
n(\mathbf{r}, \mathbf{\hat{\Omega}}, E, t)` and :math:`d\ell = v \, dt` where
:math:`d\ell` is the differential unit of track length, we then obtain
.. math::
:label: track-length-integral
V \phi = \int d\mathbf{r} \int dE \int d\ell N(\mathbf{r}, E, t).
Equation :eq:`track-length-integral` indicates that we can use the length of a
particle's trajectory as an estimate for the flux, i.e. the track-length
estimator of the flux would be
.. math::
:label: track-length-flux
\phi = \frac{1}{W} \sum_{i \in T} w_i \ell_i
where :math:`T` is the set of all the particle's trajectories within the desired
volume and :math:`\ell_i` is the length of the :math:`i`-th trajectory. In the
same vein as equation :eq:`collision-estimator`, the track-length estimate of a
reaction rate is found by multiplying equation :eq:`track-length-flux` by a
macroscopic reaction cross section:
.. math::
:label: track-length-estimator
R_x = \frac{1}{W} \sum_{i \in T} w_i \ell_i \Sigma_x (E_i).
One important fact to take into consideration is that the use of a track-length
estimator precludes us from using any filter that requires knowledge of the
particle's state following a collision because by definition, it will not have
had a collision at every event. Thus, for tallies with outgoing-energy filters
(which require the post-collision energy), scattering change-in-angle filters,
or for tallies of scattering moments (which require the scattering cosine of
the change-in-angle), we must use an analog estimator.
.. TODO: Add description of surface current tallies
----------
Statistics
----------
As was discussed briefly in :ref:`methods_introduction`, any given result from a
Monte Carlo calculation, colloquially known as a "tally", represents an estimate
of the mean of some `random variable`_ of interest. This random variable
typically corresponds to some physical quantity like a reaction rate, a net
current across some surface, or the neutron flux in a region. Given that all
tallies are produced by a `stochastic process`_, there is an associated
uncertainty with each value reported. It is important to understand how the
uncertainty is calculated and what it tells us about our results. To that end,
we will introduce a number of theorems and results from statistics that should
shed some light on the interpretation of uncertainties.
Law of Large Numbers
--------------------
The `law of large numbers`_ is an important statistical result that tells us
that the average value of the result a large number of repeated experiments
should be close to the `expected value`_. Let :math:`X_1, X_2, \dots, X_n` be an
infinite sequence of `independent, identically-distributed random variables`_
with expected values :math:`E(X_1) = E(X_2) = \mu`. One form of the law of large
numbers states that the sample mean :math:`\bar{X_n} = \frac{X_1 + \dots +
X_n}{n}` `converges in probability`_ to the true mean, i.e. for all
:math:`\epsilon > 0`
.. math::
\lim\limits_{n\rightarrow\infty} P \left ( \left | \bar{X}_n - \mu \right |
\ge \epsilon \right ) = 0.
.. _central-limit-theorem:
Central Limit Theorem
---------------------
The `central limit theorem`_ (CLT) is perhaps the most well-known and ubiquitous
statistical theorem that has far-reaching implications across many
disciplines. The CLT is similar to the law of large numbers in that it tells us
the limiting behavior of the sample mean. Whereas the law of large numbers tells
us only that the value of the sample mean will converge to the expected value of
the distribution, the CLT says that the distribution of the sample mean will
converge to a `normal distribution`_. As we defined before, let :math:`X_1, X_2,
\dots, X_n` be an infinite sequence of independent, identically-distributed
random variables with expected values :math:`E(X_i) = \mu` and variances
:math:`\text{Var} (X_i) = \sigma^2 < \infty`. Note that we don't require that
these random variables take on any particular distribution -- they can be
normal, log-normal, Weibull, etc. The central limit theorem states that as
:math:`n \rightarrow \infty`, the random variable :math:`\sqrt{n} (\bar{X}_n -
\mu)` `converges in distribution`_ to the standard normal distribution:
.. math::
:label: central-limit-theorem
\sqrt{n} \left ( \frac{1}{n} \sum_{i=1}^n X_i - \mu \right ) \xrightarrow{d}
\mathcal{N} (0, \sigma^2)
Estimating Statistics of a Random Variable
------------------------------------------
Mean
++++
Given independent samples drawn from a random variable, the sample mean is
simply an estimate of the average value of the random variable. In a Monte Carlo
simulation, the random variable represents physical quantities that we want
tallied. If :math:`X` is the random variable with :math:`N` observations
:math:`x_1, x_2, \dots, x_N`, then an unbiased estimator for the population mean
is the sample mean, defined as
.. math::
:label: sample-mean
\bar{x} = \frac{1}{N} \sum_{i=1}^N x_i.
Variance
++++++++
The variance of a population indicates how spread out different members of the
population are. For a Monte Carlo simulation, the variance of a tally is a
measure of how precisely we know the tally value, with a lower variance
indicating a higher precision. There are a few different estimators for the
population variance. One of these is the second central moment of the
distribution also known as the biased sample variance:
.. math::
:label: biased-variance
s_N^2 = \frac{1}{N} \sum_{i=1}^N \left ( x_i - \bar{x} \right )^2 = \left (
\frac{1}{N} \sum_{i=1}^N x_i^2 \right ) - \bar{x}^2.
This estimator is biased because its expected value is actually not equal to the
population variance:
.. math::
:label: biased-variance-expectation
E[s_N^2] = \frac{N - 1}{N} \sigma^2
where :math:`\sigma^2` is the actual population variance. As a result, this
estimator should not be used in practice. Instead, one can use `Bessel's
correction`_ to come up with an unbiased sample variance estimator:
.. math::
:label: unbiased-variance
s^2 = \frac{1}{N - 1} \sum_{i=1}^N \left ( x_i - \bar{x} \right )^2 =
\frac{1}{N - 1} \left ( \sum_{i=1}^N x_i^2 - N\bar{x}^2 \right ).
This is the estimator normally used to calculate sample variance. The final form
in equation :eq:`unbiased-variance` is especially suitable for computation since
we do not need to store the values at every realization of the random variable
as the simulation proceeds. Instead, we can simply keep a running sum and sum of
squares of the values at each realization of the random variable and use that to
calculate the variance.
Variance of the Mean
++++++++++++++++++++
The previous sections discussed how to estimate the mean and variance of a
random variable using statistics on a finite sample. However, we are generally
not interested in the *variance of the random variable* itself; we are more
interested in the *variance of the estimated mean*. The sample mean is the
result of our simulation, and the variance of the sample mean will tell us how
confident we should be in our answers.
Fortunately, it is quite easy to estimate the variance of the mean if we are
able to estimate the variance of the random variable. We start with the
observation that if we have a series of uncorrelated random variables, we can
write the variance of their sum as the sum of their variances:
.. math::
:label: bienayme-formula
\text{Var} \left ( \sum_{i=1}^N X_i \right ) = \sum_{i=1}^N \text{Var} \left
( X_i \right )
This result is known as the Bienaymé formula. We can use this result to
determine a formula for the variance of the sample mean. Assuming that the
realizations of our random variable are again identical,
independently-distributed samples, then we have that
.. math::
:label: sample-variance-mean
\text{Var} \left ( \bar{X} \right ) = \text{Var} \left ( \frac{1}{N}
\sum_{i=1}^N X_i \right ) = \frac{1}{N^2} \sum_{i=1}^N \text{Var} \left (
X_i \right ) = \frac{1}{N^2} \left ( N\sigma^2 \right ) =
\frac{\sigma^2}{N}.
We can combine this result with equation :eq:`unbiased-variance` to come up with
an unbiased estimator for the variance of the sample mean:
.. math::
:label: sample-variance-mean-formula
s_{\bar{X}}^2 = \frac{1}{N - 1} \left ( \frac{1}{N} \sum_{i=1}^N x_i^2 -
\bar{x}^2 \right ).
At this point, an important distinction should be made between the estimator for
the variance of the population and the estimator for the variance of the
mean. As the number of realizations increases, the estimated variance of the
population based on equation :eq:`unbiased-variance` will tend to the true
population variance. On the other hand, the estimated variance of the mean will
tend to zero as the number of realizations increases. A practical interpretation
of this is that the longer you run a simulation, the better you know your
results. Therefore, by running a simulation long enough, it is possible to
reduce the stochastic uncertainty to arbitrarily low levels.
Confidence Intervals
++++++++++++++++++++
While the sample variance and standard deviation gives us some idea about the
variability of the estimate of the mean of whatever quantities we've tallied, it
does not help us interpret how confidence we should be in the results. To
quantity the reliability of our estimates, we can use `confidence intervals`_
based on the calculated sample variance.
A :math:`1-\alpha` confidence interval for a population parameter is defined as
such: if we repeat the same experiment many times and calculate the confidence
interval for each experiment, then :math:`1 - \alpha` percent of the calculated
intervals would encompass the true population parameter. Let :math:`x_1, x_2,
\dots, x_N` be samples from a set of independent, identically-distributed random
variables each with population mean :math:`\mu` and variance
:math:`\sigma^2`. The t-statistic is defined as
.. math::
:label: t-statistic
t = \frac{\bar{x} - \mu}{s/\sqrt{N}}
where :math:`\bar{x}` is the sample mean from equation :eq:`sample-mean` and
:math:`s` is the standard deviation based on equation
:eq:`unbiased-variance`. If the random variables :math:`X_i` are
normally-distributed, then the t-statistic has a `Student's t-distribution`_
with :math:`N-1` degrees of freedom. This implies that
.. math::
:label: t-probability
Pr \left ( -t_{1 - \alpha/2, N - 1} \le \frac{\bar{x} - \mu}{s/\sqrt{N}} \le
t_{1 - \alpha/2, N - 1} \right ) = 1 - \alpha
where :math:`t_{1-\alpha/2, N-1}` is the :math:`1 - \alpha/2` percentile of a
t-distribution with :math:`N-1` degrees of freedom. Thus, the :math:`1 - \alpha`
two sided confidence interval for the sample mean is
.. math::
:label: two-sided-ci
\bar{x} \pm t_{1 - \alpha/2, N-1} \frac{s}{\sqrt{N}}.
One should be cautioned that equation :eq:`two-sided-ci` only applies if the
*underlying random variables* are normally-distributed. In general, this may not
be true for a tally random variable --- the central limit theorem guarantees
only that the sample mean is normally distributed, not the underlying random
variable. If batching is used, then the underlying random variable, which would
then be the averages from each batch, will be normally distributed as long as
the conditions of the central limit theorem are met.
Let us now outline the method used to calculate the percentile of the Student's
t-distribution. For one or two degrees of freedom, the percentile can be written
analytically. For one degree of freedom, the t-distribution becomes a standard
`Cauchy distribution`_ whose cumulative distribution function is
.. math::
:label: cauchy-cdf
c(x) = \frac{1}{\pi} \arctan x + \frac{1}{2}.
Thus, inverting the cumulative distribution function, we find the :math:`x`
percentile of the standard Cauchy distribution to be
.. math::
:label: percentile-1
t_{x,1} = \tan \left ( \pi \left ( x - \frac{1}{2} \right ) \right ).
For two degrees of freedom, the cumulative distribution function is the
second-degree polynomial
.. math::
:label: t-2-polynomial
c(x) = \frac{1}{2} + \frac{x}{2\sqrt{x^2 + 2}}
Solving for :math:`x`, we find the :math:`x` percentile to be
.. math::
:label: percentile-2
t_{x,2} = \frac{2\sqrt{2} (x - 1/2)}{\sqrt{1 - 4 (x - 1/2)^2}}
For degrees of freedom greater than two, it is not possible to obtain an
analytical formula for the inverse of the cumulative distribution function. We
must resort to either numerically solving for the inverse or to an
approximation. Approximations for percentiles of the t-distribution have been
found with high levels of accuracy. OpenMC uses the `following approximation`_:
.. math::
:label: percentile-n
t_{x,n} = \sqrt{\frac{n}{n-2}} \left ( z_x + \frac{1}{4} \frac{z_x^3 -
3z_x}{n-2} + \frac{1}{96} \frac{5z_x^5 - 56z_x^3 + 75z_x}{(n-2)^2} +
\frac{1}{384} \frac{3z_x^7 - 81z_x^5 + 417z_x^3 - 315z_x}{(n-2)^3} \right )
where :math:`z_x` is the :math:`x` percentile of the standard normal
distribution. In order to determine an arbitrary percentile of the standard
normal distribution, we use an `unpublished rational approximation`_. After
using the rational approximation, one iteration of Newton's method is applied to
improve the estimate of the percentile.
.. only:: html
.. rubric:: References
.. _following approximation: https://doi.org/10.1080/03610918708812641
.. _Bessel's correction: https://en.wikipedia.org/wiki/Bessel's_correction
.. _random variable: https://en.wikipedia.org/wiki/Random_variable
.. _stochastic process: https://en.wikipedia.org/wiki/Stochastic_process
.. _independent, identically-distributed random variables: https://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables
.. _law of large numbers: https://en.wikipedia.org/wiki/Law_of_large_numbers
.. _expected value: https://en.wikipedia.org/wiki/Expected_value
.. _converges in probability: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability
.. _normal distribution: https://en.wikipedia.org/wiki/Normal_distribution
.. _converges in distribution: https://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
.. _confidence intervals: https://en.wikipedia.org/wiki/Confidence_interval
.. _Student's t-distribution: https://en.wikipedia.org/wiki/Student%27s_t-distribution
.. _Cauchy distribution: https://en.wikipedia.org/wiki/Cauchy_distribution
.. _unpublished rational approximation: https://web.archive.org/web/20150926021742/http://home.online.no/~pjacklam/notes/invnorm/
.. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf