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

Updated theory/methods section with changes from thesis. Closes #84 on github.

Browse source
This commit is contained in:
Paul Romano 2012-11-30 11:31:31 -05:00
parent 55b9b08893
commit 2c5eb197f8
10 changed files with 1867 additions and 894 deletions
Download patch file Download diff file Expand all files Collapse all files

792
docs/img/uniongrid.svg Normal file
View file

@ -0,0 +1,792 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:cc="http://creativecommons.org/ns#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns="http://www.w3.org/2000/svg"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
width="901.13098"
height="360.56549"
id="svg2"
version="1.1"
inkscape:version="0.48.3.1 r9886"
sodipodi:docname="New document 1">
<defs
id="defs4">
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend"
style="overflow:visible">
<path
id="path4037"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)"
inkscape:connector-curvature="0" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-2"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-2"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-26"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-3"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-4"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-7"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-9"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-4"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-1"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-0"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-0"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-5"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
<marker
inkscape:stockid="Arrow1Lend"
orient="auto"
refY="0"
refX="0"
id="Arrow1Lend-23"
style="overflow:visible">
<path
inkscape:connector-curvature="0"
id="path4037-9"
d="M 0,0 5,-5 -12.5,0 5,5 0,0 z"
style="fill-rule:evenodd;stroke:#000000;stroke-width:1pt"
transform="matrix(-0.8,0,0,-0.8,-10,0)" />
</marker>
</defs>
<sodipodi:namedview
id="base"
pagecolor="#ffffff"
bordercolor="#666666"
borderopacity="1.0"
inkscape:pageopacity="0.0"
inkscape:pageshadow="2"
inkscape:zoom="0.66342"
inkscape:cx="447.84509"
inkscape:cy="213.81917"
inkscape:document-units="px"
inkscape:current-layer="layer1"
showgrid="true"
fit-margin-top="0"
fit-margin-left="0"
fit-margin-right="0"
fit-margin-bottom="0"
inkscape:window-width="1366"
inkscape:window-height="712"
inkscape:window-x="0"
inkscape:window-y="27"
inkscape:window-maximized="1">
<inkscape:grid
type="xygrid"
id="grid2985"
empspacing="5"
visible="true"
enabled="true"
snapvisiblegridlinesonly="true"
originx="20.5655px"
originy="-299.43451px" />
</sodipodi:namedview>
<metadata
id="metadata7">
<rdf:RDF>
<cc:Work
rdf:about="">
<dc:format>image/svg+xml</dc:format>
<dc:type
rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
<dc:title></dc:title>
</cc:Work>
</rdf:RDF>
</metadata>
<g
inkscape:label="Layer 1"
inkscape:groupmode="layer"
id="layer1"
transform="translate(20.5655,-392.36218)">
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987"
width="60"
height="60"
x="160"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="170"
y="432.36218"
id="text2989"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991"
x="170"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3882">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993">1</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-3"
width="60"
height="60"
x="250"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="260"
y="432.36218"
id="text2989-4"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-9"
x="260"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3884">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-1">2</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-33"
width="60"
height="60"
x="340"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="350"
y="432.36218"
id="text2989-1"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8"
x="350"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3">3</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-6"
width="60"
height="60"
x="430"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="440"
y="432.36218"
id="text2989-48"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-5"
x="440"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3888">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-17">4</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-9"
width="60"
height="60"
x="520"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="530"
y="432.36218"
id="text2989-3"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-93"
x="530"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3890">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-4">5</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-2"
width="60"
height="60"
x="610"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="620"
y="432.36218"
id="text2989-0"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-0"
x="620"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3892">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-5">6</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-5"
width="60"
height="60"
x="700"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="710"
y="432.36218"
id="text2989-30"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-2"
x="710"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3894">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-8">7</tspan></tspan></text>
<rect
style="fill:#8fbdf4;fill-opacity:1;stroke:none"
id="rect2987-22"
width="60"
height="60"
x="790"
y="392.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="800"
y="432.36218"
id="text2989-00"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-09"
x="800"
y="432.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3896">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-0">8</tspan></tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="70"
y="412.36218"
id="text3898"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3900"
x="70"
y="412.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">Union</tspan><tspan
sodipodi:role="line"
x="70"
y="442.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3902">Energy Grid</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908"
width="60"
height="50"
x="160"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="180"
y="542.36218"
id="text3910"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912"
x="180"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">0</tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="70"
y="518.36218"
id="text3898-8"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3900-1"
x="70"
y="518.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">Nuclide</tspan><tspan
sodipodi:role="line"
x="70"
y="548.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3902-4">Pointers</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-2"
width="60"
height="50"
x="250"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="270"
y="542.36218"
id="text3910-9"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-8"
x="270"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">0</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-23"
width="60"
height="50"
x="340"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="360"
y="542.36218"
id="text3910-8"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-5"
x="360"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">1</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-9"
width="60"
height="50"
x="430"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="450"
y="542.36218"
id="text3910-7"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-9"
x="450"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">1</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-92"
width="60"
height="50"
x="520"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="540"
y="542.36218"
id="text3910-3"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-81"
x="540"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">1</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-0"
width="60"
height="50"
x="610"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="630"
y="542.36218"
id="text3910-0"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-1"
x="630"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">2</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-5"
width="60"
height="50"
x="700"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="720"
y="542.36218"
id="text3910-1"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-2"
x="720"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">3</tspan></text>
<rect
style="fill:#ffffff;fill-opacity:1;stroke:#000000;stroke-width:1.11800003;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
id="rect3908-235"
width="60"
height="50"
x="790"
y="502.36218"
ry="3.1880078" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="810"
y="542.36218"
id="text3910-80"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3912-0"
x="810"
y="542.36218"
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">3</tspan></text>
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 190,452.36218 0,50"
id="path4028"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 280,452.36218 0,50"
id="path4028-0"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 370,452.36218 0,50"
id="path4028-9"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 460,452.36218 0,50"
id="path4028-8"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 550,452.36218 0,50"
id="path4028-3"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 640,452.36218 0,50"
id="path4028-7"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 730,452.36218 0,50"
id="path4028-2"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1.29099441px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;marker-end:url(#Arrow1Lend)"
d="m 820,452.36218 0,50"
id="path4028-6"
inkscape:connector-curvature="0" />
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-9"
width="60"
height="60"
x="340"
y="582.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="350"
y="622.36218"
id="text2989-1-0"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-7"
x="350"
y="622.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-2">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-9">1</tspan></tspan></text>
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-1"
width="60"
height="60"
x="610"
y="582.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="620"
y="622.36218"
id="text2989-1-1"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-9"
x="620"
y="622.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-1">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-93">2</tspan></tspan></text>
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-0"
width="60"
height="60"
x="700"
y="582.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="710"
y="622.36218"
id="text2989-1-3"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-4"
x="710"
y="622.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-5">E</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-1">3</tspan></tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="70"
y="602.36218"
id="text3898-2"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3900-5"
x="70"
y="602.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">Nuclide</tspan><tspan
sodipodi:role="line"
x="70"
y="632.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3902-1">Energy Grid</tspan></text>
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-9-5"
width="60"
height="60"
x="340"
y="662.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="350"
y="702.36218"
id="text2989-1-0-9"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-7-8"
x="350"
y="702.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-2-3">σ</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-9-6">1</tspan></tspan></text>
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-9-5-0"
width="60"
height="60"
x="610"
y="662.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="620"
y="702.36218"
id="text2989-1-0-9-9"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-7-8-5"
x="620"
y="702.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-2-3-2">σ</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-9-6-2">2</tspan></tspan></text>
<rect
style="fill:#98f48f;fill-opacity:1;stroke:none"
id="rect2987-33-9-5-5"
width="60"
height="60"
x="700"
y="662.36218"
ry="2.5504062" />
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="710"
y="702.36218"
id="text2989-1-0-9-8"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan2991-8-7-8-7"
x="710"
y="702.36218"><tspan
style="font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3886-2-3-9">σ</tspan><tspan
style="font-size:65.00091553%;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:start;line-height:125%;writing-mode:lr-tb;text-anchor:start;baseline-shift:sub;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan2993-3-9-6-9">3</tspan></tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;text-align:center;line-height:125%;letter-spacing:0px;word-spacing:0px;text-anchor:middle;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="70"
y="682.36218"
id="text3898-2-6"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3900-5-6"
x="70"
y="682.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter">Nuclide</tspan><tspan
sodipodi:role="line"
x="70"
y="712.36218"
style="font-size:24px;font-style:normal;font-variant:normal;font-weight:normal;font-stretch:normal;text-align:center;text-anchor:middle;font-family:Bitstream Charter;-inkscape-font-specification:Bitstream Charter"
id="tspan3902-1-1">Cross Sections</tspan></text>
<text
xml:space="preserve"
style="font-size:40px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="251.72591"
y="769.73578"
id="text4762"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan4764"
x="251.72591"
y="769.73578"></tspan></text>
<rect
style="fill:none;stroke:#000000;stroke-width:2;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:4,2;stroke-dashoffset:0"
id="rect4766"
width="900"
height="280"
x="-20"
y="472.36218"
ry="2.5504062" />
</g>
</svg>
Side by side

After

Width:  |  Height:  |  Size: 35 KiB

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

@ -0,0 +1,94 @@
.. _methods_cross_sections:
============================
Cross Section Representation
============================
The data governing the interaction of neutrons with various nuclei are
represented using 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. The use of a standard cross section format allows for a
direct comparison of OpenMC with other codes since the same cross section
libraries can be used.
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.
Unionized Energy Grid
---------------------
The most naïve method to reduce the number of energy grid searches is to
construct a new energy grid that consists of the union of the energy points of
each nuclide and use this energy grid for all nuclides. This method is
computationally very efficient as it only requires one energy grid search at
each collision as well as one interpolation between cross section values since
the interpolation factor can be used for all nuclides. However, it requires
redundant storage of cross section values at points which were added to each
nuclide grid. This additional burden on memory storage can become quite
prohibitive. To lessen that burden, the unionized energy grid can be thinned
with cross sections reconstructed on the thinned energy grid. This method is
currently used by default in the Serpent Monte Carlo code.
Unionized Energy Grid with Nuclide Pointers
-------------------------------------------
While having a unionized grid that is used for all nuclides allows for very fast
lookup of cross sections, the burden on memory is in many circumstances
unacceptable. The OpenMC Monte Carlo code utilizes a method that allows for a
single energy grid search to be performed at every collision while avoiding the
redundant storage of cross section values. Instead of using the unionized grid
for every nuclide, the original energy grid of each nuclide is kept and a list
of pointers (of the same length as the unionized energy grid) is constructed for
each nuclide that gives the corresponding grid index on the nuclide grid for a
given grid index on the unionized grid. One must still interpolate on cross
section values for each nuclide since the interpolation factors will generally
be different. The figure below illustrates this method. All values within the
dashed box would need to be stored on a per-nuclide basis, and the union grid
would need to be stored once. This method is also referred to as *double
indexing* and is available as an option in Serpent (see paper by Leppanen_).
.. figure:: ../../img/uniongrid.svg
:width: 600px
:align: center
:figclass: align-center
Mapping of union energy grid to nuclide energy grid through pointers.
.. _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: http://dx.doi.org/10.1016/j.anucene.2009.03.019

34
docs/source/methods/criticality.rst → docs/source/methods/eigenvalue.rst
View file

@ -1,20 +1,20 @@
.. _methods_criticality:
.. _methods_eigenvalue:
========================
Criticality Calculations
========================
=======================
Eigenvalue Calculations
=======================
A criticality calculation is a transport simulation wherein the source of
neutrons includes a fissionable material. Some common criticality calculations
include the simulation of nuclear reactors, spent fuel pools, nuclear weapons,
and other fissile systems. The term criticality calculation is also synonymous
with the term eigenvalue calculation. The reason for this 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. Criticality simulations using Monte Carlo methods are becoming
increasingly common with the advent of high-performance computing.
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 criticality
This section will explore the theory behind and implementation of eigenvalue
calculations in a Monte Carlo code.
.. _method-successive-generations:
@ -23,7 +23,7 @@ calculations in a Monte Carlo code.
Method of Successive Generations
--------------------------------
The method used to converge on the fission source distribution in a criticality
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,
@ -48,7 +48,7 @@ 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 achiable parallel scaling. This topic is discussed at length
large impact on the achievable parallel scaling. This topic is discussed at length
in :ref:`fission-bank-algorithms`.
-------------------------
@ -74,7 +74,7 @@ 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 criticality calculations is to use
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.

224
docs/source/methods/geometry.rst
View file

@ -11,12 +11,13 @@ 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, intersection, or difference 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
*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. The region for which :math:`f(x,y,z)
< 0` can be called the negative half-space (or simply the "negative side") and
the region for which :math:`f(x,y,z) > 0` can be called the positive half-space.
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
@ -41,28 +42,77 @@ One can confirm that any point inside this sphere will correspond to
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. For example, to refer to the negative half-space of a sphere (the volume
inside the sphere) with unique ID 35, the reference would be -35. These
references to half-spaces are used in created regions in space of homogeneous
material, known as "cells".
sign. The following illustration shows an example of an ellipse with unique ID 1
dividing space into two half-spaces.
.. figure:: ../../img/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, known as *cells*. While some codes allow regions
to be defined by intersections, unions, and differences or half-spaces, OpenMC
is currently limited to cells defined only as intersections of
half-spaces. Thus, the specification of the cell must include a list of
half-space references whose intersection defines the region. The region is then
assigned a material defined elsewhere. The following illustration shows an
example of a cell defined as the intersection of an ellipse and two planes.
.. figure:: ../../img/union.svg
:align: center
:figclass: align-center
In OpenMC, any second-order surface of the form
The shaded region represents a cell bounded by three surfaces.
.. math::
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.
f(x,y,z) = Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
.. table:: Surface types available in OpenMC.
can be modeled in OpenMC. For example, the equation for a sphere centered at
:math:`(\bar{x},\bar{y},\bar{z})` and of radius :math:`R` can be written as
+----------------------+------------+------------------------------+-------------------------+
| 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:`x - x_0 = 0` | :math:`y_0` |
| to :math:`y`-axis | | | |
+----------------------+------------+------------------------------+-------------------------+
| Plane perpendicular | z-plane | :math:`x - x_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` |
+----------------------+------------+------------------------------+-------------------------+
.. _universes:
@ -75,13 +125,12 @@ 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 closed volume, or cell, in OpenMC can either be filled with a normal
material or with a universe. If the cell is filled with a univese, 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.
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
@ -91,7 +140,7 @@ 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
specified universe is automatically part of the *base universe* whose
coordinates are just the normal coordinates in Euclidean space.
Lattices
@ -103,7 +152,7 @@ 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
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
@ -127,17 +176,19 @@ 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.
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,y,z)` traveling in the
direction :math:`u,v,w`. To find the distance :math:`d` to a surface
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 + du, y + dv, z + dw) = 0
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
@ -147,6 +198,15 @@ 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
@ -165,10 +225,6 @@ 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.
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.
Plane Perpendicular to an Axis
------------------------------
@ -303,6 +359,8 @@ 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.
.. TODO: Need to add derivation for x-cone, y-cone, and z-cone.
.. _find-cell:
----------------------------
@ -314,7 +372,7 @@ 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 do a loop over each cell within
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
@ -331,30 +389,27 @@ is found that contains the specified point.
Determining if a Coordinate is in a Cell
----------------------------------------
One aspect of being able to determine what cell a particle is in is determining
if a particle's coordinates lie within a given cell. The current geometry
implementation in OpenMC limits all cells to being simple cells, i.e. they are
defined only with intersection of half-spaces and not unions, differences,
etc. This makes the job of determining if a point is in a cell quite simple.
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
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.
is inside the cell. Note that this algorithm works only for *simple cells*
defined as intersections of half-spaces.
Let us illustrate this idea with a concept. Let's say we have a cell defined as
It may help to illustrate this algorithm using a simple example. Let's say we
have a cell defined as
.. code-block:: xml
<cell id="1" surfaces="-1 2 -3" />
<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
@ -368,9 +423,9 @@ satisfy the following equations
x - (-3) > 0 \\
x - 2 < 0
So in order to determine if a point is inside the cell, we would plug its
coordinates into equation :eq:`cell-contains-example` and if the inequalities
are satisfied, than the point is indeed inside the cell.
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
@ -390,9 +445,9 @@ 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 except 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.
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:
@ -401,15 +456,15 @@ 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 contain the
surface in their specification 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.
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:
@ -432,9 +487,9 @@ 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 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
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
@ -442,10 +497,10 @@ equation :eq:`reflection-v` to
\mathbf{\Omega'} = \mathbf{\Omega} - 2 (\mathbf{\Omega} \cdot
\hat{\mathbf{n}}) \hat{\mathbf{n}}
The direction of the surface normal will be the gradient to 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
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
@ -471,8 +526,8 @@ series of equations:
w' = w - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2}
\frac{\partial f}{\partial z}
We can now use this form to develop rules for how to transform a particle's
direction for different types of surfaces.
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
------------------------------
@ -523,7 +578,8 @@ 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
.. 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}
@ -532,13 +588,15 @@ A cylinder parallel to, for example, the x-axis has the form :math:`f(x,y,z) =
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
.. 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
.. math::
:label: reflection-cylinder-constant
\frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} =
\frac{\bar{y}v + \bar{z}w}{R^2}
@ -548,7 +606,8 @@ Substituting equations :eq:`reflection-cylinder-constant` and
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
.. math::
:label: reflection-cylinder
v' = v - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{y}}{R^2} \\
@ -561,7 +620,8 @@ 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
.. 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} \\
@ -570,13 +630,15 @@ The surface equation for a sphere has the form :math:`f(x,y,z) = (x - x_0)^2 +
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
.. 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
.. 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}
@ -585,14 +647,16 @@ 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
.. 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} \\
w' = w - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{z} }{R^2}
.. TODO: Add in derivation for cone surfaces.
.. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry
.. _surfaces: http://en.wikipedia.org/wiki/Surface

5
docs/source/methods/index.rst
View file

@ -9,9 +9,10 @@ Theory and Methodology
:maxdepth: 3
introduction
criticality
statistics
geometry
cross_sections
random_numbers
physics
tallies
eigenvalue
parallelization

20
docs/source/methods/introduction.rst
View file

@ -45,13 +45,13 @@ following steps:
- Initialize the pseudorandom number generator.
- Read ACE format cross-sections specified in the problem.
- Read ACE format cross sections 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 fixed source problem, source sites are sampled from the specified
source. In a criticality problem, source sites are sampled from some initial
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.
@ -64,15 +64,15 @@ proceed. The life of a single particle will proceed as follows:
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
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.
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
cross section is :math:`\Sigma_t`, this can be shown to be
.. math::
@ -88,7 +88,7 @@ proceed. The life of a single particle will proceed as follows:
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
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::
@ -97,7 +97,7 @@ proceed. The life of a single particle will proceed as follows:
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
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
@ -106,8 +106,8 @@ proceed. The life of a single particle will proceed as follows:
P(x) = \frac{\sigma_x}{\sigma_t}.
9. If the sampled reaction is elastic or inelastic scattering, the outgoing
energy and angle is sampled from the appropriate distribution. If the
reaction is (n,xn), it's also treated as scattering and the weight of the
energy and angle is sampled from the appropriate distribution. Reactions
of type :math:`(n,xn)` are treated as scattering and the weight of the
particle is increased by the multiplicity of the reaction. The 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
@ -121,7 +121,7 @@ be performed before the run is finished. This include the following:
- All tallies and other results are written to disk.
- If requested, a source file is written to disk
- If requested, a source file is written to disk.
- All allocatable arrays are deallocated.

815
docs/source/methods/physics.rst
View file

File diff suppressed because it is too large Load diff

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

@ -0,0 +1,74 @@
.. _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
[LEcuyer]_).
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 :eq:`lcg-skipahead` has the same general form as \eqref{eq:lcg}, so
the idea is to determine the new multiplicative and additive constants in
:math:`O(\log_2 N)` operations.
----------
References
----------
.. [LEcuyer] P. L’Ecuyer, "Tables of Linear Congruential Generators of
Different Sizes and Good Lattice Structures," *Math. Comput.*, **68**, 249
(1999).
.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl_rn_arb-strides_1994.pdf
.. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator

369
docs/source/methods/statistics.rst
View file

@ -1,369 +0,0 @@
.. _methods_statistics:
==========
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 approximation from
[George]_:
.. 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.
------------------------
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 a properties that
make it appear to be random, e.g. being uniformly distributed and not exhibiting
correlation between successive terms. Since the numbers are not truly "random"
in the strict sense, they are typically referred to as pseudo-random numbers,
and the techniques used to generate them are pseudo-random number generators
(PRNGs). Numbers sampled on the unit interval can then be used transformed for
the purpose of sampling other continuous or discrete probability distributions.
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 some random number seed :math:`\xi_0` and a sequence
of random numbers is 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 said 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 binary 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}`.
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
just 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 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 form as equation :eq:`lcg`
so the idea is to determine the new multiplicative and additive constants in
:math:`O(\log N)` operations.
----------
References
----------
.. [George] E. E. Olusegun George and Meenakshi Sivaram, "A modification of the
Fisher-Cornish approximation for the student t percentiles," Communication
in Statistics - Simulation and Computation, 16 (4), pp. 1123-1132 (1987).
.. _linear congruential generator: http://en.wikipedia.org/wiki/Linear_congruential_generator
.. _Brown: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/anl_rn_arb-strides_1994.pdf
.. _Bessel's correction: http://en.wikipedia.org/wiki/Bessel's_correction
.. _random variable: http://en.wikipedia.org/wiki/Random_variable
.. _stochastic process: http://en.wikipedia.org/wiki/Stochastic_process
.. _independent, identically-distributed random variables: http://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables
.. _law of large numbers: http://en.wikipedia.org/wiki/Law_of_large_numbers
.. _expected value: http://en.wikipedia.org/wiki/Expected_value
.. _converges in probability: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability
.. _normal distribution: http://en.wikipedia.org/wiki/Normal_distribution
.. _converges in distribution: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
.. _confidence intervals: http://en.wikipedia.org/wiki/Confidence_interval
.. _Student's t-distribution: http://en.wikipedia.org/wiki/Student%27s_t-distribution
.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution
.. _unpublished rational approximation: http://home.online.no/~pjacklam/notes/invnorm/

334
docs/source/methods/tallies.rst
View file

@ -29,16 +29,16 @@ equation :eq:`tally-integral`). For example, if the desired tally was 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, (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, and
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, and
post-collision energy.
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
sections can subsequently be used in deterministic solvers such as coarse mesh
finite difference (CMFD) diffusion.
------------------------------
@ -87,7 +87,9 @@ 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.
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
-------------------
@ -145,7 +147,7 @@ start with an expression for the volume integrated flux, which can be written as
: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).
\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}}`
@ -159,7 +161,7 @@ where :math:`n` is the angular neutron density, we can rewrite equation
: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))
\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
@ -168,7 +170,7 @@ n(\mathbf{r}, \mathbf{\hat{\Omega}}, E, t)` and :math:`d\ell = v \, dt` where
.. math::
:label: track-length-integral
V \phi = \int d\mathbf{r} \int dE \int d\ell N(\mathbf{r}, E, t)
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
@ -188,7 +190,7 @@ 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)
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
@ -197,8 +199,314 @@ had a collision at every event. Thus, for tallies with outgoing-energy filters
(which require the post-collision energy) or for tallies of scattering moments
(which require the scattering cosine), we must use an analog estimator.
---------------
Surface Current
---------------
.. 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 approximation from
[George]_:
.. 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.
----------
References
----------
.. [George] E. E. Olusegun George and Meenakshi Sivaram, "A modification of the
Fisher-Cornish approximation for the student t percentiles," Communication
in Statistics - Simulation and Computation, 16 (4), pp. 1123-1132 (1987).
.. _Bessel's correction: http://en.wikipedia.org/wiki/Bessel's_correction
.. _random variable: http://en.wikipedia.org/wiki/Random_variable
.. _stochastic process: http://en.wikipedia.org/wiki/Stochastic_process
.. _independent, identically-distributed random variables: http://en.wikipedia.org/wiki/Independent_and_identically_distributed_random_variables
.. _law of large numbers: http://en.wikipedia.org/wiki/Law_of_large_numbers
.. _expected value: http://en.wikipedia.org/wiki/Expected_value
.. _converges in probability: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_probability
.. _normal distribution: http://en.wikipedia.org/wiki/Normal_distribution
.. _converges in distribution: http://en.wikipedia.org/wiki/Convergence_of_random_variables#Convergence_in_distribution
.. _confidence intervals: http://en.wikipedia.org/wiki/Confidence_interval
.. _Student's t-distribution: http://en.wikipedia.org/wiki/Student%27s_t-distribution
.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution
.. _unpublished rational approximation: http://home.online.no/~pjacklam/notes/invnorm/
.. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf